We enrolled 60 patients with hematologic malignancies subjected to stem cell therapy. Patients were at least 15 years old, displayed adequate cardiac, hepatic and renal functions, and Karnofsky performance scores [13] of 70 or higher. Patient characteristics are presented in Table 1.

Subjects received 3.2 mg/kg/day iBu following either of two treatment regimens, specifically, 4 times daily iBu×4 days (BU4 arm) or once-daily iBu×4 days (BU1 arm) as conditioning therapy for stem cell transplantation. Patients were randomly assigned to the two treatment groups of 30 each. Block randomization method was employed, including stratification according to the conditioning regimen (busulfan-cyclophosphamide [BuCy] versus busulfan-fludarabine-antithymocyte globulin [BuFluATG] versus busulfan only [Bu]). Randomization was carried out centrally by the pharmacist using computer-generated random number tables. The treatment allocation was concealed from the investigators until 1 week before the administration of the study drug.

Hematopoietic cell grafts were infused on day 0 (for bone marrow) or days 0 and 1 (for granulocytecolony stimulating factor [G-CSF] mobilized peripheral mononuclear cells) without T cell depletion. For the BuCy regimen, intravenous busulfan (3.2 mg/kg/day) was administered on days -7 to -4 and cyclophosphamide (60 mg/kg/day) on days -3 and -2. The time between the last dose of busulfan and the first dose of cyclophosphamide was 14 hours in the BU4 arm and 27 hours in the BU1 arm. For the BuFluATG regimen, we administered intravenous busulfan (3.2 mg/kg/day) for 2 days (days -7 and -6), fludarabine (30 mg/kg) for 6 days (days -7 to -2), and antithymocyte globulin on days -4 to -2 with a matched sibling donor and an unrelated donor or -4 to -1 with haplo-identical familial donor. For the Bu regimen, intravenous busulfan (3.2 mg/kg/day) was administered on days -6 to -3. All patients received an intravenous loading dose of phenytoin (15 mg/kg) the day before the first busulfan administration, and oral dosing was continued to maintain therapeutic levels (10 to 20 mg/L) until the day after the last dose of busulfan.

Patients in the BU4 arm received iBu (0.8 mg/kg) every 6 h in 2 h infusions, while those in the BU1 arm received 3.2 mg/kg iBu every 24 h in 3 h infusions. Busulfan was diluted in normal saline to 0.5 mg/ml, and introduced using an infusion pump through a central venous catheter. Doses of busulfan were calculated using selected body weight (SBW) which was: (1) actual body weight (ABW) if less than or equal to ideal body weight (IBW), (2) IBW if ABW was more than IBW but within 120% of IBW or (3) Adjusted ideal body weight (AIBW)="IBW+0.40·(ABW-IBW)" if ABW exceeded IBW by more than 120% [14]. IBW was estimated using the following equation, measuring height in inches and weight in kilograms: (1) IBW (men)=50+2.3·(height-60) or (2) IBW (women)=45+2.3·(height-60) [15]. Venous blood samples (5 ml) were obtained from all patients at 5 time-points after the first dose of busulfan therapy using the limited sampling strategy adopted from a previous study [16]. Venous blood was drawn at 2.5, 3, 4, 5, and 6 h after the start of the infusion from patients in the BU4 arm, and at 3.5, 5, 6, 7 and 22 h after infusion from those in the BU1 arm. Samples were obtained via a peripheral venous catheter. One ml of blood was discarded before the collection of blood specimen and sterile 0.9% saline (1 ml) was injected into the catheter after each blood sampling procedure. The blood sample was introduced into pre-chilled heparin tubes, and within 30 minutes, plasma was separated by centrifugation at 1,286 G over 10 min at 4℃. Plasma samples were stored at -40℃ until analysis. The blood samples were analyzed for glutathione S-transferase (GST) genetic variants GSTM1 (null allele) and GSTT1 (null allele) as described by Arand et al. [17]. The protocol was approved by the Institutional Review Board (IRB) of the Asan Medical Center (IRB registration number 2004-068). All subjects gave their written informed consent before participating in the study.

The plasma concentration of busulfan was measured by validated liquid chromatography with tandem mass spectrometry (LC-MS/MS) in a method similar to that described by dos Reis et al. [18] performed on an API 3000® triple quadruple mass spectrometer equipped with an electrospray ion source (MDS SCIEX, South San Francisco, CA, USA). An aliquot of the sample (20 µl) was delivered to the electrospray ion source using HPLC (Agilent 1100 series, Agilent Technologies Inc., Santa Clara, CA, USA) with a C18 Capcell Pak® MG column (2.0×50 mm, 3.0 µm particle size). The mobile phase comprised acetonitrile, tetrahydrofurane and distilled water (65:5:30). For validation procedures, plasma calibration curves, each comprising six levels of busulfan (30~6,000 ng/ml) and a fixed concentration of internal standard (metronidazole 500 ng/ml), were prepared and assayed. To assess the intra- and inter-day precision and accuracy of the method, five replicates of the plasma standards at three concentrations (40, 400 and 4,000 ng/ml) were analyzed. Calibration curves were linear throughout the concentration range of the study, with correlation coefficients greater than 0.998 for all cases. Based on a signal-to-noise level of 10, the quantification limit for busulfan was calculated as 30 ng/ml. The intra-day CV was ≤10.14%, and intra-day accuracy ranged from 94.10% to 107.80%, while the inter-day CV was ≤6.29%, with accuracy ranging from 95.59% to 101.44%. The bio-analysis was done at the Pharmacokinetics laboratory, Clinical Reseearch Center, Asan Medical Center (Seoul, Korea).

In total, 295 measurements from 60 patients were analyzed by mixed-effect modeling using NONMEM® (Version VI, GloboMax LLC, Ellicott City, MD, USA). The pharmacokinetic parameters were estimated with NONMEM subroutines ADVAN1 TRANS2, using the First Order Conditional Estimation with interaction (FOCEi) method. The parameters for a specific subject were described using the following equation:

where P_TV is the typical value of the parameter, and ηⁱ is a normally distributed variable with zero mean.

The residual error model was characterized with the proportional error model described using the following equation:

where ε represents a zero-mean normally distributed variable.

Various structural pharmacokinetic and error models were assessed, guided by the graphical assessment of optimum fit properties and statistical significance criteria. A likelihood ratio test was applied to discriminate between the reduced and full models at a significance level of p≤0.05, equivalent to a change of 3.84 in the objective function value. Standard diagnostic plots, including the observed values of the dependent variable (DV) versus individual predicted values (IPRE) and IPRE versus individual weighted residuals, were used for the diagnosis of optimum fit capabilities. Standard errors of parameter estimates of the pharmacokinetic model were employed as a diagnostic.

Potential covariates affecting the clearance (CL) and volume of distribution (V_d) were explored. A regression model for each structural model parameter was constructed in three steps using the original dataset. Individual covariates were initially screened. The full model was defined as incorporating all significant covariates. The final model was elaborated by backward elimination from the full model. For each analysis, the improvement in fit obtained upon the addition of a covariate selected from step 2 to the base model was assessed by changes in the NONMEM® objective function value (OFV). To discriminate between the reduced and full models, a significance level of p≤0.05, equivalent to a change of 3.84 in the objective function value, was applied.

In the first step, the following covariates for iBu pharmacokinetics were screened: sex, age, actual body weight (ABW), ideal body weight (IBW), adjusted ideal body weight (AIBW), selected body weight (SBW), height, body mass index (BMI), and body surface area (BSA) calculated from five different equations, serum creatinine, creatinine clearance, aspartate transaminase (AST), alanine transaminase (ALT), alkaline phosphatase (ALP), albumin, total bilirubin, and genetic polymorphisms in GSTM1 and GSTT1. Serum creatinine was measeured on the day of the study and creatinine clearance was estimated from serum creatinine using the Cockroft and Gault equation [19]. The following five different BSA equations were used to find out which would better correlate with pharmacokinetics of iBu:

Mosteller [20] formula, BSA (m²)=((Height (cm)·Weight (kg))/3,600)^½

Du Bois and Du Bois [21] formula, BSA (m²)=0.20247·Height (m)^0.725·Weight (kg)^0.425

Haycock et al. [22] formula, BSA (m²)=0.024265·Height (cm)^0.3964·Weight (kg)^0.5378

Gehan and George [23] formula, BSA (m²)=0.0235·Height (cm)^0.42246·Weight (kg)^0.51456

Boyd [24] formula, BSA (m²)=0.0003207·Height (cm)^0.3·Weight (g)^{(0.7285-0.0188·log}₁₀^(weight(g)).

The distribution of empirical Bayesian parameter estimates was explored, and the relationships between covariates and individual pharmacokinetic parameter estimates evaluated. Data were subjected to a stepwise (single term addition/deletion) procedure using the generalized additive model (GAM) [25] in which Akaike's information criterion (AIC) [26] was applied for model selection.

After the final covariate model was built, population shrinkage of interindividual random variability (η) was calculated as follows:

Where ω is the estimate of inter-individual variability and SD() represents the standard deviation of the empirical Bayesian estimate of η for each individual. Shrinkage is smaller when data are more informative.

Random permutation tests [27] of 2,000 samples were performed to examine the statistical significance of the covariates. The tested covariate was considered significant if the OFV from the original data set was below the 2.5^th percentile of OFVs from randomized datasets.

Two thousand datasets were simulated from the final pharmacokinetic model using NONMEM VI, and the 95% prediction interval compared visually and numerically with actual plasma busulfan concentration data. The simulated datasets were re-fitted by the final model, and posterior distribution of the parameter estimates compared with the original final parameter estimates. The observed and simulated concentrations were compared visually using mirror plots. Bootstrapping, a resampling technique with replacement, was performed for the bias and stability of parameter estimates. In total, 2,000 bootstrap runs were performed, and 95% confidence intervals of the parameter estimates obtained as 2.5^th and 97.5^th percentiles from the resultant parameter distributions. Cook and weisberg [28] score and covariance ratio [29] were estimated to detect influential subjects.

The modeling process was facilitated by Xpose® [26] (version 4.0) run on the R statistical software package (version 2.6.0, The R Foundation for Statistical Computing, Vienna, Austria, URL http://www.R-project.org) and Asan Software Tool for NONMEM, an interface for NONMEM® based on text editor and R.

An equation for dose calculation on a once-daily basis was derived from the population estimate of iBu clearance and target AUC from time 0 to infinity (AUC_0-inf). The median of AUC_0-inf, evaluated by dividing the drug dose by individual busulfan clearance estimates from patients in the BU1 arm, was applied as the therapeutic target.

Two datasets of 2,000 hypothetical patients were created using the Monte Carlo simulation, one using a 3.2·SBW mg q24h scheme, and the other with a newly proposed scheme based on the final population PK model. Busulfan concentration profiles were simulated assuming a single dose administration with the same plasma concentration measurement schedule as the BU1 arm. AUC_0-inf values were generated from each simulated patient using clearance estimates. The probabilities of attaining an AUC_0-inf range of target AUC_0-inf±10% and ±20% were calculated.

Plasma log-concentration versus time curves after iBu infusion disclosed a linear relationship (Fig. 1). A one-compartment model with exponential inter-individual variability and a proportional error model were optimal to describe the time-concentration curve. In terms of GAM analysis, ABW, sex, GSTM1, GSTT1, and total bilirubin were selected for CL, and ABW, sex, and ALP for V_d. The population model with covariates was built using the NONMEM® program on the basis of GAM analysis data. ABW was a predictor of both CL and V_d with an objective function value difference of more than 3.84 (p=0.05) between each model in which ABW was introduced alone and the basic model of each pharmacokinetic parameter without ABW. The best fit was obtained with CL and V_d modeled as power functions of ABW, and with V_d modeled with a sexual difference term. Initial estimation (standard error) of power terms for CL and V_d were 0.36 (2.56) and 0.48 (0.154). Simplified power model with the power terms fixed at 0.5 was fitted to the data and estimation of power terms as unknown parameters was found not to improve the model compared to fixing the power term. The final covariate model was as follows:

CL=(θ₁·ABW^0.5)·(exp(η₁))

V_d=(θ₂·ABW^0.5·(1+SEX·θ₃))·(exp(η₂))

with CL presented in L/hr, ABW in kilograms, V_d in L and SEX coded as female=0 and male=1.

Parameter estimates of the final covariate model are presented in Table 2.

Plots of observed versus predicted concentrations for the final covariate model are shown in Fig. 2. The shrinkage of η₁ and η₂ were 1% and 9%. A random permutation test showed that the 2.5^th percentiles of OFVs from randomized datasets for ABW and sex exceeded the OFVs from the original dataset, confirming that ABW and sex are significant covariates. The 5^th and 95^th percentiles (prediction intervals) of simulated dose-normalized concentrations were calculated and plotted against the observed concentrations (Fig. 3). The model predicted the simulated concentrations fairly accurately, but the observed concentrations were slightly more variable. A numerical predictive check was performed to evaluate the stability of the final model. For each observation, 2,000 predictions were generated, and the corresponding 50%, 90%, and 95% prediction intervals defined. Ideally, 25%, 5%, and 2.5% of the observations would be above and below the 50%, 90%, and 95% prediction intervals, respectively. Out of the 295 observed busulfan concentrations, 2.72%, 5.70%, and 25.85% were above, and 1.70%, 5.78%, and 22.79% were below prediction intervals respectively. Two thousand datasets were simulated using the estimated parameters and the final covariate model, and each dataset was re-fitted using the final covariate model. In a posterior predictive check, posterior distributions of parameters derived from the simulated datasets were evaluated. Distribution patterns of pharmacokinetic parameters from a single run with the original dataset were comparable to those of parameters derived by fitting the final covariate model to simulated datasets. Simulated datasets were compared visually with the original dataset. Dispersion patterns around the lowess lines were similar between the original observations and simulated concentrations. The final population pharmacokinetic model obtained from the previous step was fitted repeatedly to 2,000 bootstrapped samples. In all runs, the minimization and covariance steps were successful. The parameter estimates of the final model using the original data and the mean parameter estimates from the 2,000 bootstrap replicates are presented in Table 3. The mean bootstrap parameter estimates were within 1% of those obtained with the original dataset. Case deletion diagnostics revealed that no individual had a high Cook score together with a low covariance ratio (<0.5), indicative of an influential subject.

The known therapeutic window for 4 times daily oral busulfan, AUC_0-6 target of 900~1,500 mol/l·min [4-7], is not applicable for once-daily administration of iBu. The median AUC_0-inf after the first dose of the 4 times daily regimen of iBu was 1,481 mol/l·min, comparable to target AUC_0-6 of 1,200 mol/l·min. Since the pharmacokinetic profiles and post-transplant complications were similar for the BU1 and BU4 arms in this study [12], the median AUC_0-inf of patients in the BU1 arm, 6,378 mol/l·min (=26.18 mg/l·hr), was used as the therapeutic target. Based on the CL estimate of 0.947 l/hr ·ABW^0.5, the appropriate dose to achieve target AUC_0-inf was calculated as CL·AUC_0-inf=24.79·ABW^0.5 mg q24h. ABW is measured in kilograms.

Two datasets of 2,000 hypothetical patients were created using a Monte Carlo simulation, one with a 3.2·SBW mg q24h scheme, and the other with a proposed scheme of 24.79·ABW^0.5 mg q24h. The probabilities of attaining a AUC_0-inf range of target AUC_0-inf±10% (23.56~28.80 mg/l·hr) and ±20% (20.94~31.42 mg/l·hr) were 69.2% and 96.9% with the conventional scheme, and 89.1% and 99.5% with the new scheme, respectively.