Abstract
Purpose: Doserelated toxicity of cyclophosphamide may be reduced and therapeutic efficacy may be improved by pharmacokinetic sampling and dose adjustment to achieve a target area under the curve (AUC) for two of its metabolites, hydroxycyclophosphamide (HCY) and carboxyethylphosphoramide mustard (CEPM). To facilitate realtime dose adjustment, we developed opensource code within the statistical software R that incorporates individual data into a population pharmacokinetic model.
Experimental Design: Dosage prediction performance was compared to that obtained with nonlinear mixedeffects modeling using NONMEM in 20 cancer patients receiving cyclophosphamide. Bayesian estimation of individual pharmacokinetic parameters was accomplished from limited (i.e., five samples over 016 hours) sampling of plasma HCY and CEPM after the initial cyclophosphamide dose. Conditional on individual pharmacokinetics, simulations of the AUC of both HCY and CEPM were provided for a range of second doses (i.e., 0100 mg/kg cyclophosphamide).
Results: The results compared favorably with NONMEM and returned accurate predictions for AUCs of HCY and CEPM with comparable mean absolute prediction error and root mean square prediction error. With our method, the mean absolute prediction error and root mean square prediction error of AUC CEPM were 11.0% and 12.8% and AUC HCY were 31.7% and 44.8%, respectively.
Conclusions: We developed dose adjustment software that potentially can be used to adjust cyclophosphamide dosing in a clinical setting, thus expanding the opportunity for pharmacokinetic individualization of cyclophosphamide. The software is simple to use (requiring no programming experience), reads individual patient data directly from an Excel spreadsheet, and runs in less than 5 minutes on a desktop PC.
 Cyclophosphamide
 population pharmacokinetics
 dose adjustment
The alkylating agent cyclophosphamide is used for treating a wide range of malignancies in both pediatric and adult patients and as part of myeloablative conditioning regimens before hematopoietic cell transplant. Cyclophosphamide dosing based on body weight (mg/kg) or body surface area (mg/m^{2}) results in highly variable area under the concentrationtime curve (AUC) of the parent drug between patients and an even greater variability in the AUC of its metabolites (1, 2).
Cyclophosphamide is a prodrug with complex metabolic features that include autoinduction and inhibition of its metabolism. Approximately 10% to 30% of a cyclophosphamide dose is excreted unchanged in the urine, whereas 5% of a dose is inactivated to dechloroethylcyclophosphamide. The remaining 65% to 85% of a dose is activated by several cytochrome P450 enzymes to hydroxycyclophosphamide (HCY); this oxidative pathway is autoinducible. HCY or its tautomer aldophosphamide is metabolized via at least five separate pathways, the most relevant of which to this work is its conversion to carboxyethylphosphoramide mustard (CEPM) by aldehyde dehydrogenase 1A1 (3–5). Increased CEPM systemic exposure, as measured by its AUC, has been associated with an increased risk of sinusoidal obstruction syndrome and increased mortality in patients undergoing cyclophosphamide/total body irradiation conditioning for hematopoietic cell transplant (2). In addition, a pharmacodynamic relationship has been found between cyclophosphamide and its metabolites and the risk of recurrence in children receiving 300 to 1,000 mg/m^{2} cyclophosphamide for treatment of Bcell nonHodgkin's lymphoma (6). A higher recurrence rate was found in children with low cyclophosphamide clearance (i.e., <3.5 L/h·m^{2}) and detectable plasma concentrations of two cyclophosphamide metabolites, CEPM and dechlorocyclophosphamide. Thus, dose adjustment of cyclophosphamide based on the AUC of its metabolites, specifically HCY and CEPM, may improve efficacy and/or reduce toxicity of cyclophosphamide.
Bayesian dose adjustment based on population pharmacokinetics has been deployed in various therapeutic contexts (7–9), including applications to anticancer drugs (10). In fact, several Bayesian models (11–14) have been developed for the pharmacokinetics of cyclophosphamide and its metabolites. In this context of dose individualization, “Bayesian” is taken to mean maximum a posteriori (MAP) estimation, where the population pharmacokinetics is represented by a nonlinear mixedeffects model and the individual pharmacokinetic parameters are estimated conditional on the population parameters. During a previous study of 147 patients who received a regimen of cyclophosphamide/total body irradiation in preparation for hematopoietic cell transplant (14), we developed a population pharmacokinetic model describing cyclophosphamide elimination and formation of HCY by noninducible and inducible routes, the latter of which was described by a hypothetical enzyme compartment (see Fig. 1 ). Prediction of CEPM AUC was clinically accurate and precise (mean absolute prediction error, −3.5%; root mean square prediction error, 12.2%) based on plasma concentration time data limited to five to six time points obtained over 16 hours after initialization of the first cyclophosphamide dose. However, the AUC of HCY was generally slightly overestimated (mean absolute prediction error, 16.923.6%; ref. 14).
Concurrent with this model development, a phase I trial was conducted in 20 adults receiving the cyclophosphamide/total body irradiation conditioning regimen for whom cyclophosphamide doses were adjusted in realtime via a regression model to achieve a target AUC of CEPM and HCY (15). Pharmacokinetic plasma samples were obtained over 16 hours after initiation of the first cyclophosphamide dose. Plasma CEPM and HCY concentrations were immediately assayed followed by noncompartmental analysis to estimate the optimal cyclophosphamide dose for achieving respective target AUC of CEPM and HCY of 325 ± 25 and >50 μmol/L·h. The target AUCs were based on a pharmacokinetic/pharmacodynamic analysis of data from a previous population study, in which CEPM AUC <325 μmol/L·h was associated with the lowest risk of nonrelapse mortality, mainly due to sinusoidal obstruction syndrome, and the minimal HCY AUC was 50 μmol/L·h (2). Although dose adjustment of cyclophosphamide was feasible using noncompartmental analysis, a retrospective analysis revealed that Bayesian modeling of cyclophosphamide metabolism (Fig. 1) using HCY and CEPM plasma concentrations from 0 to 16 hours after cyclophosphamide dose 1 could lead to more accurate dose tailoring (15).
Implementation of this integrated mechanismbased model for clinical dosage adjustment requires computational software that is operationally straightforward and efficient, such that the analysis can be accomplished in realtime by laboratory personnel within minutes after plasma CEPM and HCY data are available to achieve target AUCs. The complexity of the pharmacokinetic model and the necessity of rapid calculations motivated the development of a customized software tool to perform MAP Bayesian calculation of pharmacokinetic parameters. We wrote the computer code within the opensource statistical software R (16) using packages RODBC (17), MVTNORM (18), and ODEsolve (19), which allows for rapid dose adjustment of cyclophosphamide within a clinical setting. Our code does not require the availability of complex population pharmacokinetic modeling programs, such as NONMEM, reads individual patient data directly from an Excel spreadsheet, and can be distributed under an opensource license. To our knowledge, this is the first time that a MAP Bayesian dose adjustment software, developed as open source, is made available to the scientific and clinical community. This software can be adapted for application beyond cyclophosphamide dose adjustment and can expand the availability of Bayesian modeling techniques to individualize dosing of chemotherapy drugs in hopes of improving their therapeutic efficacy and/or reducing their doserelated toxicity. For the remainder of this article, we will refer to this novel method as BaRD or “BAyesian R Dose adjustment.”
Materials and Methods
Study population. Pharmacokinetic data were obtained from 20 patients enrolled in a previous phase I trial to evaluate the feasibility of realtime dose adjustment of cyclophosphamide in hematopoietic cell transplant patients undergoing a cyclophosphamide/total body irradiation preparative regimen (15). All patients were given an initial i.v. infusion of 45 mg cyclophosphamide per kg body weight followed by serial blood sampling over 16 hours. The second cyclophosphamide dose was administered at 24 hours after initiation of dose 1. During the original study, the second dose of cyclophosphamide was calculated based on a regression equation relating the measured CEPM and HCY AUCs from 0 to 16 hours (i.e., the pharmacokinetic sampling period) to that from 0 to 48 hours, the pharmacokinetic indicator of nonrelapse mortality risk (15). At the completion of patient accrual, our previous study (14) evaluated whether the Bayesian approach using plasma concentrations of CEPM and HCY only (without incorporation of cyclophosphamide concentrations due to time required for its quantitation) accurately estimated individual pharmacokinetic parameters and CEPM and HCY concentrations using the actual second cyclophosphamide dose as chosen during the protocol (i.e., by the regression equation). As we reported previously (15), this retrospective evaluation showed that the Bayesian approach implemented in NONMEM (University of California at San Francisco and Globomax ICON, Hanover, MD), a nonlinear mixedeffects regression modeling software, was more accurate in predicting the AUC of CEPM compared with the regression method and only marginally better for the prediction of HCY AUC. Moreover, plasma concentrationtime data measured between 16 and 48 hours did not inform the model, and were used for comparison between predicted and observed AUCs. To assess the performance of our new Rbased MAP Bayesian parameter estimation software, we compared its predictions for the 20 patients in the phase 1 trial to those obtained using NONMEM.
Pharmacokinetic model. The population pharmacokinetics of cyclophosphamide and its metabolites can be described with a fourcompartment model represented schematically in Fig. 1 and as a system of nonlinear differential equations:where A_{1}(t), A_{3}(t), and A_{4}(t) (mg/kg) are amounts of cyclophosphamide, HCY, and CEPM, respectively, and A_{2}(t) (dimensionless) represents the amount of ICY, the autoinducible enzyme pool for the 4hydroxylation of cyclophosphamide as proposed previously (20). ICY is a nonlinear function of plasma cyclophosphamide concentrations, which were estimated conditional on CEPM and HCY concentrations. CL_{NON} and CL_{IND} are noninducible and initial inducible clearances for cyclophosphamide, respectively (L/kg/h), whereas KL_{ENZ} is a rate parameter describing enzyme synthesis and elimination and E_{MAX} and EC_{50} are enzyme synthesis parameters. V_{CY} is the apparent volume of distribution of cyclophosphamide (liters) and K_{HCY} and K_{CEPM} are the elimination rate constants (in reciprocal hours) of HCY and CEPM, respectively. All initial conditions are zero, with the exception of the baseline enzyme amount A_{2}(0), which is scaled to be equal to 1.
Plasma concentrations are predicted by dividing the amounts for cyclophosphamide, HCY and CEPM, i.e., A_{1}(t), A_{3}(t) and A_{4}(t), by their respective apparent volumes of distribution. Model parameter V_{CY} represents the cyclophosphamide volume of distribution. The HCY and CEPM volumes (V_{HCY} and V_{CEPM}) were assumed to be equal at 0.013 L/kg. Further model details are available in our earlier publication (14). The modelpredicted concentrations for HCY and CEPM are given bywhere t_{1}, …, t_{n} are the available sampling times.
Our intent was to determine an individual parameterization of the pharmacokinetic model (Eq. A) given the fixed or population effects that describe both central parameter tendencies and their expected variation in the patient population and conditional only on parameter priors derived from the original 147patient group and limited HCY and CEPM plasma samples taken from the new patient within the first 16 hours after initiation of the cyclophosphamide infusion. This individualized pharmacokinetic model would then be used to determine the second dose for that patient to achieve the target 0 to 48hour HCY and CEPM AUCs (i.e., AUC CEPM 325 ± 25 μmol/L h and AUC HCY >50 μmol/L·h). One present challenge is that plasma cyclophosphamide measurements were not available during realtime dosage adjustment, whereas the original pharmacokinetic model was developed based on plasma concentration data for cyclophosphamide and its two metabolites.
Statistical framework. In general, given a statistical model for each data point over a time course, we assume an additive error scheme: y = f(b) + ε, where y = [y_{1}, y_{2}, …, y_{n}] is the time vector of observed values (data points at times t_{1}, …, t_{n}), f(b) is the time vector of corresponding modelpredicted values dependent on the unknown parameter vector b, and ε ∼ N (0,σ^{2}) is the corresponding vector of random measurement errors assumed to be Gaussian with zero mean. Then, the expected value of the observations with respect to the measurement error is given by E[y] = f(b) and
The joint likelihood of all available data given the parameters is the product of the individual likelihoods as follows:where n is the total number of data points. In multivariable fitting, such as when maximizing likelihood conditioned on both HCY and CEPM data, the vector y is assumed to contain all available values (this only causes the notation to become slightly more complicated, but the above rationale still holds with no modifications).
The above description assumes that the measurement error (and thus the data) is characterized by constant variance (homoscedastic measurements). If the measurement error variance changes instead with time and/or the measurement value, then the data are heteroscedastic. An example of heteroscedasticity is given by a statistical model comprising both a proportional and an additive (constant variance) component, such as y = (1 + ε_{2})f + ε_{1}, where ε_{1} ∼ N(0,σ_{1}^{2}) and ε_{2} ∼ N(0,σ_{2}^{2}). Then, E[y] = f(b) again, but
The negative log likelihood (nLL; or extended least squares estimator; ref. 21) from Eq. D can be expressed in vectormatrix notation and is given bywhere R(b) is a diagonal matrix of variances with elements like Eq. C corresponding to homoscedastic data and elements like Eq. E corresponding to heteroscedastic (proportional and additive) data.
In our case, the data vector y consists of measured plasma concentrations of both HCY and CEPM at times i = 1, …, n. Similarly, the vector f(b) consists of modeled plasma concentrations of both HCY and CEPM at times i = 1, …, n and is computed conditional on the pharmacokinetic model of cyclophosphamide and its metabolites (Eq. A).
Parameter prior information. MAP estimation provides a way to individualize pharmacokinetic parameters based on limited individual information and population information about the model parameters (22). If the model parameters b ∼ N(b̄;,ω^{2}) are independent and normally distributed, the likelihood of the parameters can be written aswhere b̄;_{j}, ω_{j}^{2}, the prior mean and variance of parameter b_{j}, are assumed uncorrelated, so ω^{2} is a diagonal matrix. In ref. 14, a lognormal distribution was instead assumed for all betweensubject variation in model parameters. That is, log b ∼ N(log b̄;,D) where diagonal matrix D is the covariance of the log b. The likelihood of the parameters can be written similarly as
Combining parameter and data likelihoods. Using Bayes' rule, we combine the likelihood of the data given the parameters (Eq. D) with the prior for the parameters (Eq. G) to obtain the joint likelihood of the parameters and data. The negative log of this joint likelihood is:where b is the vector of random effects to be determined by minimization of the objective (loss) function given by (Eq. H). Note that the rate parameters that were characterized by betweensubject variation (see Table 1 and Qiu et al. for more details) were CL_{NON}, CL_{IND}, K_{ENZ}, K_{HCY}, K_{CEPM}, and K_{34}. The other model parameters do not vary among subjects and are thus not individualized. b̄; is the vector of populationwide expected values (means; see Table 1) of the elements of the vector b that vary among subjects. y is the vector of measured HCY and CEPM concentrations in the plasma at specified times. f(b) is the vector of modeled HCY and CEPM concentrations (Eq. B) in the plasma at specified times (same as for y above). D is the prior covariance of the logarithm of the random effects; a diagonal matrix of known values (each corresponding to an element of b; see Table 1) determined by population pharmacokinetic analysis. Note that the elements of D have dimensions of coefficient of variation squared as the lognormal density is approximated (to the first order) by a proportional Gaussian. R(b) is the covariance of measurement error. Based on previous results (14), the measurement error of HCY was assumed to be heteroscedastic (a combined proportional and additive error), whereas that associated to CEPM was assumed to be homoscedastic (constant variance). Thus, R(b) is a diagonal matrix of the formwhere the first n elements of the diagonal provide the HCY concentration measurement variance (with λ_{1} = Var[ε_{1}] and λ_{2} = Var[ε_{2}] from Eq. E); the final n elements provide the CEPM concentration measurement variance (λ_{3} = Var[ε] from Eq. C); and the metabolite dynamics is defined by Eq. A in terms of metabolites A_{1}(t), …, A_{4}(t) representing concentrations of cyclophosphamide, ICY, HCY, and CEPM, respectively. Note that although only HCY and CEPM concentrations are measured, simultaneous modeling of all four metabolite time courses is necessary.
The diagonal nature of matrices R(b) and D allows us to rewrite the joint log likelihood (Eq. H) as a summation instead of a vectormatrix product, simplifying subsequent computations. As mentioned before, model parameter values can be partitioned among parameters that do and do not vary among subjects. More specifically, the values of b̄;, D, λ_{1}, λ_{2}, and λ_{3} as well as three model parameters with no associated random effect (V_{CY}, E_{MAX}, and EC_{50}) were taken from the previous population study (14). The elements of b̄; and the corresponding elements of D (the covariance of the logarithm of the elements of b) are listed in Table 1. Note also that the typical value of CL_{NON}, the noninducible clearance of cyclophosphamide, actually depends on the individual patient's age (in years) via a power function. The variance parameters in the measurement error of HCY and CEPM are also presented in Table 1.
Strategy. We individualized the model (Eq. A) by determining the parameters b, which minimize the joint nLL (Eq. H), where y equals the HCY and CEPM plasma concentration data from the first 16 hours after cyclophosphamide infusion initiation and f(b) provides the corresponding modeled concentrations. It is necessary to resolve the system of ordinary differential equations (Eq. A) in each iteration of the optimization method. This is accomplished using the function lsoda available in the ODEsolve package (19) in R. Individual simulations of plasma HCY and CEPM were then modeled over 48 hours with a range of second doses applied at 24 hours and using final metabolite amounts from the first 24 hours as initial conditions for the second, extrapolated, 24hour period. As a measure of metabolite exposure, AUCs with confidence intervals were computed for all second dose simulations.
Parameter optimization. Of the standard optimization methods available in R to minimize the value of the nLL with respect to the individual parameters, we chose the BroydenFletcherGoldfarbShanno method (a quasiNewton method that uses function and gradient values to build up a picture of the optimization surface) based on brief testing for both stability and accuracy (data not shown).
The populationwide parameter values were used as starting values for the individual optimization and the default stopping tolerance was employed. Because the routine is to be used for realtime dose adjustment, we included an emergency runtime option to employ a specifiable number of NelderMead method iterations for the purpose of improving starting values before employing the BroydenFletcherGoldfarbShanno method. This option proved useful when, to test the robustness of the code, optimization was initiated with unreasonable starting parameter estimates.
Early during parameter search, some negative values of the model parameters were attempted by the optimization routine. A smooth penalty formulation was instituted to discourage this unphysiologic behavior. Negative parameters are replaced by 10^{−18} in that iteration to safeguard the ordinary differential equation and adjoint ordinary differential equation solvers and the nLL function from unpredictable results and possible runtime crashes. A smooth penalty of the form 10^{5} b_{j}^{2}/b̄;_{j} was added to the nLL value for each negative parameter b_{j} in that iteration. It is noteworthy that the penalties were never in effect at convergence.
Derivative calculations. Optimization of individual parameters requires the calculation of the gradient of the nLL (Eq. H). To compute derivatives of the log likelihood, we need partial derivatives of f(b) with respect to all elements of b. Equivalently, we need partial derivatives of predicted amounts of HCY and CEPM, A_{3}(t) and A_{4}(t), with respect to parameter vector b. These are found as solutions to the adjoint ordinary differential equations, which are formulated by lettingwhere A(t,b) = [A_{1}(t)A_{2}(t)A_{3}(t)A_{4}(t)]. The new system of ordinary differential equations becomeswhere the function g is the righthand side of the original system of ordinary differential equations (Eq. A) andSolving the system of eight ordinary differential equations (Eq. J) simultaneously provides solution to the original system of equations (Eq. A) as well as the derivative of the model functions with respect to the parameter b_{j}:
Confidence intervals. To compute confidence intervals for the AUCs, let b̂ = b̃ + Δb, where b̂ are the parameter estimates, b̃ are the unknown actual parameter values, and Δb are the errors (discrepancies between true and estimated parameter values). Then, Var(b) = E[ΔbΔb^{T}] = H^{−1}, where H is the Hessian of the nLL evaluated at b̂. To compute H, we used centered differences on the entries of the gradient of the nLL.
If h(b) is some other function of the parameters b (e.g., the metabolite AUCs), then the variance of h(b̂) can be approximated by
To compute ∇h(b̂) where h(b) = ∫_{0}^{48}A_{i}(b,t)dt is the AUC for metabolite i, thenwhere Z_{k}^{j} are the solutions to the adjoint equations (Eq. J) as described previously.
We also found it potentially useful to compute a joint confidence (probability) that each putative second dose would meet prespecified AUC targets. Given a specified second cyclophosphamide dose and the resultant predicted AUC pair (x̄,ȳ), where x = AUC HCY and y = AUC CEPM, the joint confidence that x ≥ 50 μmol/L·h and y ≤ 350 μmol/L·h iswhere the probability density function is.
The joint confidence is computed in the BaRD program using functions from the MVTNORM package (18) in R.
Full Bayesian approach. An alternative to the MAP method might be to employ a full Bayesian, Markov chain Monte Carlo approach, possibly using WinBugs (23), as was done for two population pharmacokinetic studies in ref. 24. For our individual pharmacokinetic problem, where a range of second doses is investigated based on metabolite AUCs, the solution would result in posterior distributions on the AUCs for each potential second dose. The mode and credible region of the AUCs would be used in the clinic in the same manner as the AUC means and joint confidence regions (Eq. K) to select the second cyclophosphamide dose.
An advantage of the full Bayesian method would be the direct computation of the AUC posterior distributions. Disadvantages in our clinical setting would likely include longer and more complicated computer runs and the potential necessity of expert analysis to interpret the results. Testing would be necessary to determine feasibility of this approach for our fourcompartment (four ordinary differential equation) model and clinical setting as well as to determine if AUC targeting results or clinical efficacy is improved. Implementation and testing of a full Bayesian approach would comprise an interesting topic for future research.
Results
Figures 2 and 3 compare the BaRDpredicted HCY and CEPM concentrations, respectively, over 48 hours to measured plasma concentrations for the 20 phase 1 trial patients (15). The model was conditioned only on the 0 to 16hour concentration data (circles). Further concentration data over the 16 to 48hour period was plotted (triangles) for comparison with predicted concentrations (where the solid line represents the fit of the model to metabolite concentration data from the first dose and the dashed line represents the predicted concentration profile from the additional dose). Agreement between measured and predicted concentration values seems acceptable in all cases.
Table 2 provides the individual parameter estimates for the cyclophosphamide model (Eq. A) in comparison with those obtained using the POSTHOC option available in NONMEM. The model was informed only by the population pharmacokinetic parameter estimates from the earlier trial (14) and the individual concentration measurements through 16 hours after the first cyclophosphamide dose. Individual parameter values compared well with values (mean ± SD) for CL_{NON}, CL_{IND}, K_{ENZ}, K_{HCY}, K_{CEPM}, and K_{34} estimated by NONMEM being 1.42 ± 0.91 10^{−2} L/h/kg, 4.38 ± 2.87 10^{−2} L/h/kg, 1.50 ± 0.46 10^{−2} L/h, 1.38 ± 0.22 10^{2} L/h, 1.13 ± 0.24 L/h, and 2.00 ± 0.41 L/h and those estimated by BaRD being 1.46 ± 0.94 10^{−2} L/h/kg, 4.21 ± 2.50 10^{−2} L/h/kg, 1.50 ± 0.43 10^{−2} L/h, 1.40 ± 0.20 10^{2} L/h, 1.13 ± 0.26 L/h, and 2.08 ± 0.46 L/h, respectively. The linear correlation (as measured by r^{2}) between NONMEM and BaRDprovided values was 0.999, 0.998, 0.988, 0.998, 0.974, and 0.917 for CL_{NON}, CL_{IND}, K_{ENZ}, K_{HCY}, K_{CEPM}, and K_{34}, respectively.
The resulting nLL values are provided as a way of comparing the parameter optimizations between software implementations (a smaller value is better in this case). The nLL was evaluated for both NONMEM (20.23 ± 3.40, unitless) and BaRD (20.14 ± 3.16, unitless) parameterizations using the likelihood evaluation code developed in R. In principle, even slightly different parameterizations would provide different values on the nLL surface. Note that our method provided a better likelihood value (i.e., lower nLL or higher likelihood) than that provided by the parameterization chosen by NONMEM in all 20 cases, although only slightly better in most cases. However, this is likely a function of stopping tolerances only and not necessarily a measure of method efficacy. In fact, nLL values correlated extremely well (despite the relatively narrow range), with r^{2} = 0.999 between the two implementations.
Table 3 provides predicted and observed HCY and CEPM AUCs for the 20 study patients and compares predictions based on NONMEM and BaRD. An exhaustive comparison between databased and modelpredicted AUC values was conducted (15), which indicated that the population pharmacokinetic model was performing very well. Agreement between the nonlinear mixedeffectsbased population pharmacokinetic methods was very good. For CEPM, the correlation between the observed AUC and the NONMEMpredicted AUC had r^{2} = 0.756, whereas the correlation between the observed AUC and BaRD had r^{2} = 0.724. Correlation between NONMEM and BaRD AUCs was extremely good, with r^{2} = 0.993. For HCY, these correlations had r^{2} = 0.839, 0.838, and 0.974, respectively.
Tables 4 and 5 provide examples of how the BaRD tool may be used in choosing an individualized second cyclophosphamide dose. The BaRD code provides a range of second cyclophosphamide doses with resultant predicted AUCs of HCY and CEPM as well as joint confidences that the AUCs are contained in prespecified target ranges. We wish to choose the second cyclophosphamide dose so that the 0 to 48hour AUC is >50 μmol/L·h for HCY and between 300 and 350 μmol/L·h for CEPM. Joint confidences can provide additional guidance to a trained professional in selecting an appropriate second cyclophosphamide dose based on these particular targets.
Discussion
Over the past 10 years, increasing attention has been given to the inability of dosing according to weight or body surface area to account for individual variation in pharmacokinetics and clinical response to anticancer agents (25–29). Several other dosing strategies have been proposed, with dose adjustment based on individual pharmacokinetic data being one alternative in those instances where pharmacokinetic/pharmacodynamic relationships are definable (e.g., busulfan and methotrexate; refs. 30, 31). Noncompartmental methods using individual patient data have been the traditional pharmacokinetic method to estimate AUC, and thus dose, based on the relationship AUC = dose/clearance. However, there is gaining acceptance of the use of population pharmacokinetic modeling for anticancer agents in the clinical setting, for example, predicting leukopenia after epirubicin/doxetaxel (32) and realtime dose adjustment of topotecan (33). Cyclophosphamide is of particular interest because of its complex metabolic fate and the definable relationships between its metabolite concentration and occurrence of sinusoidal liver toxicity (2). A dose adjustment scheme for three anticancer agents (cyclophosphamide, thiotepa, and carboplatin) using a population pharmacokinetic model in women with breast cancer was recently published (13). Cyclophosphamide was administered as a 1hour infusion daily for 4 days; on days 1 and 2, standard dose cyclophosphamide was administered; on days 3 and 4, individualized cyclophosphamide doses were administered. The dose individualization was done based on pharmacokinetic sampling conducted only on day 1, thus allowing a workup of over 24 hours for drug metabolite assay and modeling of cyclophosphamide and HCY pharmacokinetics using NONMEM. For our hematopoietic cell transplant patients receiving cyclophosphamide/total body irradiation, a maximum of 8 hours is available between the time of the last pharmacokinetic sample and the time for administration of the second cyclophosphamide dose. The population pharmacokinetic model for cyclophosphamide and its relevant metabolites (i.e., HCY and CEPM) is complex, incorporating both enzyme autoinduction and nonlinear parameters (Fig. 1; Eq. A). This makes the implementation of a population pharmacokinetic model for cyclophosphamide particularly challenging for realtime dose individualization. The BaRD code we have written within R for cyclophosphamide dose adjustment can be modified for other medications of interest, for example, by using the published population pharmacokinetic model for topotecan (33) or epirubicin/docetaxel (32).
Population pharmacokinetic analysis is gaining increased utilization to estimate betweensubject variation in pharmacokinetic parameters. The output of a population pharmacokinetic model comprises both central tendencies (the fixed effects) and variability measures (the variances of the random effects) for the relevant parameters that shape a patient's pharmacokinetics and/or pharmacodynamics. As such, these parameters provide a probability distribution of pharmacokinetic results, and one application of knowing these probabilistic distributions is to assist in the estimation of an individual's parameters based on very sparse data.
In this work, we have used a previously constructed population pharmacokinetic model (14) to develop individualized dosing strategies for cyclophosphamide metabolites that are easily applicable in realtime in a clinical trial setting. The procedure is implemented in a userfriendly, opensource software and can in principle be modified to characterize dosing of other medications. Although pharmacokinetic dose tailoring has been practiced for many years, this is the first time to our knowledge that this procedure is applied to a medication with highly nonlinear, multimetabolite pharmacokinetics where parent drug concentrations are not available. It is easy to envision how this framework can be extended to other drugs; the modular nature of our software ensures that this can be done simply by updating or changing drugspecific modules, such as those that implement the pharmacokinetic model equations. The nonlinear regression routine would be similar for any drug.
As a validation of application of BaRD for realtime adjustment of the second cyclophosphamide dose, we compared the BaRD output with both an alternative implementation of the same model using the NONMEM population pharmacokinetic software and a regression equation of the AUC 0 to 16 hours to AUC 0 to 48 hours. Comparison in both cases was very favorable: the former shows that BaRD affords nearly identical predictions as NONMEM, whereas the latter result supports the superiority of population pharmacokinetics over the empirical regression approach.
Acknowledgments
We thank Dr. Bradley M. Bell for assistance in the development of the approaches used in this work and Scott Cole and Ami Batchelder for assistance in the clinical trial. We also wish to thank the anonymous reviewers for their helpful and constructive comments.
Footnotes

Grant support: Department of Pharmacy, University of Washington (D.H. Salinger, J.S. McCune, and D.D. Shen), NIH/National Cancer Institute grant CA18029 (J.S. McCune, J.T. Slattery, and G.B. McDonald), National Institute of Biomedical Imaging and Bioengineering grant EB001975 (P. Vicini), and Eli Lilly & Co. Foundation predoctoral fellowship (A.G. Ren).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.
 Accepted June 8, 2006.
 Received September 23, 2005.
 Revision received May 12, 2006.