Population pharmacokinetic model for once-daily intravenous busulfan in pediatric subjects describing time-associated clearance.

This study aimed to characterize the population pharmacokinetics (PK) of busulfan focusing on how busulfan clearance (CL) changes over time during once-daily administration and assess different methods for measuring busulfan exposure and the ability to achieve target cumulative exposure under different dosing adjustment scenarios in pediatric stem cell transplantation recipients. Daily serial blood sampling was performed and concentration-time data were analyzed using a nonlinear mixed-effects approach. The developed PK model was used to assess achievement of target exposure under six dose-adjustment scenarios based on simulations performed in RStudio (RxODE package)®. A total of 2491 busulfan plasma concentration-time measurements were collected from 95 patients characterizing 379 dosing days. A two-compartment model with time-associated CL best described the data with a typical CL of 14.5 L/h for an adult male with 62 kg normal fat mass (NFM; equivalent to 70 kg total body weight), typical volume of distribution central compartment (V1) of 40.6 L/59 kg NFM (equivalent to 70 kg total body weight), and typical volume of distribution peripheral compartment of 3.57 L/62 kg NFM. Model interindividual variability in CL and V1 was 14.7% and 34.9%, respectively, and interoccasional variability in CL was 6.6%. Patient size described by NFM, a maturation component, and time since start of treatment significantly influenced CL. Simulations demonstrated that using model-based exposure estimates with each dose, and either a proportional dose-adjustment calculation or model-based calculated individual CL estimates to support dose adjustments, increased proportion of subjects attaining cumulative exposure within 5% of target compared with using noncompartmental analysis (100% vs. 0%). A time-associated reduction in CL during once-daily busulfan treatment was described.

WHAT IS THE CURRENT KNOWLEDGE ON THE TOPIC?

Previous publications have reported an 8.1 to 20% reduction in busulfan clearance over a typical four‐day treatment course. , , , , , , Product information leaflets recommend using a proportional equation for dose adjustment which does not account for changes in clearance over time.

WHAT QUESTION DID THIS STUDY ADDRESS?

The population pharmacokinetic model developed in this study characterises the time‐associated reduction in busulfan clearance based on data collected across all dosing days. Simulations were performed to access different methods to estimate exposure and perform dose adjustments.

WHAT DOES THIS STUDY ADD TO OUR KNOWLEDGE?

Simulations based on non‐compartmental analysis (NCA) to determine busulfan exposure resulted in all patients achieving cumulative exposure above the target, even when samples were taken following each administered dose. Using a model‐based method to estimate exposure combined with daily sampling resulted in 100% of subjects being within ±5% of target when either proportional dose adjustment or model calculated individual CL was used to calculate next doses.

HOW MIGHT THIS CHANGE DRUG DISCOVERY, DEVELOPMENT, AND/OR THERAPEUTICS?

Busulfan exposure is better calculated using model‐based methods (preferably involving a population pharmacokinetic model characterising change in clearance over the treatment course) rather than an NCA approach. Monitoring and dose adjustment after each dose increases the proportion of patients achieving target exposure. Bayesian forecasting software would allow for model‐based individual CL estimates to be utilised for dosing in everyday clinical practice.

INTRODUCTION

Busulfan is a bifunctional, alkylating agent commonly used in hematopoietic stem cell transplant (HSCT) conditioning regimens. Several studies have reported that intravenous (i.v.) busulfan clearance (CL) decreases during the course of treatment, , , , , , , with a reduction in CL ranging from 8.1% to 20%. It has been hypothesized that this may be due to busulfan metabolites competing with free busulfan for binding to glutathione‐s‐transferase (GST) enzymes and glutathione, resulting in reduced metabolism and glutathione depletion. , , , One publication explored a semimechanistic model based on glutathione depletion to predict reduction in CL of busulfan. Investigators observed that the effect of glutathione depletion on CL of busulfan was proportional and increased in patients aged >40 years.

Toxicity with busulfan therapy is well described, with monitoring to optimize busulfan exposure and improve treatment outcomes considered standard practice. , The American Society for Blood and Marrow Transplantation practice guidelines committee attempted to develop an evidence‐based review about personalizing busulfan‐based conditioning. However, published literature was too heterogenous and lacked adequately powered controlled studies to provide a complete review, with guidance on only some aspects of personalizing busulfan dosing in children. Exposure monitoring should be performed due to an association between busulfan exposure and clinical outcomes, , and validated pharmacokinetics (PK) modeling tools should be used. PK‐guided dosing of i.v. busulfan was shown to be superior to fixed dosing based on patient body surface area in adults with acute myeloid leukemia and myelodysplastic syndrome with reduced relapse (38% vs. 56%), transplant‐related mortality (24% vs. 39%), and an overall hazard ratio of 0.64 (95% confidence interval, 0.35–0.94) observed.

Busulfan exposure is estimated using various indices across the literature and in clinical practice, including the area under the concentration‐time curve (AUC) over the dosing interval, the AUC over a 24‐h period (AUC0–24), and cumulative AUC following all doses (AUCcum). AUCcum is increasingly being used for both the reporting of outcomes and for individual dose targeting during exposure monitoring. , , , , , , , Bartelink et al. recommends an AUCcum between 78 and 101 mg · h/L for optimal survival. Of note, AUCcum in this publication was estimated using a model‐based approach rather than noncompartmental analysis (NCA). There is an expected bias toward underestimation with NCA compared with model‐based methods. ,

Busulfan dose adjustment across the typical 4‐day treatment course can vary across treatment centers. , Many centers adjust dosing based on the assumption that busulfan CL remains constant over time using product information recommendations , despite contrary evidence of reduction in CL during the treatment course. , , , , , , Busulfan concentration measurement is commonly performed after the first dose only, with subsequent dosage adjustment made on Days 2 or 3. Application of a population model describing changes in CL of busulfan over time in a Bayesian forecasting program could assist with making dosage adjustments. , ,

A recent review of published non‐linear mixed effects (NLME) models for i.v. busulfan in children identified 21 studies ; however, few focused on once‐daily therapy or included samples from all 4 days of treatment to enable characterization of temporal changes in busulfan CL. Most common covariates shown to explain PK variability were patient total body weight (TBW), age, Glutathione‐S‐Transferase‐ A1 genotype, and busulfan dosing day/time since start of therapy.

This study aimed to (1) characterize the population PK of busulfan focusing on how busulfan CL changes over a typical 4‐day treatment course by using a unique data set where sampling was performed following each administered dose and (2) assess different methods for estimating busulfan exposure and their ability to achieve target cumulative exposure under different dosing adjustment scenarios in pediatric subjects receiving once‐daily i.v. busulfan for HSCT.

METHODS

Subjects and data collection

Data were collected both prospectively and retrospectively across three hospitals in Australia and one in New Zealand from 2016 to 2021. Ethics approval for the study was obtained from the human research ethics committees from all centers involved (HREC/16/QRCH/388, HREC:2017000073/HREC/16/RCH/388, HREC:RGS0000000497, HREC:18/CEN/10). Subjects were eligible for enrollment if they were HSCT recipients aged <18 years receiving once‐daily i.v. busulfan and had adequate central venous access for blood sampling. All subjects received an initial once‐daily i.v. infusion of busulfan between 3.2 and 4.8 mg/kg over 3 h, with dose calculated based on patient total weight as per Australian Busulfex® product information except for two obese patients where an adjusted dosing weight (adjusted body weight 25) was used. Subsequent doses across 4 days of treatment were modified to target an AUCcum that was individualized for each patient based on clinical and disease factors. , , , Blood samples for busulfan plasma concentration measurement were collected following each dose, from an alternative lumen to that used for administration, at 0 h (predose), 3 h (immediately following first flush), and 3.25, 4, 5, 6, and 8 h after the start of the infusion. Concentrations equal to or below the lower limit of quantification (LLOQ) were removed from the final data set used to build the NLME model. Demographic data including date of birth, age, sex, body surface area, weight, height and clinical data including diagnosis, liver function tests, measured glomerular filtration rate, serum albumin, prior chemotherapy, HSCT conditioning regimen, immunotherapy, supportive care medications, and coadministration of interacting medications including azole antifungal agents, metronidazole, acetaminophen, and seizure prophylaxis were collected. Busulfan doses administered, timing of administration, duration of infusion (DI), and individual subjects' AUCcum targets were all determined by the treating physician.

Busulfan drug assay

Samples collected from July 2014 onward for busulfan concentration measurement were analyzed at Pathology Queensland using a validated reverse phase ultra‐high performance liquid chromatography method with tandem mass spectrometry detection (UPLC‐MS/MS). Following separation of plasma via centrifugation, the plasma samples were precipitated with methanol containing deuterated busulfan. Samples were then centrifuged, and supernatant was transferred to an autosampler vial and injected onto the UPLC‐MS/MS system. Busulfan and internal standard were separated chromatographically using a C18 reverse phase column. Multiple reaction monitoring was then carried out for each individual analyte. The gradient was then returned to initial conditions in preparation for the next sample. Intra‐ and interrun imprecision was less than 6%. Since July 2014, the assay has been shown to be linear from 0.01 to 20 mg/L with an LLOQ of 0.01 mg/L. For samples analyzed before July 2014, a similar assay was used with an LLOQ of 0.1 mg/L using precolumn derivatization with reverse‐phase high‐performance liquid chromatography. Subjects from Queensland had their samples tested immediately, whereas subjects from interstate or overseas centers had their separated plasma immediately frozen and stored and then transported on dry ice to the measurement laboratory no later than 6 months from freezing. Interstate and overseas subjects had separate samples tested in real time at their local laboratory to guide dose adjustments.

Model development and evaluation

NLME modeling PK analysis was performed using NONMEM® (Version 7.4.3, ICON) in conjunction with Perl Speaks NONMEM (Version 5.0) and Pirana (Version 2.9.9). R Studio® (Version 4.0.5; http://www.r‐project.org) and Xpose® (http://xpose.ourceforge.net) were used for data exploration and visualization. Initially, one‐, two‐, and three‐compartment models were trialed, with goodness‐of‐fit (GOF) diagnostic plots, evaluation of population fixed‐ and random‐effects parameters, prediction‐corrected visual predictive checks (pcVPC), and changes to the objective function value (OFV) used for model selection. Medians of the bootstrap estimate distributions and 95% confidence intervals from a nonparametric bootstrap analysis (N = 1000) were compared with final model estimates. Model parameters were estimated using first‐order conditional estimation with interaction with predictions for population pharmacokinetics library models. A convergence criterion of three significant digits was used to identify successful minimization.

Information on duration of busulfan infusion (DI) was initially applied as recorded in patient medical records and was later included as a parameter to estimate. Interindividual variability (IIV) was explored on all parameters using an exponential model assuming a log‐normal distribution. Additive, proportional, and combined additive and proportional residual error models were investigated.

Changes in busulfan CL over the course of treatment were explored during structural model building by including a covariate effect to examine the change in CL on each day in comparison to Day 1. A Michaelis–Menten concentration‐dependent model and an empirical continuous time‐associated model to describe the change in CL across dosing days were also tested. In the latter, the decrease in CL was described as shown in Equation (1): where CL is CL in the ith individual at a particular time t, ΔCLmax, is the maximal possible change in CL in the ith individual, t is time after the start of the infusion of the first dose in the ith individual in h, T 50, is the time at which 50% of ΔCLmax is attained in h in the ith individual, and γ describes the shape of the relationship. IIV on both ΔCLmax, and T 50, were investigated using an exponential model assuming a log‐normal distribution and removed if not supported.

The influence of body size and composition on PK parameters (SIZEparam) were incorporated following theory‐based allometry as shown in Equation (2): where SIZEparam, is the fractional difference in allometrically scaled size compared with an adult with TBW of 70 kg, Size is the individual patient size in kilograms, and Sizestd was the standard weight in kg corresponding to a TBW of 70 kg in an adult. Size was expressed in terms of either TBW or normal fat mass (NFM) for the ith individual. Sizestd was set at either an NFM of 62 kg for CL, 59 kg for typical volume of distribution central compartment (V1) and typical volume of distribution peripheral compartment (V2), and 56.1 kg for intercompartmental clearance (Q), corresponding to an allometrically scaled TBW of 70 kg or set to 70 kg if TBW was the size descriptor. Size descriptors including TBW and NFM are commonly used in published population PK models for busulfan in pediatrics. , NFM has recently been used to identify differences between oncology versus non‐oncology cohorts. In keeping with prior published pediatric models, the allometric exponent (power) in Equation (2) was fixed a priori to 0.75 for CL and Q and 1 for V1 and V2.

A maturation component was included on CL a priori. This was tested using Equation (3) , , : where F mat describes the maturation process, PMA is postmenstrual age in weeks for the ith individual, , TM50 is the PMA in weeks where the CL of busulfan is considered to be 50% of the adult value, and the Hill coefficient represents the steepness of the function. Subjects in this study were assumed to have been born at term (i.e., 40 weeks PMA). The Hill coefficient was estimated and fixed to 2.3. TM50 was fixed to 45.6 weeks as supported by two studies involving 540 and 1610 patients, respectively. , Lastly, intraoccasional variability (IOV) was modeled using an additional random‐effects parameter and was tested one by one on CL, V1, Q, and V2.

Structural and stochastic models were selected based on GOF and a reduction in OFV by ≥3.84 points for one degree of freedom at p < 0.05 between nested models by model performance, including convergence and precision of parameter estimates. Thereafter, additional covariates were tested on CL in a stepwise process based on physiological plausibility and literature review and retained in the model if statistically significant based on a χ2 test (p < 0.05). Covariates tested included use of concurrent medications for seizure prophylaxis, antifungal agents, acetaminophen, and fludarabine and an oncology versus nononcology diagnosis. The influence of interacting medications was tested on both ΔCLmax,i and time at which 50% of ΔCLmax,i is attained (T 50,i) using Equation (1).

Dose‐adjustment simulations

Six different dose‐adjustment scenarios were compared based on random simulation of a virtual population of 1000 subjects whose age, height, sex, and TBW were sampled from distributions with the same mean and standard deviation as the study population using RStudio (RxODE® package). Sex was randomized with a similar distribution as the study population. Dose 1 (D1) was calculated from TBW according to the product information. A summary of the different scenarios and their underlying assumptions is provided in Table 1. Dosing scenarios were designed to reflect common situations described in the literature , , , and seen in clinical practice. , , Busulfan concentrations were simulated based on the model and sampling times at 0, 1, 2, and 4 h post infusion as currently recommended within the product information. , These were then used to obtain daily exposure estimates (AUC0–24) either determined using the developed model by numeric integration (signified by MOD in the scenario name) or using NCA (signified by NCA in the scenario name). Sampling from D1 only versus sampling following each of the four doses indicated as D1 or D1–4 in the scenario name was compared together with the method to calculate the next dose. Dosage adjustment was either based on a proportional relationship (as suggested in the product information [signified with PI in the scenario name]) or using individual CL (CL ) estimates from the model immediately prior to the next dose (signified with MOD in the scenario name). In all scenarios, an AUCcum of 90 mg · h/L was targeted, the middle point of the optimal target range recommended by Bartelink and colleagues. The proportion of subjects who obtained exposure within 5% of the 90 mg ·h/L target and the number and proportion of subjects who obtained exposure within the optimal target range of 78–101 mg · h were reported across the different scenarios. AUCcum reported for all scenarios was determined using the RxODE package in R using the final NLME model developed in this study. Using numeric integration AUCcum was determined from AUC0–24 for Doses 1–3 and then AUC0‐∞ for Dose 4. The sample R code for key components of the simulation can be found in the Supplementary Material.

TABLE 1

Description of six scenarios and their underlying assumptions applied in the simulations

Scenario name	Dose/s with concentration sampling	Dose/s at which change could occur	Method for calculation of AUC0–24	Method to calculate dose change	Method for calculating the change in dosage	Description and assumptions made
NCA_PI_D1	1	2	NCA	Product Information–proportional	Next dose calculated per product information equation based on NCA calculated AUC0–24 after Dose 1, administered Dose 1 dosage (Current dose), and the next dose target AUC0–24: Next dose targetAUC_0−24mg·h/L=90mg·h/L−AUC_0−24Dose1mg·h/LNumber of doses remaining, Next dosemg=Current dosemg·Next dose targetAUC_0−24mg·h/LActualAUC_0−24mg·h/L, where the number of doses remaining is 3	NCA used to calculate AUC0–24 for first dose only from samples taken at 0, 1, 2, and 4 h from completion of infusion Dose adjustment calculated using assumption that CL remains the same over the 4‐day course (proportional adjustment as per product information)
NCA_PI_D1‐4	1, 2, 3, 4	2, 3, 4	NCA	Product information–proportional	Next dose calculated per product information equation based on NCA calculated AUC0–24 after each nth dose, current dose administered, and the next dose target AUC0–24: Next dose targetAUC_0−24mg·h/L=90mg·h/L−AUC_cumAlldoses administeredmg·h/LNumber of doses remaining, Next dosemg=Current dosemg·Next dose targetAUC_0−24mg·h/LActualAUC_0−24from current dosemg·h/L, where the number of doses remaining is: 4–n	Assumption that busulfan concentration results available in time for adjustment of next dose NCA used to calculate AUC0–24 for all doses from samples taken at 0, 1, 2, and 4 h from completion of infusion Dose adjustment calculated using assumption that CL remains the same over the 4‐day course (proportional adjustment as per product information)
MOD_PI_D1	1	2	Model	Product information–proportional	Next dose calculated per product information equation using a model‐based calculated AUC0–24 after Dose 1, administered Dose 1 dosage (Current dose), and the next dose target AUC0–24: Next dose targetAUC_0−24mg·h/L=90mg·h/L−AUC_0−24Dose1mg·h/LNumber of doses remaining, Next dosemg=Current dosemg·Next dose targetAUC_0−24mg·h/LActualAUC_0−24mg·h/L, where the number of doses remaining is 3	Model incorporating time‐associated CL used to calculate AUC0–24 for first dose only from samples taken at 0, 1, 2, and 4 h from completion of infusion Dose adjustment calculated using assumption that CL remains the same over the 4‐day course (proportional adjustment as per product information)
MOD_PI_D1‐4	1, 2, 3, 4	2, 3, 4	Model a	Product information–proportional	Next dose calculated per product information equation using a model‐based calculated AUC0–24 after each nth dose, current dose administered, and the next dose target AUC0–24: Next dose targetAUC_0−24mg·h/L=90mg·h/L−AUC_cumAlldoses administeredmg·h/LNumber of doses remaining, Next dosemg=Current dosemg·Next dose targetAUC_0−24mg·h/LActualAUC_0−24from current dosemg·h/L, where the number of doses remaining is: 4–n	Assumption that busulfan concentration results available in time for adjustment of next dose Model incorporating time‐associated CL used to calculate AUC0–24 for all doses from samples taken at 0, 1, 2, and 4 h from completion of infusion Dose adjustment calculated using assumption that CL remains the same over the 4‐day course (proportional adjustment as per product information)
MOD_MOD_D1	1	2	Model	Model‐calculated CL i	Next dose calculated using model‐based calculated AUC0–24 after administered Dose 1, the model‐based calculated CL i after Dose 1, and the next dose target AUC0–24: Next dose targetAUC_0−24mg·h/L=90mg·h/L−AUC_0−24Dose1mg·h/LNumber of doses remaining, Next dosemg=CL_i,n·Next dose targetAUC_0−24, where the number of doses remaining = 3, CL i,n  = CL i calculated at the end of the nth dose time interval (immediately prior to next dose), n = 1	Model incorporating time‐associated CL used to calculate AUC0–24 for first dose only from samples taken at 0, 1, 2, 4 h from completion of infusion Dose adjustment calculated using model‐based individual CL i calculated immediately prior to Dose 2
MOD_MOD_D1‐4	1,2,3,4	2,3,4	Model	Model‐calculated CL i	Next dose calculated using model based calculated AUC0–24 after each nth dose, the model‐based calculated CL i after each nth dose, and the next dose target AUC0–24: Next dose targetAUC_0−24mg·h/L=90mg·h/L−AUC_cumAlldoses administeredmg·h/LNumber of doses remaining, Next Dosemg=CL_i,n·Next dose targetAUC_0−24, where CL i,n  = CL i calculated at the end of the nth dose time interval (immediately prior to next dose) and the number of doses remaining is: 4–n	Assumption that busulfan concentration results available in time for adjustment of next dose Model incorporating time‐associated CL used to calculate AUC0–24 for all doses from samples taken at 0, 1, 2, 4 h from completion of infusion Dose adjustment calculated using daily updated model‐based individual CL i estimate, incorporating time‐associated CL

Note: For all scenarios, initial dose was calculated according to product information based on simulated patient weight (five weight‐based strata: <9 kg–4 mg/kg/dose, 9 to <16 kg–4.8 mg/kg/dose, 16–23 kg–4.4 mg/kg/dose, >23–34 kg–3.8 mg/kg/dose, >34 kg–3.2 mg/kg/dose), and all individuals had a target AUCcum of 90 mg · h/L.

Abbreviations: actual AUC0–24, area under the concentration‐time curve over a 24‐h period calculated based on the current dose administered; AUCcum, cumulative area under the concentration‐time curve following all doses; AUCcum(all doses administered), cumulative area under the time versus concentration curve for all doses administered; CL, clearance; CL , individual model‐based clearance calculated at the end of the dosing interval for the current dose; CL , individual model‐based clearance calculated at the end of the dosing interval for the nth dose; current dose, current dose administered; n, dose number of current dose administered; NCA, noncompartmental analysis; next dose target AUC0–24, remaining exposure to the target AUCcum divided by the number of doses remaining.

Model‐based AUC0–24 and AUCcum were calculated by numeric integration using Rstudio® RxOde package from 0 to 24 h post dose.

,

,

,

,

,

,

,

Next dose calculated using model‐based calculated AUC0–24 after administered Dose 1, the model‐based calculated CL after Dose 1, and the next dose target AUC0–24:

,

,

where the number of doses remaining = 3, CL  = CL calculated at the end of the nth dose time interval (immediately prior to next dose), n = 1

Dose adjustment calculated using model‐based individual CL calculated immediately prior to Dose 2

Next dose calculated using model based calculated AUC0–24 after each nth dose, the model‐based calculated CL after each nth dose, and the next dose target AUC0–24:

,

,

where CL  = CL calculated at the end of the nth dose time interval (immediately prior to next dose) and the number of doses remaining is: 4–n

Dose adjustment calculated using daily updated model‐based individual CL estimate, incorporating time‐associated CL

RESULTS

Subjects and data

Data were collected from a total of 95 pediatric HSCT subjects (49 retrospectively, 46 prospectively) who had received i.v. busulfan administered once daily across 4 days. Samples below the LLOQ (5.3%) were all predose values and were not included in the data analysis. A total of 80 patients were from the Queensland Children's Hospital in Brisbane (Center 1), whereas the remaining 15 came from three interstate or overseas centers (Perth Children's Hospital, Sydney Children's Hospital, or Auckland Starship Children's Hospital (Centers 2–4). The final PK model was developed based on 2491 busulfan plasma concentration–time measurements from 379 dosing days. One patient had eight blood samples that were compromised and needed to be discarded, one patient had only 3 days of treatment, and another had an additional dose administered to total 5 days of therapy. The demographic and clinical characteristics of the subjects together with the busulfan exposure targets and estimated busulfan exposure attained using NCA are summarized in Table 2.

TABLE 2

Demographics and clinical features of pediatric hematopoietic stem cell transplant patients involved in population non‐linear mixed effect (NLME) model development (n = 95)

	Nonmalignant (N = 15)	Malignant (N = 80)	Total (N = 95)
Sex, n (%)
Male	13 (86.7)	36 (45.0)	49 (51.6)
Female	2 (13.3)	44 (55.0)	46 (48.4)
Age, years
Median [lower 2.5%, upper 97.5%]	2.80 [0.375, 7.49]	4.50 [1.40, 17.3]	4.20 [0.735, 17.2]
Actual weight, kg
Median [lower 2.5%, upper 97.5%]	13.8 [6.03, 26.3]	17.6 [10.0, 85.5]	17.0 [7.77, 83.3]
Body mass index, kg/m2
Median [lower 2.5%, upper 97.5%]	17.5 [13.8, 19.17]	18.4 [13.39, 30.64]	18.2 [13.35, 32.39]
Treatment center, n (%)
Center 1	14 (93.3)	66 (82.5)	80 (84.2)
Centers 2–4	1 (6.7)	14 (17.5)	15 (15.8)
Chemotherapy regimen, n (%)
Bu/Flu	12 (80.0)	12 (15.0)	24 (25.3)
Bu/Flu/Mel	0 (0)	19 (23.8)	19 (20.0)
Bu/Cy	0 (0)	5 (6.2)	5 (5.3)
Bu/Flu/TT	2 (13.3)	22 (27.5)	24 (25.3)
Bu/Cy/TT	0 (0)	0 (0)	0 (0)
Bu/Mel	0 (0)	21 (26.2)	21 (22.1)
Other	1 (6.7)	1 (1.2)	2 (2.1)
Seizure prophylaxis, n (%)
Levetiracetam	13 (86.7)	58 (72.5)	71 (74.7)
Benzodiazepine	1 (6.7)	22 (27.5)	23 (24.2)
Other	0 (0)	0 (0)	0 (0)
Missing	1 (6.7)	0 (0)	1 (1.1)
Immunotherapy, n (%)
None	1 (6.7)	43 (53.8)	44 (46.3)
ATGAM®	0 (0)	5 (6.3)	5 (5.3)
Thymoglobulin®	11 (73.3)	1 (1.3)	12 (12.6)
Grafalon®	1 (6.7)	25 (31.3)	26 (27.4)
Alemtuzumab	1 (6.7)	0 (0)	1 (1.1)
Missing	1 (6.7)	6 (7.5)	7 (7.4)
Metronidazole, n (%)
No	14 (93.3)	66 (82.5)	80 (84.2)
Yes	0 (0)	3 (3.8)	3 (3.2)
Missing	1 (6.7)	11 (13.8)	12 (12.6)
Acetaminophen, n (%)
No	3 (20.0)	35 (43.8)	38 (40.0)
Yes	12 (80.0)	36 (45.0)	48 (50.5)
Missing	0 (0)	9 (11.2)	9 (9.5)
Posaconazole, n (%)
No	14 (93.3)	66 (82.5)	80 (84.2)
Yes	0 (0)	4 (5.0)	4 (4.2)
Missing	1 (6.7)	10 (12.5)	11 (11.6)
Voriconazole, n (%)
No	14 (93.3)	66 (82.5)	80 (84.2)
Yes	0 (0)	4 (5.0)	4 (4.2)
Missing	1 (6.7)	10 (12.5)	11 (11.6)
Fluconazole, n (%)
No	4 (26.7)	42 (52.5)	46 (48.4)
Yes	10 (66.7)	28 (35.0)	38 (40.0)
Missing	1 (6.7)	10 (12.5)	11 (11.6)
AUCcum target for individual, mg · h/L
Mean (SD)	82.7 (10.7)	79.3 (5.48)	79.8 (6.62)
Median [lower 2.5%, upper 97.5%]	90.0 [65.1, 90]	75.0 [65.3, 90]	78.0 [65.5, 90]
AUCcum attained for individual, mg · h/L
Mean (SD)	80.2 (10.2)	77.8 (7.29)	78.2 (7.82)
Median [lower 2.5%, upper 97.5%]	81.7 [63.2, 91.4]	77.7 [63.1, 91.5]	78.0 [63.3, 91.3]
Missing, n (%)	0 (0)	1 (1.2)	1 (1.1)

Abbreviations: AUCcum, cumulative area under the concentration‐time curve following all doses; Bu, busulfan; Cy, cyclophosphamide; Flu, fludarabine; Mel, melphalan; SD, standard deviation; SOS, sinusoidal obstructive syndrome; TT, thiotepa.

PK model

A two‐compartment model with first‐order elimination best described the concentration‐time data (change in OFV [ΔOFV] −140.3 compared with a one‐compartment model). A combined residual error was superior to a proportional error (ΔOFV −652 compared with a proportional error model). Estimating DI was superior to using the recorded duration listed in the medical notes (ΔOFV −43). DI was later fixed to 2.43 h and confirmed using a sensitivity analysis. Addition of time‐associated CL based on a continuous model (Equation 1) improved the model fit further (ΔOFV −493.2). In terms of the stochastic model, inclusion of IIV on CL and V1 with correlation (ΔOFV −232.2), inclusion of IOV on CL (ΔOFV −1147), and inclusion of IOV on V1 (ΔOFV −156.3) all improved the model fit.

During covariate screening, inclusion of allometrically scaled NFM on CL, Q, V1, and V2 (Equation 2) together with a maturation function on CL (Equation 3) significantly improved the fit of the model (ΔOFV −108.8 and ΔOFV −6.06, respectively) and reduced IIV on CL from 43% to 14.8% and IIV on V1 from 69.2% to 36%. NFM was a better descriptor of size than TBW; however, due to limited numbers of overweight and obese patients, the fraction of fat mass contributing to NFM on PK parameters for busulfan was fixed a priori (Table 3) as per McCune and colleagues. The Hill coefficient within the maturation equation for CL was estimated as 2.43 and later fixed to 2.3 in line with literature. , When both TM50 and Hill parameters were estimated for the maturation component, the model failed to converge. A simplified time‐associated CL equation involving the γ term fixed to 1 showed no worsening of model fit and was hence adopted for further covariate testing. No concomitant medication had a significant effect on busulfan CL during covariate screening. A trend toward a difference in CL, T 50,i and ΔCLmax,i with the use of an azole antifungal was observed but requires further investigation with additional subjects and hence was not included in the final model. There was no difference in CL observed in subjects with and without malignant disease.

TABLE 3

Parameter estimates for the final population pharmacokinetics model together with the bootstrap analysis

Parameter	Final model			Bootstrap analysis
	Estimate	RSE (%)	Shrinkage (%)	Median	5th CI	95th CI
CL (L/h/62 kg)	14.5	5.86	–	14.7	13.1	16.9
V1 (L/59 kg)	40.6	4.42	–	40.5	37.3	43.7
Q (L/h/56.1 kg)	1.92	22.9	–	2.02	1.23	3.17
V2 (L/59 kg)	3.57	15.9	–	3.70	2.73	4.81
DI (h)	2.43 (fixed)	–	–	2.43 (fixed)	–	–
ΔCLmax	−0.198	16.9	–	−0.193	−0.283	−0.146
T 50 (h)	50.6	34.6	–	49.6	24.1	106
IIV on CL (CV%)	14.7%	12.4	1.96	14.3%	10.9%	17.6%
IIV on V1 (CV%)	34.9%	20.1	0.997	34.5%	20.6%	46.0%
IOV on CL (CV%)	6.61%	8.77	13.9	6.47%	5.35%	7.53%
IOV on V1 (CV%)	9.71%	9.61	29.8	9.64%	8.11%	11.2%
Prop RUV (CV%)	24.3%	5.12	–	24.3%	23.2%	25.3%
Add RUV (mg/L)	0.0300	22.9	–	0.0297	0.0182	0.0405
Ffat(CL)	0.509 (fixed)	–	–	0.509 (fixed)	–	–
Ffat(V1)	0.203 (fixed)	–	–	0.203 (fixed)	–	–
Ffat(Q)	0 (fixed)	–	–	0 (fixed)	–	–
Ffat(V2)	0.203 (fixed)	–	–	0.203 (fixed)	–	–
Hill (maturation)	2.3 (fixed)	–	–	2.3 (fixed)	–	–
TM50 (maturation)	45.6 (fixed)	–	–	45.6 (fixed)	–	–

Note: CL, Q, V1, and V2 parameter estimates were allometrically scaled using a standard NFM of 62 kg for CL, 56.1 kg for Q, and 59 kg for V1 and V2 corresponding to an allometrically scaled total body weight of 70 kg. Coefficient of Variation (CV%) are calculated as the Square root of variance (OMEGA from NONMEM®) × 100; RSE of parameter estimates are calculated as 100 × (SE/typical value); RSE of between‐subject variability magnitude are calculated as 100 × (SE/variance estimate)/2. Proportional RUV is presented as standard deviation. Shrinkage (%) is calculated as 100 × (1 − SD (ETA from NONMEM®)/sqrt(variance)). Overall residual variability shrinkage was estimated to be 14.5%. The correlation coefficient between CL and V1 was estimated as 0.0295. Hill = the steepness of the function within CLmat equation.

Abbreviations: ΔCLmax, maximal possible change in CL relative to baseline for the individual; Add RUV, additive residual unexplained variability; CI, confidence interval; CL, clearance; DI, duration of infusion; Ffat(CL), the fraction of fat mass (Ffat) contributing to NFM for CL parameter; Ffat(V1), the fraction of fat mass (Ffat) contributing to NFM for V1; Ffat(Q), the fraction of fat mass (Ffat) contributing to NFM for Q; Ffat (V2), the fraction of fat mass (Ffat) contributing to NFM for V2; IIV, interindividual variation; IOV, intraoccasion variation; NFM, normal fat mass; Prop RUV, proportional residual unexplained variability; Q, intercompartmental clearance; RSE, relative standard error; T 50, time at which 50% of ΔCLmax is attained; TM50, Post menstrual age in weeks where the CL of busulfan is considered to be 50% of the adult value; V1, volume of distribution central compartment; V2, volume of distribution peripheral compartment.

FINAL MODEL

Central compartment

Peripheral compartment

Parameters

where dc/dt = change in concentration over time, CONC = concentration busulfan in central compartment (V1) in mg/L, PERI = concentration busulfan in peripheral compartment (V2) in mg/L, and Time1 = time in h since beginning of infusion of Dose 1, CLsize = Equation 2 for CL, CLmat = Equation 3.

Population PK parameter estimates for the final model are summarized in Table 3. Standard errors for the final model parameter estimates showed that parameters were well estimated and were below 25% for all except for T 50, which had a relative standard error (RSE) of 35%. The pcVPCs for the final model showed good agreement between the observed and predicted data across all days (Figure 1). Results from the nonparametric bootstrap analysis demonstrated that final model estimates were robust (Table 3). GOF plots are shown in Figure S1, demonstrating adequate model fit.

FIGURE 1

Prediction‐corrected visual predictive check plots based on 1000 simulations from the final population PK model across the treatment course. Observed data (blue dots), median (red line), 95th and 5th percentiles of observed data (red dashed line), median (black line), 95th and 5th percentiles of simulated data (black dashed line), and 95% confidence interval of simulated data (gray area).

Simulations: six dosing scenarios for the virtual population

Table 4 shows the predicted AUCcum for the same virtual population of 1000 simulated subjects for each of the Scenarios 1–6. When an AUC0–24 was estimated following D1 only using NCA to determine the next doses (Doses 2–4) using the equation from the product information (NCA_PI_D1), a median AUCcum of 129 mg · h/L (range, 105–241) (Table 4, Figure 2) was obtained and 100% of patients were above both the individualized target AUCcum of 90 mg · h/L ±5% and the target range of 78–101 mg · h/L (Table 4). When sampling was performed and considered following each dose (Scenario NCA_PI_D1‐4) the median AUCcum of 125 mg · h/L (range, 100–231) remained high, and the proportion of subjects achieving target was considerably reduced compared with scenarios using model‐based exposure estimates (MOD_PI_D1‐4 and MOD_MOD_D1‐4). Under Scenario NCA_PI_D1‐4, all subjects had above target exposure (Table 4). When AUC0–24 was estimated using model‐based methods following D1 only to calculate the next doses (Doses 2–4) from the equation within the product information (MOD_PI_D1), median AUCcum achieved was lower compared with Scenario NCA_PI_D1 at 96.9 mg· h/L (range, 95.7–98.6) (Table 4, Figure 2), with 100% of patients achieving exposure within the optimal exposure target range (78–101 mg · h/L; Table 4, Figure 2) but 100% above target exposure of 90 mg · h/L ±5%. When sampling was performed following each dose (MOD_PI_D1‐4), the median AUCcum achieved was 90.5 mg · h/L (range, 90.4–91.2), and all subjects achieved target exposure (100% within 90 mg · h/L ±5%) as seen in Figure 2. When model‐based AUC0–24 and individual CL estimates (immediately prior to Dose 2) were used to calculate Doses 2–4 from sampling following D1 only (MOD_MOD_D1), the median AUCcum achieved was 92.4 mg · h/L (range, 91.9–92.9). When sampling was performed following D1 only, this scenario performed best with all patients attaining the target AUCcum of 90 mg · h /L ±5% (Table 4). When samples were taken with each dose allowing subsequent dose adjustments (MOD_MOD_D1‐4), the median AUCcum was 90 mg · h/L (range, 90–90.8), with all subjects within the target range (Figure 2).

TABLE 4

Results of simulation AUCcum targets achieved and doses administered based on six different dose‐adjustment scenarios (N = 1000 virtual simulated subjects)

	NCA_PI_D1 (N = 1000)	NCA_PI_D1‐4 (N = 1000)	MOD_PI_D1 (N = 1000)	MOD_PI_D1‐4 (N = 1000)	MOD_MOD_D1 (N = 1000)	MOD_MOD_D1‐4 (N = 1000)
Subjects who achieved AUC(CUM) targets, n (%)
AUCcum achieve 90 mg · h/L ±5%
Below	0 (0)	0 (0)	0 (0)	0 (0)	0 (0)	0 (0)
Within	0 (0)	0 (0)	0 (0)	1000 (100)	1000 (100)	1000 (100)
Above	1000 (100)	1000 (100)	1000 (100)	0 (0)	0 (0)	0 (0)
AUCcum achieve within 78–101 mg · h/L
Below	0 (0)	0 (0)	0 (0)	0 (100)	0 (0)	0 (0)
Within	0 (0)	2 (0.2)	1000 (100)	1000 (100)	1000 (100)	1000 (100)
Above	1000 (100)	998 (99.8)	0 (0)	0 (0)	0 (0)	0 (0)
Dose administered each day (mg)
Dose 1
Mean (SD)	98.5 (52.3)	98.5 (52.3)	98.5 (52.3)	98.5 (52.3)	98.5 (52.3)	98.5 (52.3)
Median [min, max]	92.1 [21.2, 351]	92.1 [21.2, 351]	92.1 [21.2, 351]	92.1 [21.2, 351]	92.1 [21.2, 351]	92.1 [21.2, 351]
Dose 2
Mean (SD)	232 (122)	232 (122)	162 (88.3)	162 (88.3)	153 (83.4)	153 (83.4)
Median [min, max]	208 [50.4, 1280]	208 [50.4, 1280]	145 [35.7, 654]	145 [35.7, 654]	137 [33.7, 616]	137 [33.7, 616]
Dose 3
Mean (SD)	232 (122)	87.0 (46.7)	162 (88.3)	148 (80.6)	153 (83.4)	148 (80.6)
Median [min, max]	208 [50.4, 1280]	81.0 [18.5, 314]	145 [35.7, 654]	132 [32.6, 595]	137 [33.7, 616]	132 [32.6, 595]
Dose 4
Mean (SD)	232 (122)	345 (191)	162 (88.3)	139 (75.5)	153 (83.4)	145 (78.8)
Median [min, max]	208 [50.4, 1280]	309 [74.0, 2130]	145 [35.7, 654]	124 [30.5, 557]	137 [33.7, 616]	129 [31.9, 583]
AUCcum following each dose (mg · h/L)
AUCcum h24 (mg · h/L)
Mean (SD)	15.4 (2.59)	15.4 (2.59)	15.4 (2.59)	15.4 (2.59)	15.4 (2.59)	15.4 (2.59)
Median [min, max]	15.1 [9.05, 28.5]	15.1 [9.05, 28.5]	15.1 [9.05, 28.5]	15.1 [9.05, 28.5]	15.1 [9.05, 28.5]	15.1 [9.05, 28.5]
AUCcum h48 (mg · h/L)
Mean (SD)	53.9 (6.51)	53.9 (6.51)	41.8 (1.68)	41.8 (1.68)	40.3 (1.73)	40.3 (1.73)
Median [min, max]	52.6 [42.6, 88.5]	52.6 [42.6, 88.5]	41.6 [37.7, 50.3]	41.6 [37.7, 50.3]	40.1 [36.1, 49.1]	40.1 [36.1, 49.1]
AUCcum h72 (mg · h/L)
Mean (SD)	93.8 (12.9)	68.8 (7.28)	69.1 (0.728)	66.7 (0.810)	66.1 (0.836)	65.2 (0.864)
Median [min, max]	90.6 [73.6, 164]	67.8 [53.0, 107]	69.0 [67.3, 72.7]	66.6 [64.7, 70.8]	66.0 [64.0, 70.3]	65.1 [63.1, 69.6]
AUCcum h96 (mg · h/L)
Mean (SD)	134 (19.6)	129 (18.6)	96.9 (0.271)	90.5 (0.0569)	92.4 (0.104)	90.1 (0.0633)
Median [min, max]	129 [105, 241]	125 [100, 231]	96.9 [95.7, 98.6]	90.5 [90.4, 91.2]	92.4 [91.9, 92.9]	90.0 [90.0, 90.8]

Abbreviations: AUCcum, cumulative exposure of busulfan over the entire course (4 days); h, hour; L, liter; mg, milligram; max, maximum; min, minimum; SD, standard deviation.

Scenario descriptions are in Table 1.

FIGURE 2

Boxplot showing cumulative area under the concentration‐time curve following all doses (AUCcum) of four once‐daily doses of busulfan. Black dot‐dash line = target AUCcum of 90 mg · h/L, gray‐shaded area = target AUCcum of 90 mg · h/L ±5% (within 85.5–94.5 mg · h/L), gray dashed lines = optimal target range described by Bartelink et al. (within 78–101 mg · h/L).* Scenario descriptions can be found in Table 1 and outlying values not shown in the figure. CL , individual model‐based clearance calculated at the end of the dosing interval for the current dose; Model, final model as per article; NCA, noncompartmental analysis; PI (proportional), proportional equation according to product information leaflet.

Scenarios using NCA to estimate exposure (NCA_PI_D1 and NCA_PI_D1‐4) were associated with the highest mean and median doses administered, the highest estimated daily AUC0–24 and AUCcum, and all patients were above target exposure. When model‐based exposure estimates were used, all patients achieved busulfan exposure within the optimal target range of 78–101 mg · h/L. In these four scenarios, using the proportional method to calculate next dose for Doses 2–4 (MOD_PI_D1) with sampling performed following D1 only led to more subjects attaining >5% above target AUCcum compared with when model‐based CL estimates were used (MOD_MOD_D1 [100% vs. 0%; Table 4]). Higher doses were suggested when using the proportional dose equation compared with model‐based CL (Table 4), corresponding to higher AUC0–24 exposure, particularly following Dose 2 (Table 4, Figure S2). Scenarios MOD_PI_D1‐4, MOD_MOD_D1, and MOD_MOD_D1‐4 were associated with smaller dose adjustments and reduced doses during the 4‐day treatment course compared with Scenarios NCA_PI_D1, NCA_PI_D1‐4, and MOD_PI_D1.

DISCUSSION

This is the largest study to date to characterize the PK of busulfan based on data from all 4 days of treatment in pediatric HSCT subjects. The final population PK model adequately predicted individual busulfan concentrations. Busulfan demonstrated a time‐associated reduction in CL during the course of treatment. The average reduction in CL was 11.6%, with 8.1% reduction occurring within 48 h following initiation of therapy. This finding is in keeping with other studies that have reported reductions in CL of between 8.1% and 20% , , , , , , (Table S1). RSE on T 50 was 35%, and bootstrap results had a large range. This parameter may have some predictable variation between patients that is yet to be described and could be related to concomitant azole administration. Azole antifungals each have different potential to induce and/or inhibit different cytochrome P450 enzyme subfamilies as well as compete for metabolism as a substrate, potentially impacting metabolism pathways for busulfan.

An empirical Emax model was explored to describe time‐associated busulfan CL as it allows for a reduction in CL over time without CL achieving zero or negative values, which are physiologically implausible. This model has been used previously for busulfan and described the data well. The biological rationale for the observed reduction in busulfan CL with continued therapy requires further exploration. It is possible that busulfan metabolites (tetrahydrothiophene and γ‐glutamyldehydroalanylglycine) may interfere with busulfan biotransformation , , , , , or alternatively glutathione depletion may play a role. , NFM was the preferred size descriptor, which is in keeping with one of the largest published NLME models to date. The PK parameters in the final model were similar to previously published values as reported in a recent review. Typical population CL was estimated to be 14.5 L/h/62 kg (IIV 14.8%; IOV 6.6%), which is slightly higher compared with the typical CL of 11.4 L/h/62 kg from the published model by McCune et al. using allometrically scaled NFM on parameters and a maturation component on CL.

Using dose‐adjustment strategies that assume CL remains the same during the course of treatment (equation from PI using proportional adjustment) could result in a higher‐than‐expected AUC0–24 following Doses 2–4 and thus a larger AUCcum (Table 4, Figure 2). This can be mitigated by first using a model‐based estimation of AUC0–24, which incorporates time‐associated CL estimates (Scenarios MOD_PI_D1 and MOD_PI_D1‐4) and, second, by performing sampling after each dose (MOD_PI_D1‐4). Simulations comparing the method for estimating AUC0–24 highlight that NCA underestimates AUC0–24, resulting in an adjustment of doses higher than required (comparing Scenarios NCA_PI_D1 and NCA_PI_D1‐4 using NCA calculated exposure with Scenarios MOD_PI_D1 and MOD_PI_D1‐4 using model‐based exposure). The use of daily sampling does not mitigate the impact of using NCA‐based AUC0–24 estimates for dose adjustments with 100% of subjects still attaining actual AUCcum estimates above target AUCcum (Figure 2). It has been reported that performing dose adjustments to target a certain exposure increases the risk of toxicities for subjects, , and this is likely due to the use of practices in line with Scenarios NCA_PI_D1 and NCA_PI_D1‐4 using NCA exposure estimates when the targets in the literature have been estimated using model‐based methods. Therefore, it is imperative that the method for exposure estimation is considered when determining the exposure target. , The method used to estimate exposure for busulfan and methods used to calculate dose adjustments should be listed and standardized in publications to allow for the correct interpretation of results and implementation into clinical practice.

Scenarios MOD_MOD_D1 and MOD_MOD_D1‐4 used the CL estimated using the final model to guide dose adjustments. When sampling can only be performed following D1, using both model‐based AUC0–24 estimation and CL (Scenario MOD_MOD_D1) resulted in 100% of subjects attaining the target AUCcum ±5% (compared with 0% for Scenarios NCA_PI_D1 and MOD_PI_D1). Whereas when daily sampling and dose adjustment are available, model‐based AUC0–24 with dose adjustment using either method resulted in 100% of subjects within the target AUCcum ±5%. The success of Scenario MOD_MOD_D1 highlights the benefits of using model‐based methods for both estimating exposure and calculating next doses. Studies involving reduced sampling scenarios and model‐based methods using Bayesian forecasting software have reported adequate precision and accuracy when estimating AUCcum. It is clear that model‐based estimation of exposure is important to ensure subjects attain target AUCcum. Use of a population PK model that allows for the reduction of CL during the course of treatment within a Bayesian forecasting software program should be explored in future prospective trials (Scenario MOD_MOD_D1‐4; Figure 2).

This study has some limitations. Residual unexplained variability and IOV on both CL and V1, although low (6.6% and 9.7% respectively), were not implemented in the simulations. The data set was collected across four separate hospitals, introducing variable practice in busulfan administration, sample collection, and record keeping. The number of subjects on certain concomitant medications was low, limiting analysis opportunities. The study included retrospective data that relies on the accuracy of documentation.

Time‐associated reduction in CL during a typical 4‐day course of once‐daily i.v. busulfan was described. It is strongly recommended that the antiquated practice of determining busulfan exposure using NCA be retired and replaced by model‐based methods using numeric integration, preferably implementing a model that describes a change in CL across the dosing period. Future studies should focus on the use of dose‐adjustment strategies, such as that represented in Scenarios MOD_MOD_D1‐4 and MOD_PI_D1‐4 within Bayesian forecasting software to increase success in achieving individual AUCcum exposure targets and desired patient outcomes. Sample collection following each dose increases the numbers of patients achieving the cumulative exposure target; however, limited sampling strategies could be useful when model‐based methods are implemented for both estimation of exposure and calculating next doses. Lastly, the method used to estimate exposure for outcome analysis and used for dose adjustments in patients should be considered an essential component for reporting in future articles. Standardization of methods for busulfan exposure monitoring and dose adjustment would allow for improved interpretation of results from the literature and implementation into clinical practice.

AUTHOR CONTRIBUTIONS

R.L., S.H., C.E.S., C.J.F., S.R., L.T., R.M., and T.O. wrote the manuscript. R.L., S.H., C.E.S., C.J.F., S.R., L.T., R.M., and T.O. designed the research. R.L. performed the research. R.L. and S.H. analyzed the data.

CONFLICT OF INTEREST

The authors declared no competing interests for this work. As an Associate Editor for CPT: Pharmacometrics and Systems Pharmacology, Stefanie Hennig was not involved in the review or decision process for this article.

Appendix S1:

Click here for additional data file.