Journal of Pharmacokinetics and Pharmacodynamics, Vol. 32, Nos. 3–4, August 2005 (© 2005) DOI: 10.1007/s10928-005-0026-2
Optimal Design for Model Discrimination and Parameter Estimation for Itraconazole Population Pharmacokinetics in Cystic Fibrosis Patients Timothy H. Waterhouse,1,∗ Stefanie Redmann,2 Stephen B. Duffull,2 and John A. Eccleston1 Received March 21, 2005—Final June 2, 2005 Optimal sampling times are found for a study in which one of the primary purposes is to develop a model of the pharmacokinetics of itraconazole in patients with cystic fibrosis for both capsule and solution doses. The optimal design is expected to produce reliable estimates of population parameters for two different structural PK models. Data collected at these sampling times are also expected to provide the researchers with sufficient information to reasonably discriminate between the two competing structural models. KEY WORDS: itraconazole; optimal design; model discrimination; pharmacokinetic modelling; population approach.
INTRODUCTION Background and Significance Allergic bronchopulmonary aspergillosis (ABPA) is a complication of cystic fibrosis (CF), occurring in approximately 10% of patients mostly after the age of 6 years (1). It presents as a hypersensitive reaction to Aspergillus fumigatus (AF) antigens in some CF patients colonised with AF in the lower respiratory tract (2) and requires treatment. Itraconazole has a broad spectrum as an oral antifungal drug. It is active against most fungi including Aspergillus sp. A study has shown that patients 1 School of Physical Sciences, The University of Queensland, Australia. 2 School of Pharmacy, The University of Queensland, Australia. ∗ To whom correspondence should be addressed. Department of Mathematics, The Univer-
sity of Queensland, Brisbane, QLD 4072, Australia; e-mail:
[email protected] 521 1567-567X/05/0800-0521/0 © 2005 Springer Science+Business Media, Inc.
522
Waterhouse, Redmann, Duffull, and Eccleston
achieved an improvement in the inflammatory response when itraconazole was given in combination with systemic glucocorticosteroids (3). Although relatively little is presently known about the pharmacokinetics of itraconazole in patients with CF, there is evidence which suggests that it is different to the PK of the drug in non-CF patients (4–6). It is also known that there is often large variability in the pharmacokinetic parameters (e.g. absorption rate constant, volume of distribution, and clearance) of drugs administered to CF patients (7,8). It has been observed that high doses of itraconazole are required in CF patients for the treatment of ABPA (2). Itraconazole (Fig. 1) is a highly lipophilic drug (log octanol-water partition coefficient = 5.66, pH 8.1), and is considered ‘insoluble’ in water (less than 5 µg/ml) and has a molecular weight of 705.64 (9). It is a weak base and is ionised only at low pH (pka = 3.7) (10). Itraconazole is widely distributed in the body, with a reported volume of distribution of 10.7 l/kg following intravenous administration (9). Up to 99.8% of the drug is bound to human plasma proteins. Itraconazole is extensively metabolised by the hepatic cytochrome P450 isoenzymes and has at least 30 metabolites, each representing less than 1–5% of the dose. The main metabolite is considered to be hydroxyitraconazole (also shown in N Cl
N N
Cl
H2C
Itraconazole
O
O H
CH3
O C H2
O N
N
CH3
N
N
N N Cl N
Hydroxyitraconazole Cl
H2C O
O
CH3
O H
CH3
O N
N
N
N N
Fig. 1. Chemical structures of itraconazole and hydroxyitraconazole.
OH
Optimal Design for Model Discrimination and Parameter Estimation
523
Fig. 1), which also exhibits antifungal activity against a variety of species (9,11,12). Several studies have shown that the absorption of the capsule formulation is enhanced by a low gastric pH, and taking the capsules after food or with an acidic beverage (e.g. cola) can improve the absorption (9,13,14). The absorption of the oral solution is reported not to be influenced greatly by the degree of gastric acidity (15), and the absorption is more rapid in the fasting state (16–19). It has been shown in healthy volunteers that the oral availability of itraconazole oral solution taken under these conditions can be up to 60% higher than that of capsules (15,17,20). The pharmacokinetics of itraconazole and its main metabolite has been shown to exhibit the behaviour of two linked compartmental models. Figure 2 shows the general structure of the models, including the rates of mass transfer between each compartment. The amounts of itraconazole in the gut, central and peripheral compartments are given by A1 , A2 and A3 , respectively; the amount of its metabolite (hydroxyitraconazole) is denoted by A4 . The elimination from the central compartment to a compound other than hydroxyitraconazole (shown in Fig. 2) has been described by
Fig. 2. Compartmental model of the pharmacokinetics of itraconazole and hydroxyitraco nazole.
524
Waterhouse, Redmann, Duffull, and Eccleston
some as linear (elimination rate F20 C L 2 /V2 )(10), and by others as nonlinear (elimination rate F20 Vmax /(K m + A2 /V2 )) (9,14,19,21). Structural Models The differential equations describing the two models (linear and nonlinear) from Fig. 2 are shown in Eqs. (1–8). Model 1. The linear model: d A1 = −F12 ka A1 − F14 ka A1 , dt
(1)
C L 24 C L2 d A2 Q Q = F12 ka A1 + A3 − A2 − F24 A2 − F20 A2 , dt V3 V2 V2 V2
(2)
Q d A3 Q = A2 − A3 , dt V2 V3
(3)
C L 24 d A4 C L4 = F14 ka A1 + F24 A2 − A4 . dt V2 V4
(4)
Model 2. The nonlinear model: d A1 = −F12 ka A1 − F14 ka A1 , dt
(5)
C L 24 Vmax d A2 Q Q = F12 ka A1 + A3 − A2 − F24 A2 − F20 A2 , dt V3 V2 V2 K m + A2 /V2 (6)
Q d A3 Q = A2 − A3 , dt V2 V3
(7)
C L 24 d A4 C L4 A2 − A4 . = F14 ka A1 + F24 dt V2 V4
(8)
Optimal Design for Model Discrimination and Parameter Estimation
525
Note that since the differential equations describe the rates of change of the amount of drug in each compartment (rather than the concentration), the last term in Eq. (6) describing the Michaelis–Menten elimination is slightly different to the form which may be familiar to the reader. The divisor still involves the concentration (i.e. A2 /V2 ), but we multiply by A2 to give the rate of change in the amount. The units of Vmax and K m are still of the form ‘mass/(volume × time)’ and ‘mass/volume’, respectively. In both sets of differential equations, the initial conditions are A1 (0) = dose,
A2 (0) = A3 (0) = A4 (0) = 0.
(9)
Aim This paper describes the development of an optimal design in real time for a population pharmacokinetic study of itraconazole when given as capsules and solution in adult patients with CF.
THEORY We are interested in modelling the behaviour over time of the amount of the parent drug in the central compartment and its metabolite ( A2 and A4 ), and hence their concentrations (C2 = A2 /V2 and C4 = A4 /V4 ). It is assumed that the fractions F12 , F14 , F24 and F20 are known fixed constants and that their values are available from previous experiments. The parameters to be estimated which are common to both models are therefore the volumes V2 , V3 and V4 ; and the clearances C L 4 , C L 24 and Q. The linear model has one additional parameter, C L 2 , while the nonlinear model also has parameters Vmax and K m . Considering now that we are interested in the pharmacokinetics of itraconazole in both capsule and solution form, there are actually four models under consideration: the linear model for capsules, the linear model for solutions, the nonlinear model for capsules and the nonlinear model for solutions. The models for capsules and solutions have the same structure, they only differ in terms of their parameter values. Hence there are four sets of parameter values:
1,c 1,c T 1,c β 1,c = (ka1,c , V21,c , V31,c , V41,c , C L 1,c 4 , C L 24 , Q , C L 2 ) ,
(10)
2,c 2,c 2,c 2,c T β 2,c = (ka2,c , V22,c , V32,c , V42,c , C L 2,c 4 , C L 24 , Q , Vmax , K m ) ,
(11)
526
Waterhouse, Redmann, Duffull, and Eccleston 1,s 1,s T 1,s β 1,s = (ka1,s , V21,s , V31,s , V41,s , C L 1,s 4 , C L 24 , Q , C L 2 ) ,
β 2,s =
2,s (ka2,s , V22,s , V32,s , V42,s , C L 2,s 4 , C L 24 ,
2,s 2,s T Q 2,s , Vmax , Km ) ,
(12) (13)
where for β u, f the superscripts u and f indicate the structural model u(u = 1 for linear, u = 2 for nonlinear) and the form of dose f ( f = c for capsule, f = s for solution), and “T” denotes the vector/matrix transpose. For the time being, suppose that only one structural model and only one form of dose is under consideration. The superscripts u and f are dropped, and the p-vector of parameters is denoted by β. Let θi be the vector of parameter values for the ith individual. Then for a vector of n i sampling times ξi = (ti1 , . . . , tin i )T , suppose that the n i -vector of observations is modelled by yi = η(θi , ξi ) + i .
(14)
Employing a nonlinear mixed effects model, it is assumed that each of the parameters varies randomly between individuals. Let β be the vector of fixed effects and let bi be the vector of random effects for individual i. So the random effects are written as either θi = β + bi (for additive random effects) or θi = β exp(bi ) (for exponential random effects). Thus the structural model can be written in an alternative form, η(β, bi , ξi ). It is assumed that bi ∼ N (0, ), where = diag(ω1 , . . . , ω p ). It is also assumed that, as in Retout and Mentr´e (22), the variances of the normally distributed error terms i (conditional on the random effects) have additive and proportional components, i.e. i |bi ∼ N (0, diag(σ A + σ P η(θi , ξi )2 ) for some σ A and σ P . Denote the entire vector of model parameters by = (β T , ω1 , . . . , ω p , σ A , σ P )T . Optimal Design Model discrimination aside (momentarily), the objective of the experiment is to obtain estimates of the relevant parameters with the greatest possible efficiency. For this purpose, the Fisher information matrix is employed. For a population design = (ξ1 , . . . , ξ N )T this is given by ∂ 2 l( ; y) M F ( , ) = E − , (15) ∂ ∂ T where l( ; y) is the log-likelihood of the vector of observations y = (y1 , . . . , y N )T for the population parameters . The population design is
Optimal Design for Model Discrimination and Parameter Estimation
527
limited to Q groups, such that the Nq subjects of the qth group each have the same set of n q sampling times ξq . So we have M F ( , ) =
Q
Nq M F ( , ξq ).
(16)
q=1
The Cramer–Rao Lower Bound states that the variance-covariance matrix of any unbiased estimator of is bounded below by M F−1 ( , ). It follows that any population design which maximises |M F ( , )| minimises the volume of the joint asymptotic confidence ellipsoid of the parameters. Such a design is called D-optimal. See Atkinson and Donev (23) for an introduction to D-optimal design. For an elementary design ξ , Retout and Mentr´e (22) have shown that by using a first-order Taylor series expansion of η(θ, ξ ) around the expectation of the random effects, an approximation to the Fisher information matrix for an elementary design ξ is given by 1 A(E, S) C(E, S) , (17) M F ( , ξ ) ≈ 2 C T (E, S) B(E, S) where ∂ E T −1 ∂ E S , m, n = 1, . . . , p, (18) ∂βm ∂βn ∂ S −1 ∂ S −1 , m, n = 1, . . . , p + 2, (19) (B(E, S))mn = tr S S ∂λm ∂λn ∂ S −1 ∂ S −1 , m = 1, . . . , p + 2, n = 1, . . . , p. (C(E, S))mn = tr S S ∂λm ∂βn (20) (A(E, S))mn = 2
Here E and S represent the expectation and variance of the observations, respectively: E = E(y) ≈ η(β, 0, ξ ), (21) T ∂η(β, 0, ξ ) ∂η(β, 0, ξ ) S = S(y) ≈ + diag(σ A + σ P η(β, 0, ξ ))2 . T ∂b ∂b T (22)
We assume that S is independent of β, as in Retout et al. (24). We also make the assumption that the fixed and random effects are uncorrelated, i.e. that C(E, S) is a matrix of zeros.
528
Waterhouse, Redmann, Duffull, and Eccleston
Parameter Estimation in Multiresponse Situations Consider for a moment that the observations are in fact modelled by a simpler fixed effects model with multiple response types. Suppose that the outcomes of an experiment on a single subject with r simultaneous responses can be modelled as (a)
yi
(a)
= η(a) (θ, ti ) + i
(a = 1, . . . , r ; i = 1, . . . , n),
(23)
where the θ = (θ1 , . . . , θ p )T is now a p-vector of fixed effects parame(a) 2 2 ters; E(i(a) ) = 0; E(i(a) (b) j ) = 0 for i = j; E((i ) ) = σa = σaa ; (a) (b)
and E(i i ) = ρab σa σb = σab for a = b. So the covariance matrix of responses for the i th sampling time is given by σ11 σ12 · · · σ1r σ21 σ22 · · · σ21r (24) = {σab }a,b=1,...,r = . .. . . .. .. . . . σr 1 σr 21 · · · σrr
where σab = σba . Let −1 = {σ ab }a,b=1,...,r . Again, denote the elementary design by ξ = (t1 , . . . , tn )T . Draper and Hunter (25) have shown that the information matrix for multiple responses is given by M F (θ, ξ ) =
r r
σ ab Fa T Fb ,
(25)
a=1 b=1
where Fa =
(a) (a) ∂ηu (θ, t1 ) ∂ηu (θ, tn ) ,..., , ∂θ ∂θ
a = 1, . . . , r.
(26)
The D-optimal design then maximises the determinant of M F (θ, ξ ) as given in Eq. (25). It follows that if the information matrix in Eq. (25) is simplified to M F (θ, ξ ) =
r r
Fa T Fb ,
(27)
a=1 b=1
i.e. if it is assumed that the σ ab are all equal, then this is equivalent to maximising |F T F|,
Optimal Design for Model Discrimination and Parameter Estimation
529
where F=
r
Fa ,
(28)
a=1
∂ ra=1 η(a) (θ, tn ) ∂ ra=1 η(a) (θ, t1 ) ,..., = . ∂θ ∂θ
(29)
This will greatly simplify computation, as all that is required for the criterion is to sum the r responses in the model, and then calculate the Fisher information matrix as in the single response case. While this assumption may not be realistic for this study, to best of the authors’ knowledge there are no published values for the covariance of itraconazole and hydroxyitraconazole, so this algebraically convenient (and hence computationally simple) simplification is reasonable in these circumstances. This idea is extended to the more complex nonlinear mixed effects models with heteroscedastic error. To accomplish this, the Fisher information matrix described in Eqs. (17–22) is computed, replacing η(β, 0, ξ ) with ra=1 η(a) (β, 0, ξ ). In the itraconazole example, the only responses of interest are the concentrations in the 2nd and 4th compartments. Hence the sum is only over a = 2, 4. Model Discrimination There are two competing structural models to describe the pharmacokinetics of itraconazole. The first involves linear elimination of the parent drug, the second involves a nonlinear Michaelis–Menten process. Atkinson and Fedorov (26) introduced the T -optimality criterion for model discrimination, which could be applicable here. However, T -optimality is quite expensive computationally, and is not efficient for parameter estimation (see Waterhouse et al. (27) for examples). A viable alternative is to use the product of the determinants of the information matrices, raised to the power of the reciprocal of the number of population parameters (or equivalently, the product of efficiencies, as suggested by Atkinson and Cox (28)). Waterhouse et al. (27) also suggest that in some circumstances this method can produce designs which are useful for both goals: parameter estimation and model discrimination. This is demonstrated by the simulations in Section ‘Design Evaluation’. Combining Capsule and Solution Responses As mentioned in the section Structural Models there are four models under consideration. We are, however, only interested in discriminating
530
Waterhouse, Redmann, Duffull, and Eccleston
between the linear and nonlinear structural models. To do this, a capsule dose is administered, followed by a solution dose at a time after which it is assumed that there is no residual amount of the capsule dose. It is also assumed that there is no effect of the order of the doses. An elementary design ξ includes sampling times for the capsule (ξ c ) and solution (ξ s ). The responses of subject i to both doses are combined to give the overall response. So the parameters for capsules and solutions for a given structural model must also be combined, i.e. for u = 1 and 2 (linear and nonlinear models, respectively), θ u = ((θ u,c )T , (θ u,s )T )T . In the context of nonlinear mixed effects models with additive random effects, we now have the vector of parameters in model u for subject i given by θiu = ((θiu,c )T , (θiu,s )T )T , = =
((β + biu,c )T , (β u,s β u + biu , u,c
(30) + biu,s )T )T ,
(31) (32)
where β u = ((β u,c )T , (β u,s )T )T and biu = ((biu,c )T , (biu,s )T )T , with the random effects for dose form f ( f = c, s) for each subject being normally distributed with zero mean and with covariance matrix u, f = u, f u, f diag(ω1 , . . . , ω pu, f ). The combined covariance matrix for model u is then u = diag(u,c , u,s ). Equations (31) and (32) may be adjusted for exponential random effects by replacing β + bi with β exp(bi ) (for any superscripts on β and bi ). So the entire vector of population parameters under model u is then u = (( u,c )T , ( u,s )T )T
(33)
with u, f
u, f
u, f
u, f
u, f = ((β u, f )T , ω1 , . . . , ω pu, f , σ A , σ P )T .
(34)
In this notation, the combined ath response (for compartment a, where a = 2 or 4) for capsule and solution can be written ηu(a) (β u , biu , ξi ) = (ηu(a) (β u,c , biu,c , ξic )T , ηu(a) (β u,s , biu,s , ξis )T )T .
(35)
Final Criterion Hence the overall product criterion to be maximised is given by D P ( 1 , 2 , ) = |M F1 ( 1 , )|1/ p1 |M F2 ( 2 , )|1/ p2 ,
(36)
Optimal Design for Model Discrimination and Parameter Estimation
531
where M Fu ( u , ) =
Q
Nq M Fu ( u , ξq ),
u = 1, 2.
(37)
q=1
u is the vector of population parameters for the uth model. For the uth model, the elementary information matrix is given by 1 A(E u , S u ) C(E u , S u ) u u , (38) M F ( , ξq ) ≈ 2 C T (E u , S u ) B(E u , S u ) with A(·, ·), B(·, ·) and C(·, ·) as described in Eqs. (18–19), and E u = E u (y) ≈ ηu(a) (β u , 0, ξq ), a=2,4
S = S (y) ≈ u
u
(a) ∂ a=2,4 ηu (β u , 0, ξq )
+ diag σa + σ p
∂b T
u
(39)
(a) ∂ a=2,4 ηu (β u , 0, ξq ) ∂b T
ηu(a) (β u , 0, ξq ) ,
T
(40)
a=2,4
where the sums are over the responses for compartments 2 and 4.
METHODS In this section we outline the constraints of the sampling design, procedures used in selecting appropriate parameter and constant values, a number of simplifications to the models (to fit design constraints), and a brief overview of computational methods used to optimise the design criterion. Design Constraints A number of limitations on the experimental design were imposed for reasons of time and expense. A maximum of 30 patients were to be enrolled and studied on two occasions. The number of blood samples was also constrained to be four samples per occasion. A single 200 mg dose of the drug in capsules was to be given to each patient on the first occasion, followed by a single 200 mg dose in solution on the second occasion. A time constraint of approximately 2 months was placed on the design of the experiment, from exploration of the literature to computation of the final design.
532
Waterhouse, Redmann, Duffull, and Eccleston
Patients Thirty adults with CF who are admitted to hospital for treatment of a chest exacerbation will receive two doses of itraconazole (one 200 mg dose in capsules and one 200 mg dose as an oral solution) for this study. The participating patients are divided randomly into three groups, with each group being assigned one of three elementary sampling designs. Accordingly patients will be required to take two Spor onax Capsules as one dose (200 mg itraconazole), which is a standard adult dose. Four blood samples (0.5 ml) will be then drawn by finger prick according to the optimum sampling times of the design. After a wash-out period of at least three days (after the initial dose), the same dose of 200 mg (20 ml of 100 mg/10 ml solution) itraconazole would be taken as Spor onax Oral Solution. Four blood samples (0.5 ml will again be drawn by finger prick according to the optimum sampling times of the design. Simplification of Models The design optimisation problem under consideration is quite computationally intensive. Any justifiable reduction in the number of parameters to be estimated (which leads to a reduction in the dimension of the Fisher information matrix) is greatly beneficial. We believe that the different forms of dose will effect the input process only. So rather than treat the model parameters for capsule and solution doses as being entirely separate, as in Eqs. (10–13), we suppose that there are only two sets of parameters to be estimated, the first for the linear model and the second for the nonlinear model: β 1 = (ka1 , V21 , V31 , V41 , C L 14 , C L 124 , Q 1 , C L 12 )T ,
(41)
2 2 T , Km ) . β 2 = (ka2 , V22 , V32 , V42 , C L 24 , C L 224 , Q 2 , Vmax
(42)
The parameters kau represent the absorption rate for capsules (i.e. they are strictly kau,c ). The distinction between absorption rates is instead addressed by a multiplicative factor Fka = 0.9452/0.4159 ≈ 2.273, the ratio of the absorption rate for solutions and the absorption rate for capsules (not to be estimated). So the response for the solution under model u is given by using the absorption rate Fka kau . The smaller of the two absorption rates is used for this design, as it will be the most difficult to estimate. In addition to the computational intensity of the problem, the design constraints (outlined in Section ‘Design Constraints’) mean that not all of the parameters (of which there is still a large number) will be estimable.
Optimal Design for Model Discrimination and Parameter Estimation
533
Some preliminary investigation of expected standard errors of the estimates (given by the square roots of the diagonals of the inverse of the Fisher information matrices) show that, under both linear and nonlinear models, the parameters C L 24 , Q and V3 will have very imprecise estimates under the constraints of the sampling design. Further to this, K m is not estimable with any precision under the nonlinear model. These four parameters are therefore assumed to be fixed for the purposes of this design. They are included in the calculation of the Fisher information matrices, but their corresponding rows and columns are deleted for the calculation of the optimality criterion. This method produces subtly different information matrices to those calculated by leaving the parameters out of the calculation altogether. The differences in the methods are due to the approximations used in calculating the matrices. For a ‘true’ information matrix, the inclusion of other parameters should not effect the ∂ 2 l( ;y) existing elements of the matrix. For example, E − ∂β ∂β T should not i j be effected by whether the derivatives with respect to βk are included elsewhere in the matrix. In retrospect, these fixed parameters should have been left out of the calculation of the information matrix, but as the difference between the two methods is minimal, this was not considered a major issue. To further simplify computation, the residual amounts of drug and metabolite in each compartment at the time of the second dose was assumed to be insignificant. This was found to have a negligible effect on the criterion.
Parameter and Constant Values Finding the D-optimal design for any model which is nonlinear in its parameters requires prior specification of the model parameters. Such designs are known as ‘locally’ D-optimal. The parameter values (and values of other constants) used in the construction of our design were taken from a variety of sources. To find the local optimum design, parameters which are common to both linear and nonlinear elimination models (e.g. V2 , C L 4 ) are given the same value for both model types, despite our treatment of the parameters as separate entities when constructing the design criterion. For example, V21 and V22 are not the same population parameters, but for convenience we assign them the same value in the construction of the optimum design. A recently published population pharmacokinetic study by Koks et al. (29) gave us most of the parameter values we needed. In particular, the clearance of itraconazole from the central compartment and its
534
Waterhouse, Redmann, Duffull, and Eccleston
volume of distribution in the central compartment, and therefore the elimination rate constant from the central compartment (k20 ), but also the distribution rate constants (k23 , k32 ) and the formation rate constant (k24 ). These authors did not report any nonlinear elimination of the drug. But in several other pharmacokinetic studies nonlinear pharmacokinetic behaviour was shown for itraconazole (9,14,15,19). We consequently collected information about the Michaelis–Menten kinetic parameters (km and Vmax ) from other sources (14). Since we are planning to study the relative bioavailability of the oral solution compared to the capsule formation of itraconazole in patients with CF, we had to estimate the different absorbtion rate constants (kas and kac ) as well. This was done using ka =
loge 2 t1/2,
absorbtion
,
(43)
where t1/2, absorbtion is estimated by tmax /3. The times of peak plasma concentration for itraconazole in the oral solution and the capsule formulation from the drug information web page of the Food and Drug Administration (FDA) (30) were used. As no further information on PK parameters was available, the remaining model parameters were estimated by ‘best guess’. Table I shows, along with all previously mentioned parameters and constants, the values chosen for the elimination rate constant, the volume of distribution, and the clearance of the metabolite (k40 , V4 , and C L 4 , respectively). The best guess values of these parameters were chosen as those that produced a simulated concentration-time curve that was closely comparable to profiles provided by the authors of other clinical pharmacokinetic studies (9,32). Deterministic simulations were performed using MATLAB (these are shown in Fig. 3). Residual variability was assumed to have both additive and proportional components, with σ A = 5 × 10−3 mg/l (half of the lowest concentration detectible by the assay) and σ P = 0.15. Between subject variability was assumed to be lognormal with variance 0.1, i.e. θi = β exp(bi ) with bi ∼ N (0, 0.1I p ) for all parameters, where I p is the p × p identity matrix. Computational Methods A set of MATLAB routines, POPT Version 2.0, was used to calculate the optimal population design, optimising the criterion given in Section ‘Final Criterion’. This software allows the user to find D-optimal
Optimal Design for Model Discrimination and Parameter Estimation
535
Table I. Parameter and Constant Values Used to Calculated the Optimal Design Variable kas kac Fka k20
= = = =
k23 k32 k24 k40 C L2 C L4 Q Vmax Km V2 V4 V3 C L 24
= = = = = = = = = = = = =
F12 F14 F24 F20
= = = =
absorption rate constant, solution absorption rate constant, capsule absorption rate multiplicative factor elimination rate constant from the central compartment, linear model distribution rate constant distribution rate constant formation rate constant elimination rate constant for hydroxy-itraconazole clearance of itraconazole clearance of the metabolite inter-compartmental clearance theoretical maximum rate of the process Michaelis–Menten constant volume of central compartment volume of metabolite compartment volume in peripheral compartment clearance of itraconazole by metabolism to hydroxyitraconazole fraction of parent in the gut absolute bioavailability fraction of metabolite after first pass metabolism fraction of parent converted to metabolite fraction of parent eliminated
Value
Source
0.945 h−1 0.416 h−1 2.27 7.6 × 10−2 h−1
(A) (A) = kas /kac (B)
0.126 h−1 0.15 h−1 2.93 × 10−3 h−1 7.6 × 10−2 h−1 27.9 l/h 1.75 l/h 46 l/h 9.54 ng/(ml h) 329 ng/ml 365 l 23 l 307 l 1.07 l/h
(B) (B) (B) = k20 (B) = k40 V4 = k23 V2 (C) (C) (B) 1 V ≈ 16 2 = Q/k32 = k24 V2
0.55 4.3 × 10−2 0.5 0.5
(A) (D) (E) = 1 − F24
References: (A) = Janssen Pharmaceutica Products, L.P. (30); (B) = Koks et al. (29); (C) = Barone et al. (14); (D) = Janssen Pharmaceutica Products, L.P. (31); (E) = best guess derived from MATLAB simulation.
sampling times by specifying the number of blood samples, the upper and lower boundaries for the samples, the dose interval for taking blood samples, and the number of groups in the study. A paper by Duffull et al., Consideration of designs for population pharmacokinetic studies (manuscript submitted for publication), gives further description of the software used. POPT was modified to incorporate the product criterion and numerical solutions to the differential equations. Simulated Annealing POPT employs an algorithm based on simulated annealing for continuous variables to maximise the criterion in each case. Our implementation is based on the algorithm of Corana et al. (33) which applies simulated annealing to optimisation over continuous spaces. Also see
536
Waterhouse, Redmann, Duffull, and Eccleston 500
500
Concentration (ng/mL plasma)
(a)
(b)
400
400
300
300
200
200
100
100
0
0
12
24 Time (hr)
36
48
0
0
12
24
36
48
Time (hr)
Fig. 3. Concentration–time profiles of itraconazole (solid line) and hydroxyitraconazole (dashed line) after a 200 mg solution dose of itraconazole, with (a) linear elimination of the drug; and (b) nonlinear elimination of the drug. The values of the various parameters and constants used to produce the plots are given in Table I.
Goffe et al. (34) for examples of its use. We give a brief overview of the procedure here. Consider the design as an n-vector with elements i . At the kth iteration of the algorithm, each element of (i.e. each of the sampling times) is perturbed in sequence. The sizes of the perturbations are defined by a vector s = {si } of step sizes. The new value of the design variable i is taken from a uniform distribution on the interval i ± v i . The v i are adjusted at specified intervals so that, on average, around half of the new designs are rejected. For each new design, its criterion D P ( 1 , 2 , ) (described in Eq. (36)) is calculated and compared to the criterion for the current ‘best’ design ∗ . Designs with higher criteria are always accepted, designs with lower criteria are accepted or rejected according to the Metropolis criterion, i.e. the acceptance probability is
Optimal Design for Model Discrimination and Parameter Estimation
D P ( 1 , 2 , ) − D P ( 1 , 2 , ∗ ) exp Tk
537
(44)
for the current temperature Tk . The temperature is lowered at specified intervals by geometric cooling, Tk+1 = αTk (α = 0.75 was used here to speed up run times). The initial temperature was set at 1000. The algorithm stops when the average step size has reached a suitable value, i.e. the stopping criterion is i si /n < , for a pre-specified small tolerance . This differs slightly from the algorithm of Corana et al. (33), who instead rely on the criterion reaching a stable value, regardless of the step sizes. Sufficiently small step sizes ensure that the criterion does in fact reach a stable value. Numerical Solutions to Differential Equations In the case of the nonlinear model, there does not exist a closed form solution to Eqs. (5–8), so we cannot analytically isolate the responses of interest, A2 and A4 . Instead numerical solutions to the equations are used. For this purpose, we used MATLAB’s built-in ODE solver ode45, which is applicable to nonstiff problems and has a medium order of accuracy. A relative tolerance (RelTol) of 10−4 was used in these calculations. Numerical Derivatives Once the solutions to Eqs. (5–8) are found, we require their partial derivatives with respect to each parameter (as in Eq. (40)) . These derivatives were approximated numerically using the central difference method, ∂ f (x) f (x(1 + h 1 ) + h 2 ) − f (x(1 − h 1 ) − h 2 ) ≈ ∂x 2(xh 1 + h 2 )
(45)
for a suitably small h 1 and h 2 . This method is preferred over forward and backward difference methods as the error associated with the central difference method is proportional to the square of the step size (xh 1 + h 2 ) rather than the step size itself. This leads to more stable approximations to the derivatives. Since we are using numerical approximations to both the ODE solutions and the partial derivatives, it is important that the ODE solutions are more precise than the numerical derivatives. For example, if the ODE solutions are correct to the 4th decimal place, we must choose h 1 and h 2 such that the changes in f (·) must be larger than 10−4 , otherwise the results may be inaccurate. An alternative solution would be to use the ‘direct method’ (see Atkinson and Bogacka (35) for an example of its use).
538
Waterhouse, Redmann, Duffull, and Eccleston
RESULTS Optimal Design The set of optimal sampling times ∗ are given in Table II. The times have been rounded to the nearest minute. As the sampling times have been chosen from a continuous interval, these schedules are difficult to implement in practice. Indeed, in some cases they are actually impossible, with several blood samples being taken concurrently. It is therefore necessary to extend the population design to a sequence of intervals, or windows, in which the samples may be taken. We chose sampling windows of 3 hr (smaller for the earlier sampling times), roughly centered around the optimal times. These windows are also shown in Table II. The efficiency of a population design in this case is given by D P ( 1 , 2 , ) , D P ( 1 , 2 , ∗ )
(46)
where ∗ is the optimal population design. Any deviation from the optimal sampling times will obviously lead to a decrease in efficiency. The sampling windows were evaluated by taking 1000 designs at random from all the given intervals and calculating their efficiencies. The 1000 designs were Table II. Optimal Population Design ∗ for the Population Pharmacokinetics Study of Itraconazole with 10 Patients in Each of the Three Elementary Design Groups Capsule
Group (q)
1
2
3
Nq
10
10
10
Solution
Elementary Design ξqc (hrs:mins)
Sampling Window (hrs)
Elementary Design ξqs (hrs:mins)
Sampling Window (hrs)
1:14 8:56 25:49 51:45
0.1 → 3.0 7.0 → 10.0 24.0 → 27.0 50.0 → 53.0
0:17 3:55 3:56 3:56
0.1 → 1.0 3.0 → 3.5 3.5 → 4.0 4.0 → 4.5
6:13 9:50 29:29 29:29
50.0 → 8.0 8.0 → 11.0 28.0 → 29.5 29.5 → 31.0
0:18 4:06 4:06 72:00
0.1 → 1.0 3.0 → 4.0 4.0 → 5.0 69.0 → 72.0
8:08 28:00 72:00 72:00
7.0 → 10.0 26.5 → 29.5 69.0 → 70.5 70.5 → 72.0
0:17 4:22 27:08 72:00
0.1 → 1.0 3.0 → 6.0 26.0 → 29.0 69.0 → 72.0
Optimal Design for Model Discrimination and Parameter Estimation
539
350
300
250
200
150
100
50
0 0.85
0.9
0.95
1
Fig. 4. Histogram of efficiencies of 1000 simulations of sampling times taken uniformly within all of the sampling windows described in Table II.
selected by taking a set of sampling times from uniform distributions on all of the sampling windows, i.e. these are joint (not marginal) sampling windows. Figure 4 shows the distribution of these efficiencies. The average efficiency was about 0.95, with 96.2% of the designs having an efficiency of 0.90 or greater. The ‘worst’ design found had an efficiency of about 0.85. The sampling windows were deemed acceptable in light of these results, and did not need to be shortened to give any further increase in efficiency. Design Evaluation Model Discrimination The optimal design was evaluated in terms of its ability to select whether the ‘true’ underlying model involves linear or nonlinear elimination from the parent compartment. To do this, we simulated 100 sets of data under each model (linear elimination and nonlinear elimination) using the parameter values given in Table I and the optimal design provided in Table II. The models were fitted to each set of data using NONMEM (version 5 with the G77 FORTRAN
540
Waterhouse, Redmann, Duffull, and Eccleston
compiler, using FOCE & INTERACTION), and the ‘best’ model was chosen empirically as the one with the smallest objective function (which is proportional to minus twice the log-likelihood). Recall that C L 24 , Q and V3 are considered to be fixed under both the linear and nonlinear models, and that K m is considered fixed in the nonlinear elimination model. After simulating under the linear model, the linear model was chosen correctly 74% of the time. After simulating under the nonlinear model, the nonlinear model was chosen correctly 100% of the time. For all estimation runs, the NONMEM run was reported to have converged successfully. These results are considered to be reasonable, especially considering that at low concentrations, the difference between linear and Michaelis–Menten elimination is minimal. Standard Errors The results from these simulations were also used to calculate standard errors of the parameters estimates. Empirical estimates of the standard errors were obtained from the standard deviation of the estimated parameter values from the 100 simulated data sets. This was performed only in those cases where the same model was used for simulation and estimation. These are given in Table III for both the linear and nonlinear models. We can compare these standard errors to those predicted by the information matrix used to calculate the optimal design. These ‘expected’ standard errors are given by the square roots of the diagonals of the inverse of each Fisher information matrix. These are given in Table III. The comparison of expected and simulated standard errors reveal quite satisfactory results, considering the difficulty in estimating between subject variability and absorption rate constants. The only parameter with a significant difference between the expected and estimated standard errors is ka . The absorption rate constant is notoriously difficult to estimate in PK models, so it is not unexpected that the estimates from the simulations were so variable. Parameter Estimation To assess the benefit of the optimal population design with respect to parameter estimation, two additional designs were calculated: one is optimal in terms of the linear model alone, and the other optimal in terms of the nonlinear model alone. The optimal design for model u, denoted by
∗u , was found by maximising the criterion |M Fu ( u , )|1/ pu ,
(47)
using the same methods as described in Section ‘Computational Methods’.
Optimal Design for Model Discrimination and Parameter Estimation
541
Table III. Standard Errors (Expressed as Coefficients of Variation) for the Population Mean and Between Subject Variability for Parameters Estimated Under the Linear and Nonlinear Models Model parameter Linear elimination of parent Population mean Expected CV (%) Estimated CV (%) Between subject variability Expected CV (%) Estimated CV (%) Nonlinear elimination of parent Population mean Expected CV (%) Estimated CV (%) Between subject variability Expected CV (%) Estimated CV (%)
ka
V2
V4
C L4
C L2
8.39 54.34
16.39 18.71
10.64 18.20
10.65 10.76
8.57 8.32
41.22 121.00
82.17 114.61
51.57 80.85
45.58 68.30
39.95 89.80
ka
V2
V4
C L4
Vmax
8.19 25.25
15.55 10.30
11.26 9.85
12.51 10.58
27.84 33.44
40.81 111.00
49.04 76.10
49.28 83.00
61.62 81.10
197.73
The expected standard errors are calculated from the Fisher information matrix, the estimated standard errors are the standard deviations of parameter estimates from 100 NONMEM simulation-estimation runs. The dash indicates that normality was rejected by the Kolmogorov– Smirnov test (the data was highly positively skewed), in which case calculation of standard errors is meaningless.
The efficiency of any population design with respect to model u is defined (in a similar fashion to Eq. (46)) to be the ratio Effu ( , ∗u ) =
|M Fu ( u , )|1/ pu , |M Fu ( u , ∗u )|1/ pu
(48)
where ∗ is the product optimal design given in Table II. For the product optimal design ∗ , the criteria given by Eq. (47) were calculated under the linear and nonlinear models separately. The values found were 37.591 and 115.1, respectively. The criteria for designs ∗1 and
∗2 , which are optimal under the linear and nonlinear model, respectively, are 38.993 and 120.01. Hence the efficiencies of ∗ in terms of the linear and nonlinear models are: Eff1 ( ∗ , ∗1 ) = and
37.591 = 0.9640 38.993
542
Waterhouse, Redmann, Duffull, and Eccleston
Eff2 ( ∗ , ∗2 ) =
115.1 = 0.9591. 120.01
This indicates that the product optimal design is suitably efficient at estimating the parameters of both models.
DISCUSSION The practical example of the itraconazole study is still underway at the time of writing. The authors are currently waiting on results from the actual trial. For comparison, an empirical ‘rich’ sampling design, in which 14 samples are taken over almost the entire dosing interval (72 hr) is given by (0.01, 0.5, 1, 1.5, 2, 4, 6, 8, 12, 18, 24, 36, 48, 70). In this case, again in a study involving thirty patients, all parameters in question would be estimable. The relative standard error of all fixed effects would all be less than 20%, in most cases less than 10%. The construction of the optimal design was subject to several constraints. A timeline of approximately 2 months was given to develop the entire design, from collating prior information on itraconazole and researching applicable methods, to computing the final design (a procedure which took over 7 days in itself). Future research in this area includes the consideration of further methodology (such as the Direct Method (35)). Under the given time constraints, computation time was a major issue, with the final design computation taking over a week to run (on a 2.4 GHz Pentium 4 with 1 GB of RAM). Cost limitations of the study and limitations imposed by the relevant ethics committee also limited the possible number of patients and the number of samples per patient, which further constrained the design in terms of the number of parameters which are estimable. Several ideas were used in this design which are novel to the field of experimental design for population pharmacokinetic studies. These include the use of the D-optimality product criterion for model discrimination and the incorporation of multiple responses into the Fisher information matrix (a problem previously considered in different contexts (36,37)). The results of the power tests and simulations show that the use of this product criterion has created a design which is efficient in terms of both parameter estimation and model discrimination. Whether the individual D-optimal designs were also adequate for model discrimination was not investigated for this design. However, we have seen evidence (27) that suggests that this is not true for nonlinear models.
Optimal Design for Model Discrimination and Parameter Estimation
543
The standard errors predicted by the inverse of the Fisher information matrix generally match those estimated by the simulations. We believe that these results show that the errors introduced by the approximations made by Retout and Mentr´e (22) are mostly inconsequential. Also, the assumptions about the covariance structure of the two responses in the model seem to be justified by these results. The sparse sampling design presented in this paper gives a simple and cost-effective data collection regimen which should provide acceptable estimates of the parameters of interest.
ACKNOWLEDGMENTS The author acknowledge helpful discussions with Ivelina Gueorguieva, Andrew Hooker, and France Mentr´e.
REFERENCES 1. R. B. Moss. Allergic bronchopulmonary aspergillosis. Clin. Rev. Allergy Immunol. 23:87–104 (2002). 2. M. Skov, N. Hoiby, and C. Koch. Itraconazole treatment of allergic bronchopulmonary aspergillosis in patients with cystic fibrosis. Allergy 57:723–728 (2002). 3. D. A. Stevens, J. Y. Lee, H. J. Schwartz, D. Jerome, and A. Catanzaro. A randomized trial of itraconazole in allergic bronchopulmonary aspergillosis. N. Engl. J. Med. 342:756–762 (2000). 4. A. Hedman, Y. Adan-Abdi, G. Alvan, B. Strandvik, and A. Arvidsson. Influence of the glomerular filtration rate on renal clearance of ceftazidime in cystic fibrosis. Clin. Pharmacokinet. 15:57–65 (1988). 5. W. J. Jusko, L. L. Mosovich, L. M. Gerbracht, M. E. Matter, and S. J. Yaffe. Enhanced renal excretion of dicloxacillin in patients with cystic fibrosis. Pediatrics 56:1038–1044 (1975). 6. E. Finkelstein and K. Hall. Aminoglycosides clearance in patients with cystic fibrosis. J. Pediatr. 94:163–164 (1979). 7. M. D. Reed. The pathophysiology and treatment of cystic fibrosis. J. Pediatr. Pharm. Pract. 2:285–308 (1997). 8. M. Spino. Pharmacokinetics of drugs in cystic fibrosis. Clin. Rev. Allergy 9:169–210 (1991). 9. J. Heykants, A. Van Peer, V. Van de Velde, P. Van Rooy, W. Meuldermans, and K. Lavrijsen. The clinical pharmacokinetics of itraconazole: an overview. Mycoses 32:67–87 (1989). 10. C. Koks, R. Sparidans, G. Lucassen, K. Crommentuyn, and J. Beijnen. Selective high performance liquid chromatography assay for itraconazole and hydroxyitraconazole in plasma from human immunodeficiency virus-infected patients. J. Chromatogr. B 767:103–110 (2002). 11. J. S. Hostetler, J. Heykants, K. Clemons, R. Woestenborghs, L. H. Hanson, and D. A. Stevens. Discrepancies in bioassay and chromatography determination explained by metabolism of itraconazole to hydroxyitraconazole: studies of interpatient variations in concentrations. Antimicrob. Agents Chemother. 37:2224–2227 (1993). 12. J. Van Cutsem. The in-vitro antifungal spectrum of itraconazole. Mycoses 32:7–13 (1989).
544
Waterhouse, Redmann, Duffull, and Eccleston
13. S. Jaruratanasirikul and A. Kleepkaew. Influence of an acidic beverage (coca-cola) on the absorption of itraconazole. Eur. J. Clin. Pharmacol. 52:235–237 (1997). 14. J. A. Barone, J. Koh, R. H. Bierman, J. L. Colaizzi, K. Swanson, M. Gaffar, B. L. Moskovitz, W. Mechlinski, and V. Van de Velde. Food interaction and steady-state pharmacokinetics of itraconazole capsules in healthy male volunteers. Antirnicrob. Agents Chemother. 37:778–784 (1993). 15. V. Van de Velde, A. Van Peer, J. Heykants, R. Woestenborghs, P. Van Rooy, K. De Beule, and G. Cauwenbergh. Effect of food on the pharmacokinetics of a new hydroxypropyl, β-cyclodextrin formulation of itraconazole. Pharmacotherapy 16:424–428 (1996). 16. S. Grant and S. Clissold. Itraconazole: a review of its pharmacodynamics and pharmacokinetic properties and therapeutic use in superficial and systemic mycoses. Drugs 37:310–344 (1989). 17. K. De Beule and J. Van Gestel. Pharmacology of itraconazole. Drugs 61:27–37 (2001). 18. T. Zimmerman, R. Yeates, H. Laufen, G. Pfaff, and A. Wildfeuer. Influence of concomitant food intake on the oral absorption of two triazole antifungal agents, itraconazole and iuconazole. Eur. J. Clin. Pharmacol. 46:147–150 (1994). 19. A. Van Peer, R. Woestenborghs, J. Heykants, R. Gasparini, and G. Gauwenbergh. The effects of food and dose on the oral systemic availability of itraconazole in healthy subjects. Eur. J. Clin. Pharmacol. 36:423–426 (1989). 20. J. A. Barone, B. L. Moskovitz, J. Guarnieri, A. Hassell, J. L. Colaizzi, R. H. Bierman, and L. Jessen. Enhanced bioavailability of itraconazole in hydroxypropyl-β-cyclodextrin solution versus capsules in healthy volunteers. Antimicrob. Agents Chemother. 42:1862–1865 (1998). 21. T. C. Hardin, J. R. Graybill, R. J. Fetchick, R. Woestenborghs, M. G. Rinaldi, and J. Kuhn. Pharmacokinetics of itraconazole following oral administration to normal volunteers. Antimicrob. Agents Chemother. 32:1310–1313 (1988). 22. S. Retout and F. Mentr´e. Further developments of the fisher information matrix in nonlinear mixed effects models with evaluation in population pharmacokinetics. J. Biopharm. Stat. 13:209–227 (2003). 23. A. C. Atkinson and A. N. Donev. Optimum Experimental Designs. Oxford University Press, Oxford, 1992. 24. S. Retout, F. Mentr´e, and R. Bruno. Fisher information matrix for non-linear mixedeffects models: evaluation and application for optimal design of enoxaparin population pharmacokinetics. Stat. Med. 21:2623–2639 (2002). 25. N. R. Draper and W. G. Hunter. Design of experiments for parameter estimation in multiresponse situations. Biometrika 53:525–533 (1966). 26. A. C. Atkinson and V. V. Fedorov. The design of experiments for discriminating between two rival models. Biometrika 62:57–70 (1975). 27. T. H. Waterhouse, J. A. Eccleston, and S. B. Duffull. On optimal design for discrimination and estimation. Research Report 106, Centre for Statistics, The University of Queensland (2004). 28. A. C. Atkinson and D. R. Cox. Planning experiments for discriminating between models. J. R. Stat. Soc. Ser. B 36:321–348 (1974). With discussion by H. P. Wynn, D. M. Titterington, P. J. Laycock, D. V. Lindley, D. H. Hill, Agnes M. Herzberg, P. A. Tukey, A. O’Hagan, V. V. Fedorov, J. Dickey, J. Kiefer and C. A. B. Smith. 29. C. H. W. Koks, A. D. R. Huitema, E. D. M. Kroon, T. Chuenyam, R. W. Sparidans, J. M. A. Lange, and J. H. Beijnen. Population pharmacokinetics of itraconazole in thai hiv-1-infected persons. Ther. Drug Monit. 25:229–233 (2003). 30. Janssen Pharmaceutica Products, L.P. Sporanox (itraconazole) oral solution (2002). URL http://www.fda.gov/cder/foi/label/2002/20657s7lbl.pdf. 31. Janssen Pharmaceutica Products, L.P. Sporanox (itraconazole) injection (2002). URL http://www.fda.gov/cder/foi/label/2002/20966s61bl.pdf.
Optimal Design for Model Discrimination and Parameter Estimation
545
32. J. A. Barone, B. L. Moskovitz, J. Guarnieri, A. Hassel, J. L. Colaizzi, R. H. Bierman, and L. Jessen. Food interactions and steady-state pharmacokinetics of itraconazole oral solution in healthy volunteers. Pharmacotherapy 18:295–301 (1998). 33. A. Corana, M. Marchesi, C. Martini, and S. Ridella. Minimizing multimodal functions of continuous variables with the “simulated annealing” algorithm. ACM Trans. Math. Softw. 13:262–280 (1987). 34. W. L. Goffe, G. D. Ferrier, and J. Rogers. Global optimization of statistical functions with simulated annealing. J. Econ. 60:65–99 (1994). 35. A. C. Atkinson and B. Bogacka. Compound and other optimum designs for systems of nonlinear differential equations arising in chemical kinetics. Chemometr. Intell. Lab. 61:17–33 (2002). 36. I. Gueorguieva. D-optimal Design for Multivariate Response: Prospective Planning of a Beta-blocker Study in Rat Involving Cassette Dosing. Poster Presented at PAGE 2003. 37. A. Hooker. Simultaneous Population D-optimal Designs of PharmacokineticPharmacodynmic Experiments. Poster Presented at PAGE 2003.
