Optimal multicompartmental sampling designs for parameter estimation: Practical aspects of the identification problem

Mathematics and Computers in Simulation XXIV (1982) 525-530 North-Holland Publishing Company

525

OPTIMAL M U L T I C O M P A R T M E N T A L SAMPLING DESIGNS FOR PARAMETER ESTIMATION: PRACTICAL ASPECTS OF THE IDENTIFICATION PROBLEM Elliot M. L A N D A W Department of Biomathematics, UCLA School of Medicine, Los Angeles, CA 90024, U.S.A.

A computer program is described for the determination of optimal multicompartmental sampling designs as a handle on numerical identifiability. A design measure approach is used which automatically chooses the number of sample times and sites of sampling. A two-input four-pool output example is presented.

i.

INTRODUCTION

The experimental identification of parameters of a dynamical multicompartmental model often requires multi-pool sampling and/or multiple known inputs (1,2,3). In this context, identifiability commonly refers to whether or not the parameters of the model can be (uniquely) determined assuming that the model structure is known and that a sufficient amount of data free from measurement (output) error can be collected from specified compartments. A number of other identifiability concepts do take into account the presence of measurement noise (4), but any particular planned experiment is still considered either identifiable or unidentifiable. In practice the existence of "noisy" data may severely limit the precision with which parameters can be estimated, even if they are identifiable in principle. This has been noted as a problem of "determinacy" or numerical identifiability (3,5). If the number K of timed samples that can be collected is specified in advance, it may be possible to choose the sample times judiciously in order to optimize parameter estimation precision (minimize estimation variability). Currently, several computer programs can find such exact K-point designs when the optimization criterion is minimization of the determinant (6,7,8) or weighted trace (9) of the asymptotic covarianee matrix of parameter estimates. These programs have generally been used for the single-input singleoutput (SISO) case, and in principle are easily extended to multi-pool sampling. However, they all require the user to specify how many samples are to be collected and which compartments and how m~ny samples per compartment. This paper will describe the extension of a previous program (i0) using a design measure approach (ii), which automatically generates the times, choice of compartments and number of samples per ompartment for the multi-pool sampling problem. 2.

NOTATION AND DEFINITIONS

2.1 Error Model

partment. Here "compartment" is used in a general sense. If, for example, there are m sites numbered l,...,m that can be sampled in each of j distinct experiments (e.g., different inputs) one can let N=jm and the i+(j-l)m'th component of ~ denote a sample from site i under experimental condition j. We assume that z(t) = z(t,~) + e(t)

(I)

where ~ is deterministic (typically some function of the state equations with known input and initial conditions) and e is a P-vector of identifiable parameters, e is a zero mean random vector representing measurement error whose covariance matrix is known up to a constant factor:

I(

0

if ti#t j

-

COV(~(ti)'~(tJ))-

(2)

°~(ti)

O

1 if ti=t j

L\ O"o (ti)/ This simplifying assumption of zero covariance across time and compartments is commonly held for many multicompartmental estimation problems in the biomedical sciences. Note the possible time dependence of ~ . In the example o~ will be made proportional to a constant plus a power function of Yi" 2.2 Sample or Design Space Let s denote a "design point," namely an ordered pair whose first component is an integer valued indicator of the compartment to be sampled and second component is the time of the sample. For example, a sample taken at 12.5 minutes in compartment 2 would be represented by s=(2,12.5). We let SC{ {l,2,...,N} X ~ } be the allowable space of design points. 2.3 Design Measure (11,12) A design measure ~(s) is a probability measure on the design space S which describes the proportion of independent samples taken at design point s. We need only consider discrete design measures:

Let z(t) be an N-vector whose i'th component is the sampled output at time t of the i'th com-

0378-4754/82/0000-0000/$02.75 © 1982 IMACS/North-Holland

526

E.M. Landaw / Optimal multicompartmentalsampling designs for parameter estimation

] Pi ~(s)

=

]

0

if s=s i

(3)

equivalent to either of the following two statements:

otherwise min d(6) = d(6*)

with Zp = i. The points s. are the points of i i support of ~ and the p. are their respective "masses" or proportionS. This approach allows the specification of a design without reference to the total number of sample points. However, any realizable design of K total samples will generally require that pjK (the number of independent replicates at si[ be rounded off to the nearest integer.

d(6*) = P Furthermore,

We let q(s,~) be the expected value and o2(s) the variance of the observation at s. Note that ~(s,%)=y.(t,@) and o2(s)=o~(t) if s=(j,t). Let f(s,~) b~ the column sensilivity P-vector

Step 2:

f(s,O) = 3~138 the normalized

(4) Step 4:

The normalized weighted response surface variance is d(s,6) = w(s)f t M-I(6)£

(6)

and we denote the maximum of d(s,6) over s e S as d(6).

p

(8)

pick a nonsingular nominal design ~0; i=0 _find s.l ~ S such that d(si,6i ) = d(6.); compute

=l(d(6.)-P)/P(d(6.)-i ) Step 3:

information matrix

M(6) = Z PiW(si)f(si,0)ft(si,0 ) (5) i where w(s) = I/o2(s). Note that if K is the total number of samples then KM(~) is the usual P X P Fisher information matrix (13), and if M is nonsingular M-I/K is the Cramer-Rao lower bound. For many parameter estimation procedures such as weighted nonlinear least squares the latter is an estimate of the covariance matrix for the parameter estimates.

=

This is the key to the basic Fedorov algorithm (11,12,14) used in this program for constructing D-optimal design measures: Step i:

and define

(7b)

at the points of support s* of 6"

d(s*,~*)

Information Matrix

2.4 Normalized

(7a)

construc~ 6i+i by a~ding s i to the points of support with mass ~ and multiply the remaining masses of 6 i by the factor i-~ i ÷ i+l; go to Step 2

It can be proved that the sequence ( i ~ * . Also, if the D-efficiency of a design ~ is computed as: e(6) = (det(M(6))/det(M(6*)))i/p

(9)

then (15)

P/d(6) ! e(6) ! 1

(10)

Thus a sim_ple criterion for convergence can be based on d(6i). Note that the term "efficiency" is used to indicate effective sample size compared to the D-optimal design. For example, a design measure with D-efficiency of 0.75 applied to a total sample size i00 will yield an information matrix whose determinant is the same as that for the D-optimal design measure applied to a total sample size of only 75.

2.5 Design Criterion We assume that 6 is (locally) information matrix identifiable (4) so that M is invertible for some design. The goal is then to specify a design 6 which minimizes in some sense M-I(6). One commonly chosen criterion is the minimization of det(M-l($)) called D-optimality. It can be proved that D-optimal designs are invariant under any rescaling of the parameters or any (locally) invertible reparameterization of the system (10,11,14). ~In addition a D-optimal design can be interpreted as minimizing the volume of the parameter vector confidence ellipsoid. The development can be extended to a wider class of criteria, e.g., minimizing trace (BM -I) for some fixed choice of nonnegative definite constant matrix B (10,11,12), but only D-optimality will be discussed. 3.

Step 2 is the time limiting task in the Fedorov algorithm but in the present context simply requires the equivalent of a one-dimensional search for maxima. 3.2 Atwood-St.

John Modification

(14)

At the beginning of each iteration loop the design point already in 6 i with the smallest value of (d(s,$i)-p) is found. If this latter value is negative a trial is made to remove the point from the design and transfer all of its mass to those remaining points sj in the design for which d(s~, 6i)-P is positive, and to distribute the miss proportional to d(s ,6 )-P. If this trial improves the design itJislaccept ed. This modification is very powerful in removing poor design points from the starting design.

COMPUTER PROGRAM 3.3 Cluster Identification

3.1 Equivalence Theorem and Fedorov Algorithm The Kiefer-Wolfowitz theorem (ii) states that the D-optimality of a design measure 6" is

The number of points of support in an optimal design need be no larger than P(P+I)/2 and often is equal to or slightly greater than

E.M. Landaw / Optimal multicompartmentalsampling designs for parameter estimation P (8,10,11). Since the Fedorov algorithm may add points which are not precisely equal to previous points in the design, the user is asked to check every k'th iteration (user chosen) on the number of design points. A user selected "resolution factor" determines how far apart two times in the design must be before they are considered distinct. A set of indistinguishable times within a compartment then forms a cluster, and each cluster is replaced by a single design point whose time is the (weighted) average of the times in the cluster and whose mass is the sum of the cluster masses.

527

p9

3.4 User Requirements Since the D-optimal design depends on the design space, error variance and nominals for some or all of the parameters (8,10) we have tended to view the development of such designs as an exploratory process best performed in an interactive environment. Unless one of the elementary SISO linear models is being used, the user must supply a subroutine for the output sensitivities and expected values. The design space, nominals, number of iterations between cluster searches, resolution factor, etc. are all supplied interactively. One programmed choice for the error variance of each compartment is the following parameterized structure: o2(s)

= a + bn ¥

,

(ii)

where a, b and y are indicated for each feasible compartment. The user is given the option to monitor the Fedorov algorithm, including the location of local maxima in the search of Step 2. This is particularly important in identifying the number of points of support, choosing the resolution factor and deciding when to terminate the program. The program is written in FORTRAN and is currently being incorporated as an option in the OSS package (8). 4.

EXAMPLE

4.1 Ketone Body Model The compartmental model shown in figure 1 has been used by Cobelli et. al. (16) to model ketone body metabolism and distribution in man. Pools 1 and 2 represent blood acetoacetate (AcAc) and B hydroxybutyrate (BOHB), and pools 3 and 4 represent equilibrating liver and extrahepatic tissue pools. Parameters pl,...,p9 are first order rate constants (min -I) and pl0"and pll are distribution volumes (ml). This model is uniquely identifiable if two su~arate experiments are performed (expt. A: 108 dpm of labeled AcAc injected as a bolus into pool i; expt. B: 7XI07 dpm labeled BOHB into pool 2) and sampling is confined to concentrations of label in pools 1 and 2.

Figure i:

Model for example

Although pools 3 and 4 are relatively inaccessible in man, it is instructive to include these sites as potential output ports, z is then 8dimensional with components IA,...,4A,IB,,,,.4B, where the number indicates the pool and the letter indicates the experiment. It is assumed that pools 1 and 2 are sampled as radioactivity concentrations whose expected values are the amount in the pool divided by the distribution volume and that the amounts in pools 3 and 4 are sampled directly (e.g., through scanning techniques). The error variance structure used by Cobelli et. al.(16) is that for Poisson-type error and is modelled by a=0 and y=l in equation (ii). Nominal values for the ii parameters are .768, .119, .025, .512, .260, .064, .097, .057, .180, 5530 and 5530. Sensitivities and expected values were computed by the standard method of appending sensitivity equations to the state differential equations and numeri c a l l y solving them once through a 4th order Runge-Kutta routine, i00 interpolation points were used for storage on a geometric scale from 1 to i00 minutes. 4.2 Comparison to Conventional Sampling Schedule A conventional schedule draws samples at times 1,2,3,4,5,6,8,10,12.5,15,20,25,30,40,50,60,70, 80 and i00 minutes in each of compartments IA, 2A,IB and 2B. Assuming b=l in error variance equation (Ii) for these 4 compartments one finds parameter estimate coefficients of variation (c.v.) for this design ranging from 2.3% to 11.3% (average for all ii parameters 5.1%). The D-optimal design measure was computed after confining the design space to be the time interval [i,i00] in these four compartments only (option chosen to make 3A,4A,3B and 4B inaccessible). This design is summarized in the first column of Table 1 labeled "no deep sampling." Note that the design is specified by giving the points of support and masses for each compartment. For example, in compartment 2A samples are drawn at times 2.95 and 18.8 minutes with masses .091 and .079 respectively. The D-efficiency of the conventional design compared to this design is 61.3%, and the D-optimal design applied to a total sample size of 76 provides parameter c.v.'s ranging from i.i to 9.5%

528

E.M. Landaw / Optimal multicompartmentalsampling destgns for parameter estimation

D-optimal Designs no "deep" sampling

time(mass)

"deep" sampling allowed with following b values (see section 4.3) b:49770

b=5530

b=l

IA

1.00(.091) 2.68(.091) 9.4 (.091) 29.2 (.041)

1.00(.090) 2.7 (.089) 5.8 (.009)

1.00(.091) 5.59(.083)

2A

2.95(.091) 18.8 (.079)

3.03(.089) 17.8 (.010)

3.1 (.086)

2.44(.087) 37.8 (.015)

1.9 (.089) 33.8 (.024)

1.47(.090) ll.4 (.O78)

10.7 (.076)

1.00(.090) l l . O (.090)

1.00(.090) 3.6 (.091) 12.1 (.083) 29.6 (.045)

1.00(.091)

Pool (1,2,3,4)

3A

and Experiment (A,B)

4A

IB

2.17(.091) 21.2 (.070)

2.4 (.081) 16.5 (.072)

13.8 (.029)

1.00(.090) 2.89(.086) 13.4 (.074)

1.00(.091) 13.1 (.079)

1.00(.091)

2B

1.00(.091) 2.87(.091) 11.7 (.082) 56.3 (.091)

2.3 (.055) 39.0 (.077)

1.46(.090) 40.4 (.071)

1.00(.086) 6.0 (.077) 49.4 (.091)

1.00(.088)

1.00(.087)

2.450X1023

3.899X1052

3B

4B det(M) d(~)

1.923X1017 11.04 TABLE i.

1.294X1019 11.04

11.03

11.03

D-optlmal Design Specifications

(average 4.2%). Rounded-off designs (e.g., rounding 76Xmass to the nearest integer) were quite similar with D-efficiencies better than 99%. Clearly the conventional design is not too bad for reaching a good degree of "determinancy," but the efficiency~shows that in one sense nearly 40% of the samples are wasted. The Doptimal design indicates that better parameter estimate precision can be obtained with more replicates taken at fewer design points. 4.3 Effect of Error Variance on Multi-pool Sampling For the second demonstration the design times were opened up to [i,i00] in all 8 compartments. The error variance (eq.(ll)) used a=0 and y=l in all compartments plus b=l in compartments IA,2A,IB,2B and b=5.53X107, 9X5530, 5530, or i in the remaining deep pools. Doptimal designs were computed starting with a

nominal design of samples at 1,5,10 and 30 minutes in each of the 8 compartments. The final designs are listed in the last 3 columns of Table l. It should be noted that for b=5.53X107 the Atwood-St. John modification quickly tossed out all "deep" samples, resulting in a design identical to the one in column one. Decreasing values of b led to more samples in the deeper pools, and at the smallest value of b only 2 points of support remain in the conventionally sampled ports. Heuristically, this indicates that for the latter degree of observational precision, the microparameter information is best obtained from the deep compartments, although single samples in IA and 2B at the earliest time after the bolus are needed to estimate pl0 and pll. 5.

DISCUSSION

The computing of optimal designs appears to be

E.M. Landaw / Optimal multicompartmental sampling designs for parameter estimation a useful tool for exploring the numerical identification problem for multi-pool sampling (7,8, i0). A major advantage of the design measure approach is that the user is not required to specify how many samples are to be collected nor w h i c h of the feasible compartments must be sampled. This approach is also applicable to a wide variety of optimality criteria based on the information matrix or subsets (11,12). The program that has been developed is extremely robust against poor choices for the initiating design and will generally converge quickly to designs that are better than 90 to 95% D-efficient. However, it is very slow in achieving 3 digit accuracy of the design specification. We are currently testing a two step procedure for the OSS package (8) in which the design measure program is used to reach a neighborhood of the optimal design and after rounding-off the design measure one of the exact design programs is used for the final few iterations. It is clear that this program alone cannot address traditional algebraic parameter identifiability questions. Optimal designs require knowledge of parameter nominals and variance structure, items not of concern in the purely algebraic problem. In addition, invertibility of the information matrix does not guarantee unique identifiability. Nevertheless, it is interesting to observe in the previous section that in one case the algorithm was able to toss out extraneous compartments from consideration and achieve what appears to be a minimally identifiable design. In contrast to exact K-point designs, the design measure approach may produce values for the masses Pi which do not translate into whole number replicates at the sample times. However, rounded-off designs (see section 2.3) are usually very near-optimal (18); in the example, rounded-off designs had D-efficiencies better than 99%. A more serious concern is that optimal designs (both exact and those based on design measures) tend to be schedules with independent replicates drawn at a minimal number of distinct sample times. Such designs are poor for model verification and testing assumptions about the error variance structure. One compro ~ mise for practical work is to use a sub-optimal design in which replicates in the design measure are spread out (or "split") over a range of times centered at the optimal points of support (10). Increasing levels of splitting result in a "menu" of designs with decreasing efficieneies but increasing usefulness for model testing. These sRlit designs also appear to have some robustness against nominal parameter misspecification and may complement other approaches to design robustness (7,19). Finally, the design measure approach appears to be of value in attacking some theoretical questions about optimal designs in compartmental models. Several theorems about the effect of the design space boundary and error variance on the points of support of a design have already

529

been presented for SISO linear models (I0). In this situation it has also been conjectured that the number of points of support for D-optimal designs generally is equal to the number of parameters if the error variance is sufficiently "smooth" (8,10,17). The appearance of computed designs with 12 or more points of support in Table 1 suggests that this conjecture may not hold when extended to the multi-pool sampling problem. Although these are numerical results, their confidence depends on the inability of the Atwood-St. John modification to toss out one or more of the points, the inclusion of extra points when a nominal ii point design is used, and exploration of d(s,¢). For a much simpler 5parameter 2-compartment linear model, several examples of 5, 6 and 7 point D-optimal designs have been found. 6.

ACKNOWLEDGMENT

This research was supported in part by NSF Grant ECS-8015965, USPHS Grant CA16042 and NIH NRSA F32-GM05519.

7.

REFERENCES

[1]

Cobelli, C. and D i Stefano, J.J., Parameter and structural identifiability concepts and ambiguities: a critical review and analysis, Am. J. of Physiology 239 (1980) R7-R24.

[2]

Di $tefano, J.J., Design and optimization of tracer experiments in ~hysiology and medicine, Fed. Proceedings 39 (1980) 84-90.

[3]

Godfrey, K.R., Jones, R.P., and Brown, R.F., Identifiable pharmacokinetfc models: the role of extra inputs and measurements, J. of Pharmacokin. Biopharm. 8 (1980) 633-648.

[4] Nguyen, V.V. and Wood, E.F., Review and unification of linear identifiability concepts, SIAM Review 24 (1982) 34-51.

[5]

Brown, R.F. determinacy application Biosciences

and Godfrey, K.R., Problems of in compartmental modeling with to bilirubin kinetics, Math. 40 (1978) 205-224.

[6] Metzler, C.M., Elfring, G.L., and McEwen, A. J., A package of computer programs for pharmacokinetic modeling, Biometrics 30 (1974) 562.

[7]

D'Argenio, D.Z., Optimal sampling times for pharmacokinetic experiments, J. of Pharmacokin. Biopharm. 9 (1981) 739-755.

[8]

Di Stefano, J.J., Experience with sequential optimal sampling schedule designs for pharmacokinetic and physiologic experiments, in Ames, W.F. and Vichnevetsky, R. (eds.) Proceedings of the lOth IMACS World Congress on Systems Simulation and Scientific Computation (1982).

530 [9]

E.M. Landaw / Optimal multicompartmentalsampling designs for parameter estimation Mori, F. and Di Stefano, J.J., Optimal nonuniform sampling interval and test-input design for identification of physiological systems from very limited data, IEEE Trans. Autom. Control 24 (1979) 893-900.

[i0] Landaw, E.M., Optimal Experimental Design for Biologic Compartmental Systems with Applications to Pharmacokinetics, Ph.D. Thesis, Dept. of Biomathematics, UCLA (University Microfilms International 1980). [ii] Fedorov, V.V., Theory of Optimal Experiment Design (Academic Press, New York, 1972). [12] Silvey, S.D., Optimal Design (Chapman Hall, New York, 1980). [13] Goodwin, G.C. and Payne, R.L., Dynamic System Identification: Experiment Design and Data Analysis (Academic Press, New York, 1977). [14] St. John, R.C. and Draper, N.R., D-optimality for regression designs: a review, Technometrics 17 (1975) 15-23. [15] Atwood, C.L., Optimal and efficient designs of experiments, Annals of Math. Stat. 40 (1969) 1570-1602. [16] Cobelli, C. et. at., Model of the kinetics of ketone bodies in humans, Am. J. of Physiology 243 (1982) R7-RI7. [17] Box, M.J., Some experiences with a nonlinear experimental design criterion, Technometrics 12 (1970) 569-589. [18] Wynn, H.P., Results in the theory and construction of D-optimum experimental designs. J. Royal Statistical Soc., Series B, 34 (1972) 133-147, discussion 170-186. [19] Cook, R.D. and Nachtsheim, C.J., Model Robust Linear-Optimal Designs, Technometrics 24 (1982) 49-54.