Ind. Eng. Chem. Res. 2004, 43, 1133-1143
1133
Combined Deterministic-Stochastic Approach for Pharmacokinetic Modeling Dong-Yup Lee,†,‡,§ Su Whan Sung,‡ Sang Yup Lee,†,‡,§ and Sunwon Park*,‡ Metabolic and Biomolecular Engineering National Research Laboratory, Department of Chemical and Biomolecular Engineering and BioProcess Engineering Research Center, and Department of BioSystems and Bioinformatics Research Center, Korea Advanced Institute of Science and Technology, 373-1 Guseong-dong, Yuseong-gu, Daejeon 305-701, Korea
We develop a systematic methodology for pharmacokinetic (PK) modeling to predict drug exposure profiles accurately. First, the number of compartments is determined by inspecting the singular values of an augmented matrix composed of the bolus (or impulse) drug response in the central compartment. The prediction error identification method is then applied to estimate model parameters based on the proposed deterministic-stochastic model structure in which stochastic dynamics and various dosing strategies as well as physiological delays are incorporated simultaneously. The deterministic part of the model is represented by a physical continuoustime compartmental model, whereas the stochastic part is a discrete-time empirical finite impulse response form. Three examples demonstrate that the proposed modeling strategy not only provides a good criterion to determine the appropriate compartment number but also describes well-combined deterministic-stochastic dynamics of the PK system with optimally estimated model parameters. 1. Introduction Recent advances in recombinant biotechnology of genomics, high throughput screening, and combinatorial chemistry have accelerated the drug discovery process. As a result, the clinical development processes have been the bottleneck and a major challenge in the research and development of new drug entities. Many researchers have been exerting much of their efforts in developing pharmacokinetic (PK) models in the past decade because understanding pharmacological profiles and physiological effects have been recognized to be essential in such clinical drug-development processes.1-3 Several available approaches for PK modeling include empirical approaches,4,5 physiological models,6-9 and compartmental models.5,10-12 In this study, we focus on the compartmental models because the models can describe the dynamic dose-exposure-effect relationship more systematically compared to simple empirical models. Furthermore, they have advantages in “parsimony principle” compared to physiological models that require many unavailable system parameters.1 Previous compartmental models can describe various physiological phenomena such as nonlinear dynamics caused by binding processes, cosubstrate depletion, product inhibition, etc.5,10,12 In the modeling of any dynamic systems, it is significant to determine the model structure as concise as possible according to the parsimony principle: If we use more compartments than needed, the subsequent model parameter estimation not * To whom correspondence should be addressed. Tel.: +8242-869-3920. Fax: +82-42-869-3910. E-mail: sunwon@ kaist.ac.kr. † Metabolic and Biomolecular Engineering National Research Laboratory. ‡ Department of Chemical and Biomolecular Engineering and BioProcess Engineering Research Center. § Department of BioSystems and Bioinformatics Research Center.
only becomes difficult but also fails because of singular value problems, and even the estimated parameters of the overparametrized model structure are physically meaningless. Moreover, the modeling errors, especially in the extrapolation region, may be unacceptable; thus, a systematic method is highly required to determine the number of compartments in the compartmental modeling. One of the easiest and most efficient previous methods is the graphical approach that determines the number of the compartments by checking the number of exponential phases in the plot of concentration versus time data in the semilog scale.11 However, it cannot discriminate similar time constants (equivalently, similar exponential phases) of the dynamic system. The previous methods also do not consider stochastic dynamics attributable to system noises, structural process/model mismatches, and physiological disturbances. Several researchers11,13-16 tried to incorporate the stochastic dynamics in estimating the model parameters. They assumed the probability density function of the residual between the process output and the model output and then estimated the model parameters by maximizing the overall likelihood function. However, all of the papers except one13 cannot be applied if the time sequence of the residual is correlated. Although D’Argenio and Park13 can theoretically make the residual uncorrelated using the Kalman filter, it requires a heavy computational load and prior information on the variance of the system noises and the measurement noises. In this work, we develop a systematic methodology for PK modeling on the basis of a proposed deterministic-stochastic model structure. Sections 2 and 3 contain the integrated model structure for more systematic compartmental modeling by incorporation of the underlying problems mentioned above. Effective methodologies for determining the number of compartments are discussed in section 4. A systematic parameter estimation method for the integrated model structure
10.1021/ie0305364 CCC: $27.50 © 2004 American Chemical Society Published on Web 01/23/2004
1134 Ind. Eng. Chem. Res., Vol. 43, No. 4, 2004
determined is introduced in section 5. The efficacy of the proposed modeling strategy is demonstrated by applying it to three examples in section 6. Subsequently, several PK modeling issues for the methodology developed are discussed in section 7. Finally, the main conclusions are summarized in section 8. 2. Deterministic Multicompartmental Modeling for PK Systems Consider the deterministic multicompartmental model for the PK system as shown in Figure 1. It is composed of the central, extravascular, and peripheral compartments. The drug concentrations of the extravascular and central compartments are controlled by single/multiple bolus and infusion drug dosing. Note that if we can successfully integrate these types of compartments and the multibolus/infusion drug dosing within one framework, more complicated multicompartmental models can be formulated in a similar manner. This will be discussed further in the Discussion section. The dynamics for each compartment can be mathematically described by following differential equations.
Central compartment Vc
dC ˆ c(t) ) -CLcC ˆ c(t) - CLcpC ˆ c(t) + CLpcC ˆ p(t) + dt Di,c(t) + CLecC ˆ e(t-d) (1)
C ˆ c(t) ) C ˆ c(t-) +
Db,c(t) Vc
if t ) tb,c(i), i ) 1, 2, ..., Nb,c (2)
tb,c ) [tb,c(1), tb,c(2), ..., tb,c(Nb,c-1), tb,c(Nb,c)] (3) Extravascular compartment Ve
dC ˆ e(t) ) -(CLec + CLe)C ˆ e(t) + Di,e(t) dt
Db,e(t) ˆ e(t-) + C ˆ e(t) ) C Ve
(4)
if t ) tb,e(i), i ) 1, 2, ..., Nb,e (5)
tb,e ) [tb,e(1), tb,e(2), ..., tb,e(Nb,e-1), tb,e(Nb,e)] (6) Peripheral compartment Vp
dC ˆ p(t) ) -CLpC ˆ p(t) - CLpcC ˆ p(t) + CLcpC ˆ c(t) (7) dt
where Vc, Ve, and Vp are the volumes of the central, extravascular, and peripheral compartments, respectively. C ˆ c, C ˆ e, and C ˆ p represent the drug concentrations of the deterministic model for the central, extravascular, and peripheral compartments, respectively. Db,c(t) and Db,e(t) denote the amounts of the drug entering the central and extravascular compartments, respectively, by the bolus dosing method, while Di,c(t) and Di,e(t) are the amounts entering by other dosing (infusion) methods. tb,c, tb,e and Nb,c, Nb,e are the time schedules for the multibolus input and the number of the bolus injections. d denotes a time delay of drug transportation from the
Figure 1. Integrated model structure for PKs based on a multicompartmental modeling approach.
extravascular compartment to the central compartment, and CL stands for clearance. t is time, and t- is defined by t - , where is an infinitesimally small positive value. The current model formulation can describe various dosing strategies in an integrated framework. First, to represent various routes of drug administration, four different dosing inputs including Db,c(t), Db,e(t), Di,c(t), and Di,e(t) are considered according to different dosing methods and different compartments. Db,c(t) and Db,e(t) represent bolus inputs to the central and extravascular compartments, respectively. The formulation of eqs 2 and 5 represents the effect of the bolus inputs by changing the drug concentration in eqs 1 and 4 instantaneously at the bolus dosing time, tb,c(i) and tb,e(i). On the other hand, Di,c(t) and Di,e(t) can be any arbitrary piecewise continuous-time signals so that they can describe other (e.g., infusion) dosing methodologies except the bolus dosing. It should be noticed that most medical drugs are administered in successive doses rather than a single dose to produce a desired effect and to induce lasting and effectual results. Equations 2, 3, 5, and 6 can express the incorporation of the bolus multidosing by the time schedule of eqs 3 and 6, which means that the bolus injection happens at tb,c(i), i ) 1, 2, ..., Nb,c, and tb,e(i), i ) 1, 2, ..., Nb,e. 3. Extended Formulation To Incorporate Stochastic Uncertainties The previous deterministic model represents only the relationship between dose and exposure. However, PK processes have stochastic dynamics that comes from unmodeled factors affecting transient changes in drugexposure profiles.5 Various physiological and pharmacological conditions or factors can affect this stochastic dynamics of the drug. For example, when the aqueous solubility of the drug is altered by changes of temperature or pH in the case of ionizable drugs, the dissolution rate of a drug is temporarily changed during absorption after oral dosing, resulting in alteration of the plasma concentration. There are also transient increases in drug-exposure profiles caused by a delay in gastric emptying time or changes in the viscosity or pH of the gastrointestinal tract. In addition to these system uncertainties, there exists measurement noise. If this stochastic dynamics is too big to be negligible compared to the deterministic one, we should include them systematically in the model formulation. In this work, to incorporate stochastic uncertainties into the previous deterministic formulation of eqs 1-7, a filtered uncorrelated noise sequence is adopted to
Ind. Eng. Chem. Res., Vol. 43, No. 4, 2004 1135
describe the dynamics of the stochastic uncertainties representing the effects of all factors other than the drug dose on the drug concentration (for details of the mathematical description, see refs 17-19). Even though the filtered noise sequence cannot describe all characteristics of stochastic noises, it has been recognized as acceptable for most practical purposes. In addition, any physiological stochastic systems with piecewise constant uncorrelated input noises can be represented by the filtered noise sequence form if the systems are stable and linear.
Cc(t) ) C ˆ c(t) + X(t)
(8)
X(t) ) e(t) + c1e(t-∆t) + c2e(t-2∆t) + ... + cnce(t-nc∆t) (9) where C ˆ c and Cc are the drug concentration of the deterministic model and the measured drug concentration, respectively, in the central compartment. X(t) denotes the stochastic dynamics driven by an uncorrelated noise e(t), which is assumed as a piecewise constant uncorrelated noise with a time interval ∆t as follows:
e(t) ) e(ti), ti e t < ti + ∆t
(10)
E[e(ti)] ) 0
(11)
Ee(ti) e(tj) ) 0, i * j
(12)
Ee2(ti) ) σ2
(13)
where eqs 10 and 11 express the piecewise constant and zero-mean noise. Equations 12 and 13 denote the uncorrelation property (i.e., whiteness) between samples and the variance of the Gaussian noise sequence. For a ˆ cdeterministic case, e(t) would be zero; thus, Cc(t) ) C (t). Otherwise, we can describe various dynamic behaviors of stochastic uncertainties by choosing appropriate coefficients ci, i ) 1, 2, ..., nc. If we set c1 ) c2 ) ... ) cnc ) 0, the stochastic model is for the uncorrelated output noises, while if some coefficients are nonzero, it is for the correlated output noises. 4. Determination of the Compartmental Model Structure There are various techniques to determine the model order such as testing ranks in covariance matrices, examining the information matrix, minimizing Akaike’s Information Theoretic Criterion (AIC), and checking the whiteness of residuals. For details, see the reference.18 In this work, we recommend subspace system identification methods20-23 to determine the number of the compartments because they are simple and do not require any iterative procedure. The number of compartments is determined by inspecting the singular value matrix under several conditions. To make use of the continuous-time subspace system identification method,20 the sampling rate of the measurement of the drug concentration should be small enough to guarantee acceptable accuracy in numerical integration, whereas the sampling rate of the drug concentration measurement should be regular to use the discrete-time sub
space system identification methods.21-23 If the magnitude of the sampling rate varies in a wide range, we have to use an iterative approach such as minimizing AIC and checking the whiteness of residuals. The number of compartments can be determined using discrete-time subspace approaches as follows. Equation 14 shows the augmented matrix composed of impulse responses of the central compartment output and its singular value decomposition (SVD),23 where
[
C ˜ (0) C ˜ (1) l C ˜ (P-2) C ˜ (P-1)
C ˜ (1) C ˜ (2) l C ˜ (P-1) C ˜ (P)
‚‚‚ ‚‚‚ l ‚‚‚ ‚‚‚
C ˜ (N-P-1) C ˜ (N-P) l C ˜ (N-3) C ˜ (N-2)
]
C ˜ (N-P) C ˜ (N-P+1) ) l C ˜ (N-2) C ˜ (N-1) USVT (14)
C ˜ (k) is the kth output. N and P denote the total number of the central compartment output data and the assumed maximum range of process order, respectively, under the conditions that the process input is impulse, the sampling time is constant, and there is no uncertainty. The matrices U and V stand for orthogonal matrices. Note that the impulse input is equivalent to the bolus input. Thus, C ˜ (k) simply corresponds to the measured drug concentration in the central compartment when bolus dosing is used and uncertainties are negligible. In the case of infusion dosing and/or nonnegligible uncertainties, C ˜ (k) can be easily reconstructed from the estimated high-order autoregressive-exogeneous-input (ARX) model. The procedure is simple and does not require any iterative procedure. For details, see ref 23. It was proven that the singular value matrix S is as follows:
[
σ1 0 S) 0 l 0
0 σ2 0 l 0
0 0 σ3 l 0
‚‚‚ ‚‚‚ ‚‚‚ l 0
0 0 0 l σP
0 0 0 l 0
‚‚‚ ‚‚‚ ‚‚‚ ‚‚‚ ‚‚‚
0 0 0 l 0
]
(15)
where σ1 g σ2 g ‚‚‚ g σP and σi ≈ 0, i ) norder + 1, ..., P. Therefore, we can easily infer the process order, i.e., the number of compartments (norder) by inspecting the first singular value that is close to zero, i.e., σnorder+1. Although we recommend the SVD method in this work, the previous graphical method11 is a good alternative because the number of compartments can be easily determined by checking the number of exponential phases in the plot of concentration versus time data in the semilog scale. The graphical method, however, may have some problems in discriminating similar time constants of the process. For example, consider the twocompartment model system with two sets of system parameters (Figure 2).
Set 1: (CL1 + CL12)/V1 ) 1.0/h, CL21/V1 ) 0.8/h, CL12/V2 ) 0.8/h, (CL2 + CL21)/V2 ) 1.0/h Set 2: (CL1 + CL12)/V1 )1.0 /h, CL21/V1 ) 0.2/h, CL12/V2 ) 0.2/h, (CL2 + CL21)/V2 ) 1.0/h
1136 Ind. Eng. Chem. Res., Vol. 43, No. 4, 2004
Figure 2. Two-compartmental system for demonstrating the graphical method and the SVD method.
are defined as the reciprocals of the eigenvalues of the matrix in eq 16. As shown in Figure 3a, the semilog plot shows an almost one-compartmental system for the similar time constants (set 2) and a two-compartmental system for the totally different time constants (set 1). This indicates that the graphical approach may fail to determine the number of compartments for the systems with similar time constants. On the other hand, the recommended SVD method predicts fairly well the number of compartments for both cases as shown in Figure 3b. Nevertheless, both approaches can be complementarily employed according to the researcher’s preference. 5. Model Parameter Estimation Once the model structure is determined as described above, we should estimate the model parameters to fit the experimental data sets. In this work, the prediction error identification method is exploited to minimize the following norm of the prediction error.17
[
min V(θˆ ) ) θˆ
N
1
∑L(Cc(ti)-Ch c(ti),θˆ ,ti)
Ni)1
]
(17)
subject to Central compartment dC ˆ c(t) ) -θˆ 1C ˆ c(t) + θˆ 2C ˆ p(t) + θˆ 3Di,c(t) + θˆ 4C ˆ e(t-θˆ 5) dt (18) C ˆ c(t) ) C ˆ c(t-) + θˆ 3Db,c(t)
if t ) tb,c(i), i ) 1, 2, ..., Nb,c (19)
Extravascular compartment dC ˆ e(t) ˆ e(t) + θˆ 7Di,e(t) ) -θˆ 6C dt C ˆ e(t) ) C ˆ e(t-) + θˆ 7Db,e(t)
(20)
if t ) tb,e(i), i ) 1, 2, ..., Nb,e (21)
Peripheral compartment dC ˆ p(t) ) -θˆ 8C ˆ p(t) + θˆ 9C ˆ c(t) dt Figure 3. Graphical method and the SVD method for determining the number of compartments: (a) semilogarithmic plot of concentration versus time for the graphical method; (b) singular values of the augmented matrix for the SVD method (solid line, set 1; dotted line, set 2).
The governing equation is as follows:
[
( [ ]
)
CL1 + CL12 CL21 V1 V1 d C1(t) ) CL CL2 + CL21 dt C2(t) 12 V2 V2 -
(
]
(22)
Stochastic uncertainties ˆ c(t) + θˆ 10eˆ (t-∆t) + θˆ 11eˆ (t-2∆t) + ... + C h c(t) ) C θnθeˆ (t-nc∆t) (23) eˆ (t-∆t) ) Cc(t-∆t) - C ˆ c(t-∆t) - θˆ 10eˆ (t-2∆t) - ... θˆ nθeˆ [t-(nc+1)∆t]K (24)
)[ ] C1(t) C2(t)
(16)
where C1(0) ) 1.0 g/L and C2(0) ) 0.0 g/L. The case of set 1 gives rise to the time constants of the system, 0.5556 and 5.0000, which are totally different, while two time constants of 0.8333 and 1.250 are similar for the case of set 2. Note that the time constants
System parameters θˆ ) [θˆ 1, θˆ 2, ..., θˆ nθ]T ) d,
[
CLc + CLcp CLpc 1 CLec , , , , Vc Vc Vc Vc
]
T CLec + CLe 1 CLp + CLpc CLcp , , , , c1, c2, ..., cnc Ve Ve Vp Vp (25)
where L(‚) is a scalar-valued positive function. Cc(t) and C h c(t) represent the measured and predicted drug concentrations in the central compartment, respectively.
Ind. Eng. Chem. Res., Vol. 43, No. 4, 2004 1137
Equations 18-22 can be obtained by rearranging eqs 1-9 according to eq 25. Equations 23 and 24 are the one-step-ahead optimal predictor for the stochastic process of eqs 8 and 9. Refer to Appendix A for the derivation of the optimal predictor. Herein, eˆ (t) is the estimated noise for e(t); C ˆ c(t) denotes the drug concentration in the central compartment, calculated by solving the deterministic model, eqs 18-22. The initial values of eˆ (t) for t < 0 are assumed as zeros. As outlined in Appendix B, the adjustable model parameters of eq 25 can be estimated by the Levenberg-Marquardt method,24,25 which is most widely adopted for the system identification.11,17,18 It uses the first and second derivatives of the objective function with respect to the model parameters to improve the convergence rate. Remarks. (1) Suppose that given the drug administration data (D(1), D(2), ..., D(N)), the generation of Cc(t) can be described by
Cc(t) ) g(t,Zt-1;θ) + (t,θ)
(26)
where θ and (t,θ) represent the system parameter vector and the prediction error (residual), respectively; Zt-1 denotes {D(1), Cc(1), D(2), Cc(2), ..., D(t-1), Cc(t1)}; and t represents the tth sampling. Here note that the conditional probability density function (CPDF) of (t,θ) for given Zt-1 is f(x,t;θ) if (t,θ) is time-uncorrelated. Otherwise, the CPDF of (t,θ) will be f(x,t,Zt-1;θ). If (t,θ) is time-uncorrelated, the likelihood function (equivalently, the conditional joint probability density function for the observations Cc(1), Cc(2), ..., Cc(N) and given D(1), D(2), ..., D(N)) is N
fCc(θ;ZN) )
f[Cc(t)-g(t,Zt-1;θ),t;θ] ) ∏ t)1 N
f[(t,θ),t;θ] ∏ t)1
(27)
The maximum likelihood method finds the system parameters maximizing the likelihood function of eq 27 as
θˆ ML(ZN) ) max θ
{
1
N
}
log fCc(θ;ZN) ) max θ
{∑ 1
N
Nt)1
log f[(t,θ),t;θ]
}
(28)
If the prediction error is a Gaussian noise (normal distribution) with zero mean values and covariance λ(t), we have
θˆ ML(ZN) ) max θ
{ ( 1
∑
Nt)1
{
) max θ
N
log
1
x2πλ(t)
log 2π
-
2
1
[
exp -
2λ(t)
N
∑log[λ(t)] -
2Nt)1
1
) min θ
{∑ } 1
N
])}
(t,θ)2
N
∑ 2Nt)1
}
(t,θ)2 λ(t)
(t,θ)2
Nt)1 λ(t)
Most previous maximum likelihood methods for PKs solve the minimization problem of eq 29 to obtain the system parameters under the assumption that (t,θ) is time-uncorrelated Gaussian noises. However, if the prediction error (t,θ) is time-correlated, eq 29 cannot be a maximum likelihood estimator because the conditional joint probability density function for the observations Cc(1), Cc(2), ..., Cc(N) and given D(1), D(2), ..., D(N) cannot be described as eq 27, whereas the proposed method renders the residuals time-uncorrelated so that it can be a maximum likelihood estimator for even timecorrelated noises. This is the main difference between previous output error approaches and the proposed stochastic approach. For example, if the positive funch c(ti), θˆ ,ti] ) -log f[Cc(ti)-C h c(ti), θˆ ,ti], tion L(‚) is L[Cc(ti)-C h c(ti), θˆ ,ti) is the probability density where f[Cc(ti)-C h c(ti), the function of the prediction error Cc(ti) - C minimization problem of eq 17 is equivalent to the maximization of the likelihood function of the prediction error. Consequently, the maximum likelihood estimation method is a subset of the aforementioned prediction error estimation method. When the prediction errors are assumed to be Gaussian with zero mean values and a fixed time-invariant covariance, the maximum likelihood method is equivalent to the prediction error h c(ti), method with the quadratic criterion of L[Cc(ti)-C h c(ti)]2.18 θˆ ,ti] ) 0.5[Cc(ti)-C (2) To obtain approximate correlation functions of the ˆ c(t) or X(t)], we can perform the output error [Cc(t) - C following procedure. First, we can identify the deterministic system parameters with the assumption of uncorrelated output noises (that is, with c1 ) c2 ) ... ) cnc ) 0), resulting in the output modeling error [Cc(t) C ˆ c(t)]. Then, the correlations of the resultant output modeling error can be evaluated by RX(τ) ≈ N (1/N)∑t)1 X(t) X(t+τ), τ ) 1, 2, .... Consequently, we can guess how much the output error is correlated and determine the number of the system parameters (nc) for the stochastic modeling. nc corresponds to τj satisfying |RX(τ)| e η for all τ g τj, where η is a user-determined small value (theoretically, zero). (3) Zeroes may be used as the initial values of the stochastic model parameters, which correspond to the uncorrelated assumption: no serious problems such as parameter divergences and convergences to wrong values for single output cases have been encountered in our experiences. If one wants to set the initial values more systematically, the following procedure can be followed. First, obtain the output modeling error [Cc(t) ˆ (t)] by identifying the deterministic system -C ˆ c(t) or X parameters under the assumption of uncorrelated output noises (that is, with c1 ) c2 ) ... ) cnc ) 0). Then, obtain approximate uncorrelated noises [eˆ (t)] after identifying a high-order ARX model. Finally, estimate the initial stochastic model parameters (ci, i ) 1, 2, ..., ˆ (t) ) c1eˆ (tnc) by applying the least-squares method to X ∆t) + c2eˆ (t-2∆t) + ... + cnceˆ (t-nc∆t) with the obtained data sets, eˆ (t) and X ˆ (t). For details, refer to the corresponding reference.18 (4) The expected covariance matrix of the model parameters can be derived as follows:
(29)
{
min V(θˆ ) ) θˆ
0.5 N
}
∑[Cc(ti) - Ch c(ti)]2
N i)1
(30)
where Cc(t) and C h c(t) denote the process output and the model output, respectively.
1138 Ind. Eng. Chem. Res., Vol. 43, No. 4, 2004
E[(θˆ - θ)(θˆ - θ)T] )
Q)
1
N N
∑∑
Ni)1j)1
(
)(
( ) [ ]
)
∂C h c(ti,θ) ∂C h c(tj,θ) ∂θ
-1
2 1 ∂ V(θ) N ∂θ2
∂θ
Q
∂2V(θ)
-T
(31)
∂θ2
T
E[(ti,θ)(tj,θ)] (32)
(ti,θ) ) Cc(ti) - C h c(ti,θ)
(33)
If the prediction error (ti,θ) is a uncorrelated noise of constant variance σ2, eqs 31 and 32 can be simplified as follows:
1
N
∑
Ni)1
(
)(
)
∂C h c(ti,θd) ∂C h c(ti,θd) ∂θd
∂θd
E[(θˆ - θ)(θˆ - θ)T] )
{ ∑[
1 1
N
N Ni)1
)(
T
σ2 ) Q
]}
∂C h c(ti,θ) ∂C h c(ti,θ) ∂θ
∂θ
(34)
T -1 2
σ (35)
For details, refer to ref 17. We can also derive the expected covariance matrix of the model parameters for other types of cost functions through the same procedure of the reference. (5) In the above model parametrization, the uniqueness of the optimal solution for the problem should be carefully considered. As an example, consider θˆ 4, θˆ 6, and θˆ 7 in eqs 18, 20, and 21. If Di,e(t) and Db,e(t) are persistently exciting, there are no uniqueness problem. However, if Di,e(t) ) 0 and Db,e(t) is persistently exciting, we can estimate only θˆ 6 and θˆ 4θˆ 7 terms uniquely through the optimization of eq 17. Therefore, to obtain a unique solution for θˆ 4, θˆ 6, and θˆ 7, we should fix one of the model parameters to a specific value or define a correlation between them. As another example, we realize that only two terms of θˆ 8 and θˆ 2θˆ 9 among θˆ 2, θˆ 8, and θˆ 9 can be uniquely determined when we consider eqs 18 and 22. If we assign a correlation θˆ 8 ) θˆ 9 (no clearance assumption from the peripheral compartment) among θˆ 2, θˆ 8, and θˆ 9, then there is a unique optimal solution for θˆ 2 and θˆ 8 ) θˆ 9. (6) Note that we cannot obtain all physiological parameters from the uniquely estimated system parameters θˆ as shown in eq 25. We should find the missing physiological parameters from the literature. Otherwise, we need drug concentration profiles of other compartments in addition to those of the central compartment in order to obtain all physiological parameters. 6. Case Study Three examples for the PK systems are used to evaluate the parameter estimation performance of the proposed method. Examples 1 and 2 exemplify the determination of the compartment number and the improvement of the model performance by considering stochastic uncertainties for simulated processes. Example 3 shows the identification results of the proposed method for real experimental data sets modified from the literature.26 Example 1. Figure 4 shows observed plasma concentrations with small stochastic uncertainties after two consecutive 100 µg intravenous bolus doses at 0 and 30 h. Data are gathered from the two-compartment PK process simulated with the physiological parameters (Vc,
Figure 4. Observed concentration-time data following two consecutive 100 µg intravenous bolus doses at 0 and 30 h. Data are generated from the simulated two-compartment PK process (Vc, 56.50 L; CLc, 5.97 L/h; CLd, 54.37 L/h; Vp, 56.83 L; c1, 0.5; c2, 0.3; σ2, 0.0001).
Figure 5. example 1.
Singular values of the augmented matrix in
56.50 L; CLc, 5.97 L/h; CLd, 54.37 L/h; Vp, 56.83 L) and the stochastic model uncertainty (c1, 0.5; c2, 0.3; σ2, 0.0001). First, to determine the compartment number, we construct the augmented matrix of eq 14 composed of the measured drug concentrations for the first intravenous bolus dosing; the singular values are obtained as shown in Figure 5 with the SVD of the augmented matrix; and it is clear that the number of the compartments is two. Thus, it is parametrized as shown in Figure 6. Six system parameters are selected from eqs 18-25 to represent the model in Figure 6 as follows:
θˆ ) [θˆ 1, θˆ 2, θˆ 3, θˆ 8, θˆ 10, θˆ 11]T )
[
]
CLc + CLd CLd 1 CLd , , , ,c ,c Vc Vc Vc Vp 1 2
T
)
[p1, p2, p3, p4, p5, p6]T Then, our objective is to estimate the system parameters from the measured drug concentrations and the amount
Ind. Eng. Chem. Res., Vol. 43, No. 4, 2004 1139 Table 1. Model Parameter Estimates during the Optimization for the Combined Deterministic-Stochastic Model in Example 2 iteration
norm of the prediction error V(θˆ )
adjustable parameters [p1, p2, p3, p4, p5, p6]
physiological parameters [Vc, CLc, Vp, CLd]
0 1 2 3 4 5 6 7 8 9 10 11 12
151.31 0.1653 0.0384 0.0558 0.0203 0.0073 0.0062 0.0057 0.0054 0.0053 0.0052 0.0052 0.0052
[1.100, 0.500, 0.700, 0.300, 0.000, 0.000] [1.082, 0.541, 0.020, 0.281, 0.004, 0.002] [0.948, 0.685, 0.022, 0.306, 0.232, 0.220] [0.820, 0.808, 0.019, 0.363, 0.273, 0.251] [0.838, 0.790, 0.020, 0.344, 0.262, 0.245] [0.849, 0.757, 0.019, 0.406, 0.279, 0.255] [0.842, 0.738, 0.019, 0.481, 0.301, 0.267] [0.837, 0.730, 0.018, 0.554, 0.327, 0.278] [0.843, 0.736, 0.018, 0.618, 0.361, 0.287] [0.867, 0.758, 0.018, 0.672, 0.400, 0.293] [0.904, 0.795, 0.018, 0.721, 0.435, 0.297] [0.945, 0.835, 0.018, 0.765, 0.455, 0.298] [0.975, 0.866, 0.018, 0.798, 0.456, 0.294]
[1.429, 0.857, 2.381, 0.714] [49.60, 26.86, 95.43, 26.83] [45.81, 12.06, 102.6, 31.37] [52.14, 0.643, 116.1, 42.12] [50.25, 2.392, 115.3, 39.72] [53.23, 4.880, 99.25, 40.29] [53.78, 5.591, 82.57, 39.71] [54.44, 5.781, 71.77, 39.76] [54.89, 5.896, 65.35, 40.38] [55.06, 5.957, 62.14, 41.76] [55.06, 5.990, 60.74, 43.79] [55.03, 6.009, 60.09, 45.97] [55.04, 6.019, 59.74, 47.64]
Figure 6. Two-compartmental model with physiological parameters and a bolus input of drug in example 1.
Figure 7. Observed (closed circles) and predicted (solid lines) plasma concentration profile(s) for example 1. Estimated parameters are as follows: Vc, 56.20; CLc, 5.95 L/h; CLd, 53.42 L/h; Vp, 58.21 L; c1, 0.0128; c2, 0.409.
of the bolus dosing, which minimizes the quadratic h c(ti),θˆ ,ti) ) criterion of the prediction error, L[Cc(ti)-C h c(ti)]2. The proposed method gives a fairly 0.5[Cc(ti) - C good model performance as shown in Figure 7. Example 2. The observed plasma concentration (closed circle) with large stochastic uncertainties is shown in Figure 8. The simulated process is the same as the one in example 1 except the magnitude of the uncertainties (σ2, 0.005). To describe the dynamic behavior of the observed plasma concentration, we estimated parameters in two models: a deterministic model and a combined deterministic-stochastic model. Figure 8 compares the model performances of the deterministic model and the combined deterministic-
Figure 8. Observed plasma concentration profile (closed circles) and the predicted profile (solid lines): (a) deterministic model (estimated parameters: Vc, 56.20/56.50 L; CLc, 5.95/5.97 L/h; CLd, 53.42/54.37 L/h; Vp, 58.21/56.83); (b) combined deterministicstochastic model (estimated parameters: Vc, 55.04/56.50 L; CLc, 6.02/5.97 L/h; CLd, 47.64/54.37 L/h; Vp, 59.74/56.83; c1, 0.456/0.5; c2, 0.294/0.3).
stochastic model. The deterministic model cannot represent the stochastic behavior of the drug concentration, while the combined deterministic-stochastic model describes it well. Table 1 shows the estimates with respect to the iteration in the Levenberg-Marquardt optimization for the combined deterministicstochastic model, and it is clear that the convergence of the estimates is very fast and stable.
1140 Ind. Eng. Chem. Res., Vol. 43, No. 4, 2004
Figure 10. example 3.
Singular values of the augmented matrix in
Compartmental Model Structure section. In this example, we use the raw data (unequally spaced data) to estimate the system parameters of the same model structure as that in example 1. Figure 9b exhibits the predicted plasma concentration profiles by the deterministic model (dotted line) and by the combined deterministic-stochastic model (solid line). The result shows the small difference between the prediction errors of the two models although the performance of the proposed approach is better than the other, which is due to the fact that example 3 has negligible stochastic dynamics. 7. Discussion
Figure 9. Observed concentration-time data after a 100 µg intravenous bolus dose of compound X in example 3. Closed circles represent the data from ref 26. (a) The solid line is generated by piecewise cubic-spline interpolation. (b) Predicted plasma concentrations profiles for the deterministic model (dotted line; estimated model parameters: Vc, 53.83 L; CLc, 6.59 L/h; CLd, 49.16 L/h; Vp, 57.07) and for the combined deterministic-stochastic model (solid line; estimated model parameters: Vc, 57.29 L; CLc, 6.94 L/h; CLd, 50.93 L/h; Vp, 53.41; c1, 0.1775; c2, 0.0530).
Example 3. Consider the experimental data sets shown in Figure 9. As was done in examples 1 and 2, we will demonstrate the model estimation for the experimental data after we determine the number of compartments using discrete-time process order selection approaches. As mentioned in the Determination of the Compartmental Model Structure section, the sampling time should be regular to use discrete-time approaches for the process order selection. Accordingly, we reconstructed the sampling data sets with a regular sampling time of 0.5 h using the cubic-spline interpolation method. The interpolated data sets are shown in Figure 9a as a solid line. Using the interpolated data sets, we constructed the augmented matrix and obtained the singular values as shown in Figure 10. It is clear that the number of compartments is two. Note that equally spaced data are required when we use the SVD approach to determine the number of compartments. If we do not want to use the interpolated data, other methods are available for determining the number of compartments as mentioned in the Determination of the
In this study, we address the underlying problems in modeling PK systems. They include (1) the determination of the PK model structure and (2) the effective estimation of the model parameters to describe the pharmacological system well. For this purpose, first we have developed an integrated PK model structure. It consists of the central, extravascular, and peripheral compartments to characterize the pharmacological and physiological phenomena. Various dosing forms (e.g., bolus and infusion), different routes of drug administration (e.g., extravascular and intravenous), physiological delay during the drug transportation, and single/multiple dosing can be represented within one integrated framework. If the drug dissolves in body fluids and diffuses through one or more membranes to enter the blood plasma, it corresponds to extravascular routes. The dynamic behavior of the drug through the extravascular routes can be formulated by introducing extravascular compartments as was done in eqs 4-6. For example, the most common extravascular route, oral administration either by solution or by rapid dissolution of solids (i.e., pills) can be represented by bolus input to the extravascular compartment Db,e(t). In addition, the time lag term is introduced between the central compartment and the extravascular compartment in eq 1 to represent physiological delay during drug transportation. In the case of oral administration, bioavailability can also be considered by the fraction of two clearances. On the other hand, intravenous routes can be expressed by input signals for intravenous bolus and infusion corresponding to Db,c(t) and Di,c(t), respectively.
Ind. Eng. Chem. Res., Vol. 43, No. 4, 2004 1141
Note that the deterministic part of the model structure can be any linear differential equation. Therefore, theoretically, the integrated mathematical description of eqs 1-7 is applicable to general linear compartmental systems with various compartment configurations and also various combinations of bolus and infusion dosing in a straightforward manner. However, practically, we should carefully solve the problems pertaining to experimental designs, parsimonious modeling, uniqueness of solutions, parameter identifiability, and so on. From the parsimonious point of view, it is very important to select the number of compartments as well as the number of model parameters reasonably: the simpler model is preferred. The selection can be systematically done by several approaches available.11,18,20-23 Choosing one among the available methods is up to researcher’s preference and resources. After the number of compartments is determined, various compartment configurations11 are possible for the given compartment number, such as the catenary form of successive compartments in series, the mammillary form in which a number of peripheral compartments are attached to one central compartment, and any combinations of the mammillary and catenary forms. However, the best model performances of any configuration for the given number of compartments are the same theoretically if the drug concentration data of the central compartment is available. For a more detailed compartment configuration, we need additional information such as drug concentration data of other compartments. Once the model structure is determined, the next stage is to estimate the model parameters using an optimization method minimizing the difference between the observed and predicted concentrations. Previously, several techniques for the parameter estimation have been proposed with focus on the criteria for the best fit. They include ordinary least squares, weighed least squares, and extended least squares, which are formulated for the objective function to be minimized.11,15 However, they cannot incorporate the stochastic dynamics systematically because their methods are based on deterministic types of models. In this study, we estimate the model parameters by minimizing the prediction error between the experimental data and the optimally predicted output of the combined deterministic-stochastic model. After a good model is obtained to describe the dynamics of the drug concentration using the proposed approaches, we can design an optimal dosage regimen through simulation and optimization to yield a treatment effect profile that is maintained as closely as possible to a target over the treatment period. The stochastic modeling would contribute to a more accurate prediction of the future drug concentration dynamics, thus resulting in tighter drug concentration control. We leave it as a future work to develop an optimal dosing strategy based on the proposed integrated model structure. 8. Concluding Remarks We have developed a systematic methodology for PK modeling. First, the integrated model structure is proposed to incorporate various dosing strategies and stochastic dynamics. Then, an efficient SVD approach is introduced for determining the number of compartments. Finally, the prediction error identification method
is applied to estimate model parameters of the integrated model structure. Three examples confirm the efficacy and efficiency of the proposed approach. Although the main focus of this study lies on the PK modeling, similar procedures for stochastic formulation and parameter estimation can be extended to pharmacodynamics. The results in the study can be used to construct the optimal dosage regimen for use in clinical trials and drug therapy. Acknowledgment This work was partially supported by the Advanced Backbone IT Technology Development Project of the Korean Ministry of Information and Communication and Ministry of Science and Technology, the BK21 project, and the Center for Ultramicrochemical Process Systems sponsored by KOSEF. Appendix A This appendix contains the derivation of the one-stepahead optimal predictor to predict the drug concentration Cc(t) from the past measured drug concentration and the past dosing inputs in an optimal manner. Consider the stochastic process of eqs 8 and 9. From eq 9, we can derive the time series equation to calculate the past noise sequence
eˆ (t-∆t) ) Cc(t-∆t) - C ˆ c(t-∆t) - c1eˆ (t-2∆t) - ... cnceˆ [t-(nc+1)∆t] (A1) where the initial values of eˆ (t) for t < 0 are assumed as zeros. If the Euclidean norms of all zeros of the polynomial H(q) ) 1 + c1q-1 + c2q-2 + ... + cncq-nc are less than 1.0, the effects of the wrong initial values disappear and the estimated noise sequence stably converges to the true value as time goes on. Even some Euclidean norms of the zeros of the actual noise polynomial are bigger than 1.0; those of the estimated noise polynomial are less than 1.0. Therefore, we do not have to worry about the divergence of eq A1 and the initial value problem. Now, we can derive the following equation because we calculated the past noise sequence using eq A1.
Cc(t) ) e(t) + c1eˆ (t-∆t) + c2eˆ (t-2∆t) + ... + cnceˆ (t-nc∆t) (A2) In eq A2, we cannot predict the exact value of e(t) because it is a future noise. Instead, we may know the most probable value of e(t) on the basis of the probabilistic properties. In this work, we choose the expectation value among the properties to determine the most probable value. Thereby, we can obtain the optimal onestep-ahead predictor from eq 11 as follows:
C h c(t) ) c1eˆ (t-∆t) + c2eˆ (t-2∆t) + ... + cnceˆ (t-nc∆t) (A3) Appendix B This appendix outlines the Levenberg-Marquardt method24,25 adopted to solve the optimization problem as done by Ljung18 and Sung and Lee,17 which repeats
1142 Ind. Eng. Chem. Res., Vol. 43, No. 4, 2004
eq B1 until the model parameters converge within the tolerance,
θˆ (j) ) θˆ (j-1) -
{
∂2V[θˆ (j-1)] ∂θˆ
2
][ -1
+ RI
}
∂V[θˆ (j-1)} ∂θˆ
(B1) where j denotes the iteration number and R is a small positive value that can be updated every iteration to compromise between the robustness and the convergence rate. The first and second derivatives of the cost function can be derived from eq 17 as follows:
∂V(θˆ ) ∂θˆ ∂2V(θˆ )
)-
∂θˆ 2
1
N
∑
)-
1
h c(ti) dL ∂C
N
∑ h Ni)1dC
(B2)
∂θˆ
c
+ ∂θˆ 2 h c(ti) 1 N d2L ∂C
∑
Ni)1dC h 2 c
[
[ ][ ] ] ∂C h c(ti)
∂θˆ
∂θˆ
∂C h c(t) h c(t) ∂C h c(t) ∂C ∂C h c(t) ∂C h c(t) ) ... ∂θˆ ∂θˆ 1 ∂θˆ 2 ∂θˆ nθ - 1 ∂θˆ nθ
∂ V(θˆ ) ∂θˆ 2
≈
1
N
∑
[ ][ ]
h c(ti) d2L ∂C
Ni)1dC h 2 c
∂θˆ
T
(B3)
T
(B4)
∂C h c(ti) ∂θˆ
T
(B5)
Herein, we need to calculate ∂C h c(t)/∂θˆ in order to obtain the derivatives of the cost function of eqs B2 and B5. To do so, we use the following equations (B6) and (B7), which are derived from eqs 23 and 24. where the
{
∂C h c(t)
) ∂θˆ i ∂C ˆ c(t) ∂eˆ (t-nc∆t) ∂eˆ (t-∆t) ∂eˆ (t-2∆t) + θˆ 10 + θˆ 11 + ... + θˆ nθ ∂θˆ i ∂θˆ i ∂θˆ i ∂θˆ i i ) 1, 2, ..., 9 ∂eˆ (t-nc∆t) ∂eˆ (t-∆t) ∂eˆ (t-2∆t) θˆ 10 + θˆ 11 + ... + θˆ nθ + eˆ [t-(i-9)∆t] ∂θˆ i ∂θˆ i ∂θˆ i i ) 10, 11, ..., nθ
{
ˆ p(t) ∂C ˆ p(t) ∂C ˆ c(t) d ∂C ) -θ8 + θ9 dt ∂θˆ 1 ∂θˆ 1 ∂θˆ 1
(B9)
In a similar way, we can derive all necessary differential equations to calculate the derivatives of the drug concentration.
( ) ( ) ( ) ( )
ˆ p(t) ∂C ˆ p(t) ∂C ˆ c(t) d ∂C ) -θ8 + θ9 dt ∂θˆ 2 ∂θˆ 2 ∂θˆ 2
If the model parameters are close to the real ones, the magnitude of the first term in eq B3 is negligible compared to the second term in eq B3. Then, we can approximate the second derivative of eq B3 using the following equation. 2
( ) ( )
ˆ c(t) ∂C ˆ c(t) ∂C ˆ p(t) d ∂C ) -θˆ 1 + θˆ 2 -C ˆ c(t) (B8) dt ∂θˆ 1 ∂θˆ 1 ∂θˆ 1
ˆ c(t) ∂C ˆ c(t) ∂C ˆ p(t) d ∂C ) -θˆ 1 + θˆ 2 -C ˆ p(t) (B10) dt ∂θˆ 2 ∂θˆ 2 ∂θˆ 2
2 h c(ti) dL ∂ C
Ni)1dC hc
to calculate the derivatives of the deterministic drug concentration with respect to each model parameter θi [∂C ˆ c(t)/∂θi], which can be derived from the differential equations (18)-(22). For example, the following equations are used to calculate ∂C ˆ c(t)/∂θ1.
(B6)
∂eˆ (t-∆t) ) ∂θˆ i ∂eˆ [t-(nc+1)∆t] ∂C ˆ c(t-∆t) ∂eˆ (t-2∆t) - θˆ 10 - ... - θˆ nθ ∂θˆ i ∂θˆ i ∂θˆ i i ) 1, 2, ..., 9 ∂eˆ [t-(nc+1)∆t] ∂eˆ (t-2∆t) -θˆ 10 - ... - θˆ nθ - eˆ [t-(i-8)∆t] ∂θˆ i ∂θˆ i i ) 10, 11, ..., nθ (B7)
h c(t)/∂θˆ , where t initial values of ∂eˆ (t)/∂θˆ , ∂C ˆ c(t)/∂θˆ , and ∂C < 0, are set to zeros. Now, to obtain ∂C h c(t)/∂θˆ , we need
(B11)
ˆ c(t) ∂C ˆ c(t) ∂C ˆ p(t) d ∂C ) -θˆ 1 + θˆ 2 + Di,c(t) (B12) dt ∂θˆ 3 ∂θˆ 3 ∂θˆ 3
∂C ˆ c(t-) ˆ c(t) d ∂C ) + Db,c(t) dt ∂θˆ 3 ∂θˆ 3
( )
if t ) tb,c(i), i ) 1, 2, ..., Nb,c (B13)
ˆ p(t) ∂C ˆ p(t) ∂C ˆ c(t) d ∂C ) -θ8 + θ9 dt ∂θˆ 3 ∂θˆ 3 ∂θˆ 3
( )
(B14)
ˆ c(t) ∂C ˆ c(t) ∂C ˆ p(t) d ∂C ) -θˆ 1 + θˆ 2 +C ˆ e(t-θˆ 5) (B15) dt ∂θˆ 4 ∂θˆ 4 ∂θˆ 4
( )
ˆ p(t) ∂C ˆ p(t) ∂C ˆ c(t) d ∂C ) -θ8 + θ9 dt ∂θˆ 4 ∂θˆ 4 ∂θˆ 4
( )
(B16)
ˆ c(t) ∂C ˆ c(t) ∂C ˆ p(t) dC ˆ e(t-θˆ 5) d ∂C ) -θˆ 1 + θˆ 2 - θˆ 4 dt ∂θˆ 5 ∂θˆ 5 ∂θˆ 5 dt (B17) ∂C ˆ (t) ∂C ˆ (t) ∂C ˆ (t) d p p c + θ9 (B18) ) -θ8 dt ∂θˆ 5 ∂θˆ 5 ∂θˆ 5
( )
( )
ˆ c(t) ∂C ˆ c(t) ∂C ˆ p(t) ∂C ˆ e(t-θˆ 5) d ∂C ) -θˆ 1 + θˆ 2 + θˆ 4 dt ∂θˆ 6 ∂θˆ 6 ∂θˆ 6 ∂θˆ 6 (B19) ˆ e(t) ∂C ˆ e(t) d ∂C ) -θ6 -C ˆ e(t) (B20) dt ∂θˆ 6 ∂θˆ 6
( ) ( )
ˆ p(t) ∂C ˆ p(t) ∂C ˆ c(t) d ∂C ) -θ8 + θ9 (B21) dt ∂θˆ 6 ∂θˆ 6 ∂θˆ 6 ˆ c(t) ∂C ˆ c(t) ∂C ˆ p(t) ∂C ˆ e(t-θˆ 5) d ∂C ) -θˆ 1 + θˆ 2 + θˆ 4 dt ∂θˆ 7 ∂θˆ 7 ∂θˆ 7 ∂θˆ 7 (B22) ∂C ˆ (t) ∂C ˆ (t) e e d ) -θˆ 6 + Di,e(t) (B23) dt ∂θˆ 7 ∂θˆ 7
( )
( )
( )
∂C ˆ e(t-) ˆ e(t) d ∂C ) + Db,e(t) dt ∂θˆ 7 ∂θˆ 7
if t ) tb,e(i), i ) 1, 2, ..., Nb,e (B24)
Ind. Eng. Chem. Res., Vol. 43, No. 4, 2004 1143
( ) ( ) ( ) ( ) ( )
ˆ p(t) ∂C ˆ p(t) ∂C ˆ c(t) d ∂C ) -θ8 + θ9 dt ∂θˆ 7 ∂θˆ 7 ∂θˆ 7 ˆ c(t) ∂C ˆ c(t) ∂C ˆ p(t) d ∂C ) -θˆ 1 + θˆ 2 dt ∂θˆ 8 ∂θˆ 8 ∂θˆ 8
ˆ p(t) ∂C ˆ p(t) ∂C ˆ c(t) d ∂C ) -θ8 + θ9 -C ˆ p(t) dt ∂θˆ 8 ∂θˆ 8 ∂θˆ 8 ˆ c(t) ∂C ˆ c(t) ∂C ˆ p(t) d ∂C ) -θˆ 1 + θˆ 2 dt ∂θˆ 9 ∂θˆ 9 ∂θˆ 9
ˆ p(t) ∂C ˆ p(t) ∂C ˆ c(t) d ∂C ) -θ8 + θ9 +C ˆ c(t) dt ∂θˆ 9 ∂θˆ 9 ∂θˆ 9
(B25) (B26)
(B27)
(B28)
(B29)
In summary, to obtain eqs B2 and B5, we calculate C h c(ti) and ∂C h c(t)/∂θˆ |t)ti from eqs 23 and B6 for given θˆ (j-1) by selecting them at every ti while continuously solving the ordinary differential equations (18)-(24) and (B8)-(B29) simultaneously. Next, it is straightforward to calculate the updated parameters θˆ (j) from eq B1. Repeat this procedure until the parameters converge. Literature Cited (1) Balant, L. P.; Gex-Babry, M. Modelling During Drug Development. Eur. J. Pharm. Biopharm. 2000, 50, 13. (2) Panchagnula, R.; Thomas, N. S. Biopharmaceutics and Pharmacokinetics in Drug Research. Int. J. Pharm. 2000, 201, 131. (3) Rooney, K. F.; Snoeck, E.; Watson, P. H. Modelling and Simulation in Clinical Drug Development. Drug Discov. Today 2001, 6, 802. (4) Wise, M. E. Negative Power Functions of Time in Pharmacokinetics and Their Applications. J. Pharmacokinet. Biopharm. 1985, 13, 309. (5) Kwon, Y. Handbook of Essential Pharmacokinetics, Pharmacodynamics and Drug Metabolism for Industrial Scientists; Kluwer Academic/Plenum Publishers: New York, 2001. (6) Jarvis, D. A. Physiological Pharmacokinetic ModelssA Review of Their Principles and Development. Anaesth. Pharmacol. Rev. 1994, 2, 214. (7) Weiss, M. Pharmacokinetics in Organs and the Intact Body: Model Validation and Reduction. Eur. J. Pharm. Sci. 1998, 50, 13. (8) Lau, C.; Andersen, M. E.; Crawford-Brown, D. J.; Kavlock, R. J.; Kimmel, C. A.; Knudsen, T. B.; Muneoka, K.; Rogers, J. M.; Setzer, R. W.; Smith, G.; Tyl, R. Evaluation of Biologically Based Dose-Response Modeling for Developmental Toxicity: A Workshop Report. Regul. Toxicol. Pharmacol. 2000, 31, 190.
(9) Mahfouf, M.; Linkens, D. A.; Xue, D. A New Generic Approach to Model Reduction for Complex Physiologically Based Drug Models. Control Eng. Pract. 2002, 10, 67. (10) Sun, Y.-N.; Jusko, W. J. Transit Compartments versus Gamma Distribution Function to Model Signal Transduction Processes in Pharmacodynamics. J. Pharm. Sci. 1998, 87, 732. (11) Gabrielsson, J.; Weiner, D. Pharmacokinetic/Pharmacodynamic Data Analysis: Concepts and Applications; Swedish Pharmaceutical Press: Stockholm, 2000. (12) Mager, D. E.; Jusko, W. J. Pharmacodynamic Modeling of Time-Dependent Transduction Systems. Clin. Pharmacol. Ther. 2001, 70, 210. (13) D’Argenio, D. Z.; Park, K. Uncertain Pharmacokinetic/ Pharmacodynamic Systems: Design, Estimation and Control. Control Eng. Pract. 1997, 5, 1707. (14) Guzy, S.; Hunt, C. A. Measures of Uncertainty for Pharmacokinetic and Pharmacodynamic Parameter Estimates: A New Computerized Algorithm. Comput. Biomed. Res. 1996, 29, 466. (15) Spilker, M. E.; Vicini, P. An Evaluation of Extended vs Weighted Least Squares for Parameter Estimation in Physiological Modeling. J. Biomed. Inform. 2001, 34, 348. (16) Mesnil, F.; Mentre´, F.; Dubruc, C.; The´not, J.-P.; Mallet, A. Population Pharmacokinetic Analysis of Mizolastine and Validation from Sparse Data on Patients Using the Nonparametric Maximum Likelihood Method. J. Pharmacokinet. Biopharm. 1998, 26, 133. (17) Sung, S. W.; Lee, I. Prediction Error Method for Continuous-time Processes with Time Delay. Ind. Eng. Chem. Res. 2001, 40, 5743. (18) Ljung, L. System Identification: Theory for the User; Prentice Hall: Englewood Cliffs, NJ, 1987. (19) Pintelon, R.; Schoukens, J.; Rolain, Y. Box-Jenkins Continuous-Time Modeling. Automatica 2000, 36, 983. (20) Sung, S. W.; Lee, S. Y.; Kwak, H. J.; Lee, I. ContinuousTime Subspace System Identification Method. Ind. Eng. Chem. Res. 2001, 40, 2886. (21) Larimore, W. E. Canonical Variate Analysis in Identification, Filtering, and Adaptive Control. 29th IEEE Conference on Decision and Control, Honolulu, HI, 1990; p 596. (22) Verhaegen, M. Identification of the Deterministic Part of MIMO State Space Model Given in Innovations Form from InputOutput Data. Automatica 1994, 30, 61. (23) Van Overschee, P.; De Moor, B. N4SID: Subspace Algorithms for the Identification of Combined Deterministic-Stochastic Systems. Automatica 1994, 30, 75. (24) Levenberg, K. A Method for the Solution of Certain Nonlinear Problems in Least Squares. Q. Appl. Math. 1944, 2, 164. (25) Marquardt, D. W. An Algorithm for Least Squares Estimation of Nonlinear Parameters. SIAM J. 1963, 11, 431. (26) Colburn, W. A. A Time-Dependent Volume of Distribution Term Used to Describe Linear Concentration-Time Profiles. J. Pharmacokinet. Biopharm. 1983, 11, 389.
Received for review June 30, 2003 Revised manuscript received December 8, 2003 Accepted December 18, 2003 IE0305364
