Semi-mechanistic population pharmacokinetic model incorporating glutathione S-transferase activity for personalized busulfan dosing in pediatric allogeneic hematopoietic cell transplantation

Background: Busulfan is known for its high inter- and intra-individual pharmacokinetics/pharmacodynamics (PK/PD) variability, especially in children. Therefore, we aimed to identify factors affecting PK variability of busulfan in pediatric allogeneic hematopoietic cell transplantation (HCT) recipients and investigate the effect of glutathione S-transferase (GST) activity on busulfan metabolism using a semi-mechanistic population PK model.

Methods: Overall, 636 whole-blood busulfan concentrations from 65 pediatric HCT recipients were analyzed using nonlinear mixed-effects modeling. A semi-mechanistic population PK model was developed to describe busulfan metabolism in response to glutathione (GSH) depletion. The effects of potential covariates were selected based on previous study and physiologically-based theoretical mechanisms. Virtual clinical trials were conducted to compare different dosing strategies, and model-based optimal dosing regimen was recommended.

Results: A two-compartment model with first-order absorption was selected to describe busulfan PK. A GSH compartment was added to represent the relative amount of GSH available at any time. The estimated mean clearance of busulfan was 9.57 L h⁻¹ (relative standard error: 10.8%). Busulfan disposition was best described by including normal fat mass (NFM) allometrically and GST enzyme activity on S_GSH exponentially. The S_GSH increased by 40.6% as GST enzyme activity increased from 0.9 nmol/min/mL to 20.7 nmol/min/mL. Patients with weights (WT) of 9–16 kg are at high risk of sinusoidal obstructive syndrome (SOS) when receiving WT-based dosing strategy.

Conclusion: NFM, age-dependent maturation function, and GST enzyme activity may contribute to busulfan PK variability. The WT-based dosing strategy showed a higher risk of SOS than the age-based dosing strategy in 9–16 kg patients.

1 Introduction

Busulfan, a bifunctional DNA-alkylating agent, is widely applied as a chemotherapeutic in combination with cyclophosphamide, cytarabine, and fludarabine before allogeneic hematopoietic cell transplantation (HCT) (Chen et al., 2018; Lawson et al., 2021). This treatment can reduce the immune response to avoid graft rejection and provide favorable conditions for donor cell engraftment. Critically, subtherapeutic drug exposure levels correlate with increased relapse rates or graft failure, while supratherapeutic concentrations are linked to a greater risk of severe toxicities and treatment-related mortality (Bartelink et al., 2016).

The clinical application of busulfan is complicated by its high inter- and intra-individual pharmacokinetics/pharmacodynamics (PK/PD) variability, particularly in children (Lawson et al., 2021). Consequently, intravenous (IV) administration is preferred in children because of the higher bioavailability and reduced PK variability compared to oral formulations (Palmer et al., 2016). Following IV infusion, busulfan undergoes rapid distribution and binds extensively to erythrocytes (approximately 47%) and plasma proteins (approximately 32%). Hepatic metabolism occurs primarily through conjugation with glutathione (GSH) mainly through glutathione S-transferases (GSTs) (Lawson et al., 2021; Scian et al., 2016), with renal excretion playing a minor role, only about 2% of busulfan is detected unmetabolized in the urine (Hassan et al., 1989).

Furthermore, owing to its narrow therapeutic index and large PK/PD variability, administering an initial busulfan IV dose based only on body weight (WT) may result in failure to reach the target therapeutic window (Ben Hassine et al., 2021; Huang et al., 2022). Crucially, clinical evidence demonstrates that PK-guided IV busulfan dosing is superior to body-size dosing in patients with myeloid leukemia and myelodysplastic syndrome, yielding reduced relapse, transplant-related mortality, and overall hazard ratio (Andersson et al., 2017). Given these limitations of weight-based dosing and the demonstrated superiority of personalized approaches, therapeutic drug monitoring (TDM) is recommended as the standard of care for optimizing individual regimens (Palmer et al., 2016).

Model-informed precision dosing utilizes population PK (popPK) models combined with maximum posterior Bayesian estimation to optimize both initial and subsequent dosing regimens based on TDM measurements (Darwich et al., 2017; Shukla et al., 2020). Currently, over 40 popPK models have been developed to characterize IV busulfan PK in pediatric patients. Among the covariates, body size, age, GST alpha 1 (GSTA1) genetic variations, and dosing schedule (day/time) are the most well-documented factors contributing to busulfan clearance variability (Lawson et al., 2021; Takahashi et al., 2023).

Body size descriptors, including body surface area (BSA), fat-free mass (FFM), and normal fat mass (NFM), significantly influence busulfan PK in pediatric patients (Huang et al., 2022; Lawson et al., 2021; Takahashi et al., 2023). Notably, most of these descriptors are typically incorporated into pediatric busulfan dosing individualization via allometric scaling, an approach grounded in fractal geometry principles and cross-species biological patterns (Anderson and Holford, 2008; Trame et al., 2011). Accurate quantification of allometric exponents requires data spanning the full maturation spectrum from neonates to adults (Gonzalez-Sales et al., 2022). Consequently, based on physiologically-based descriptions of body composition and theory-based allometric principles, Du et al. estimated the clearance (CL) for busulfan through allometry NFM, a maturation fraction (F_mat), and distribution volume (V) based on FFM (Du et al., 2022).

Time-varying CL was observed over a 4-day treatment with an every-6-h dosing regimen of busulfan (Lawson et al., 2021; Takahashi et al., 2023). Specifically, CL demonstrated a progressive decline of 8.1%–20% across treatment cycles compared to baseline (Day 1), making it challenging to obtain the desired busulfan target exposure (Lawson et al., 2021). To explain this nonlinear elimination, the empirical Michaelis-Menten equation and semi-mechanistic enzyme depletion model have been employed (Langenhorst et al., 2020; Long-Boyle et al., 2015). Central to this phenomenon, busulfan-GSH conjugate serve as main intermediate metabolite, with baseline GSH levels correlating with busulfan CL; therefore, Langenhorst et al. hypothesized that busulfan-mediated GSH depletion causes nonlinear elimination (Langenhorst et al., 2020). However, the GST enzyme activity was not considered in their model.

Genetic polymorphisms in GSTA1 are associated with 8%–27% reduction in CL (Huang et al., 2022). However, GST expression exhibits complex regulation beyond genetics, demonstrating age- and sex-dependent variations (Hines, 2008; Miyagi et al., 2009; Ten Brink et al., 2013). During pediatric development, age modulates hepatic enzyme maturation, serum protein concentrations, and body composition (water-to-fat ratio), while weight correlates with somatic growth and governs hepatic blood flow dynamics. These parameters jointly determine the evolving liver-to-body mass ratio, a critical determinant of drug metabolism capacity. GST activity decreases from infancy to early adolescence, and developmental differences in activity can markedly alter drug disposition (Gibbs et al., 1999; Hall et al., 1999).

Given this intricate interplay of physiological and pharmacological factors, comprehensive understanding of busulfan PK characteristics becomes essential for target exposure attainment. To address this, our study employs a semi-mechanistic popPK model to quantify sources of PK variability in IV busulfan exposure among pediatric HCT recipients, and mechanistically characterize GST-mediated metabolic pathways influencing drug disposition. In addition, virtual clinical trials were conducted to compare different dosing strategies, and a model-based optimal dosing regimen was recommended.

2 Materials and methods

2.1 Patients and data collection

Data were prospectively collected from 65 pediatric HCT recipients who underwent bone marrow transplantation after receiving busulfan IV during preparative chemotherapy at the Children’s Hospital of Fudan University. All patients were administered 0.8–1.2 mg/kg of busulfan via 2 h IV infusion every 6 h, depending on the patient’s WT. Patients with normal organ function were included, and those with unavailable busulfan PK data owning to difficulties in blood sampling were excluded. Demographic and pathophysiological data were prospectively obtained during routine clinical visits between August 2020 and November 2021. This study was approved by the Ethics Committee of the Children’s Hospital (ethics approval number: 2020-271) and conducted in accordance with the Declaration of Helsinki. Notably, all patients and their parents provided written informed consent to participate prior to enrolment in the study.

Overall, 636 whole-blood busulfan concentrations were available for model analysis. All patients received 12 doses of IV busulfan. Samples were obtained 2, 2.5, 3, 4, and 6 h following the infusion of dose 1, and pre-dose concentrations (C₀) were collected before doses 6 and 12. Additionally, to balance the blood capacity taken and the sampling of terminal elimination, samples were collected at 2 h, 4 h, 8 h (Group A, 33 patients) and 2 h, 6 h, 12 h (Group B, 32 patients) after the infusion of dose12 (Figure 1). Furthermore, 1 mL of whole blood was collected in EDTA tubes for each sample, and all samples were stored at −80°C until analysis.

Figure 1

Figure 1. Diagram showing the prospective study dataset and trial design.

2.2 Determination of busulfan concentration and GST enzyme activity

Quantification of busulfan plasma concentrations was performed by a validated liquid chromatography/mass spectrometry. The assay demonstrated linearity across the analytical range of 10–5,000 ng/mL, with a lower limit of detection of 10 ng/mL. Additionally, blood samples were collected for the first time using micro-quartz colorimetry to determine GST enzyme activity. GST can catalyze the binding of GSH to 1-chlorom-2,4-ditrobenzene, which can be detected at a wavelength of 340 nm. Furthermore, 20 μL serum was mixed with the detection reagent and detected twice, before and after a 5 min water bath.

2.3 Semi-mechanistic population pharmacokinetic modeling

The popPK model was established using a nonlinear mixed-effects modeling approach implemented in NONMEM® (version 7.4; ICON Development Solutions, Ellicott City, MD, United States), with Pirana 2.9 serving as the interface for Perl Speaks NONMEM (PsN; version 4.9.0) to streamline model diagnostics and bootstrapping (Keizer et al., 2013). Graphical analyses were conducted through R software (version 3.5.0; http://www.r-project.org/). The first-order conditional estimation method, including η-ε interactions (FOCE-I), was employed throughout the method-building procedure (Beal et al., 1989).

The busulfan PK profile was best characterized by a two-compartment structural model with first-order elimination kinetics. Primary estimated parameters included CL, central volume of distribution (V_c), inter-compartmental clearance (Q), and peripheral volume of distribution (V_p). Variability components were systematically quantified through between-subject variability (BSV), inter-occasion variability (IOV), and residual unexplained variability (RUV). BSV modeled via log-normal distributions for all parameters, except Q. However, IOV was assumed to be the same across dosing occasions (Karlsson and Sheiner, 1993).

Demographic and disease-specific pathophysiological indices, and concomitant medications (Table 1) were systematically screened for potential covariates. Body size is the most identified covariate in busulfan PK modeling; therefore, four body size metrics [WT, BSA, FFM, and NFM (Text S1)] were used to determine the most suitable body size descriptor (Lawson et al., 2021). The variabilities in CL and V were characterized allometrically using body size and composition (Anderson and Holford, 2009; West et al., 1997), while busulfan metabolism maturation upon CL was evaluated using an empirical sigmoid function (F_mat, Equation 1) (Gonzalez-Sales et al., 2022). Post-menstrual age (PMA) was a composite developmental biomarker integrating both gestational age and post-natal age.

$F_{mat}=\left(\frac{{}_1}{1+{\left(\frac{PMA}{{TM}_{50}}\right)}^{-\text{Hill}}}\right)(1)$

where TM50 is the PMA at which maturation achieving 50% of the adult value, and Hill defines the steepness of the sigmoid decline.

Table 1

Table 1. Patients demographics used to develop and evaluate population pharmacokinetic model.

Covariate selection was conducted through a stepwise approach (Beal et al., 1989). The influence of continuous covariates was evaluated through linear, exponential, and power function models. For categorical variables (e.g., concomitant medications), intergroup comparisons were performed by analyzing fractional change differences. The variability between dosing regimen cycles in the time-dependent CL of busulfan was estimated using a linear function model and IOV on CL.

Therefore, to investigate the influence of GST enzyme activity on busulfan metabolism, the empirical model developed was used as a base, with a dedicated compartment incorporated to dynamically quantify the relative amount of GSH available over time, based on theoretical mechanisms, as reported by Langenhorst et al. (2020) (Figure 2).

Figure 2

Figure 2. Busulfan semi-mechanistic population pharmacokinetic model structure. CL, clearance; Q, inter-compartmental clearance; V₁, central compartment; V₂, peripheral compartment. Dashed lines indicate the conjugation of busulfan metabolism and glutathione.

The GSH compartment was initialized with a baseline normalized value of 1, and the zero-order synthesis rate was constrained to equal the first-order elimination rate constant at equilibrium, ensuring mass balance. Busulfan metabolism was modeled as a GSH-dependent conjugation process, with the scaling parameter S_GSH quantifying the proportionality between busulfan metabolism and corresponding GSH depletion (Equation 2).

$\frac{{dA}_{GSH}}{dt}=\frac{S_{GSH}}{V_1}\times A_{GSH}\times k_{10}\times A_{{bu}_1}(2)$

Where $A_{{bu}_1}$ represents the amount of busulfan in the central compartment, $A_{GSH}$ represents the amount of GSH in the theoretical compartment, $V_1$ represents the central volume of distribution, and $k_{10}$ represents the busulfan elimination constant.

Subsequently, GST enzyme activity was tested as a continuous covariate on S_GSH.

The visual model fit was evaluated using standard goodness-of-fit (GOF) criteria, reductions in the objective function value (OFV) for nested models, Akaike information criteria (AIC) and Bayesian information criteria (BIC) for non-nested models, and acceptable precision of estimates (Beal et al., 1989; Donohue et al., 2011). Models with lower AIC and BIC values were considered superior. A covariate was considered significant if its inclusion decreased the OFV by > 3.84 (χ²-test, p < 0.05, df = 1) and if backward elimination of the covariate increased the OFV by > 10.83 (χ²-test, p < 0.001, df = 1). Moreover, covariates were included only if they had a clear pharmacological or biological basis. During the model development process, condition numbers were calculated and maintained at ≤ 1,000 to avoid over-parameterization (Owen and Fiedler-Kelly, 2014).

In addition to GOF plots, model adequacy was rigorously evaluated through prediction-corrected visual predictive checks (pcVPCs), employing 2,000 Monte Carlo simulations to account for parameter uncertainty (Bergstrand et al., 2011). Statistical agreement was assessed by comparing the 95% confidence intervals (CIs) of simulated trajectories (median, 5th and 95th percentiles) against observed data distributions across automatically determined time intervals. Quantitative validation included visual inspection of percentile superimposition and evaluation of CI envelope coverage to confirm model robustness.

To evaluate parameter estimate robustness and precision, a nonparametric bootstrap analysis was conducted (Ette et al., 2003). Using Perl modules, 500 resampled datasets were generated through random sampling with replacement (Ette, 1997). Empirical 95% CIs and median values of the bootstrap-derived parameters with successful convergence were compared with the final model parameter estimates.

2.4 Virtual clinical trial of dosing strategies

Busulfan exposure is associated with both survival and toxicity in HCT recipients. Therefore, optimizing the target for busulfan cumulative exposure following all doses (cAUC) of 78–101 mg h/L during myeloablative conditioning can have a significant effect on survival chances (Bartelink et al., 2016; Lawson et al., 2021). Furthermore, to reduce the risk of sinusoidal obstructive syndrome (SOS), the maximum busulfan concentration (C_max) should be <1.88 ng/mL (Lawson et al., 2021; Philippe et al., 2019). Therefore, Monte Carlo simulations were performed using parameter estimates from the established semi-mechanistic model, while the probabilities of target attainment for different dosing strategies were compared. Individuals involved in the evaluation dataset were regarded as a virtual population in this simulation clinical trial.

First, prediction-based metrics (median prediction error [MDPE], median absolute prediction error [MAPE], and percentage of |PE|% within 20% [F₂₀] and 30% [F₃₀]) were calculated to assess the model predictability (Mao et al., 2018). Second, the time-concentration profiles were simulated 200 times for each virtual individual. Busulfan doses were subsequently administered as a 2-h infusion every 6 h for 4 days (total: 16 doses). For WT-based dosing strategy, patients weighing <9 kg, 9–16 kg, 16–23 kg, 23–34 kg, and >34 kg received busulfan doses of 1 mg/kg, 1.2 mg/kg, 1.1 mg/kg, 0.95 mg/kg, and 0.8 mg/kg, respectively; however, for age-based dosing strategy, patients aged <4 years received 1 mg/kg, and those aged ≥4 years received 0.8 mg/kg (Gurlek Gokcebay et al., 2015; Nguyen et al., 2004). The cAUC was calculated using numerical integration (Bartelink et al., 2016), while the probabilities of target attainment for the two dosing strategies were compared. Finally, the busulfan dose was simulated at 0.8–1.2 mg/kg, with a step of 0.05 mg/kg for each virtual individual; subsequently, a model-based optimal dosing regimen was recommended for each involved individual.

3 Results

3.1 Patients

The demographic characteristics and clinical data of the study population are presented in Table 1. In total, 636 busulfan whole-blood samples were obtained from 65 HCT recipients. Notably, all participants were randomly divided into two groups, and 100 samples from 10 HCT recipients were used for model evaluation. Concentrations below the lower quantification limit were not included in the analysis. The median of patient postnatal age was 1.5 years (range, 0.2–14.1), with 21 patients aged <1 year. A correlation chart of patient characteristics is presented in Supplementary Figure S1.

3.2 Semi-mechanistic population pharmacokinetic model development

A two-compartment model with first-order absorption was selected as the base model to describe busulfan PK. The model combining proportional and additive models provided the best results for the RUV. The BSV of the mean CL/F in the base model was 60.8% with a relative standard error of 8.0%. The parameter estimates and associated precisions are listed in Supplementary Table S1.

Mechanistic plausibility was considered as a potential covariate, and incorporated into the base model. First, four different body size metric-based allometric candidate models were tested to determine the most suitable body size descriptor. As shown in Supplementary Table S1, the influence of patient body size on busulfan disposition was best described by allometric scaling based on NFM, with the AIC reduced by −228.1. Second, eight different models based on the NFM were compared (Supplementary Table S2), as proposed by Du et al. (2022). The AIC value of Model Ⅲ, which included NFM allometrically and busulfan metabolism maturation upon CL based on PMA physiologically, was the lowest. Therefore, Model Ⅲ was selected as the basic structural model for further analysis. Furthermore, the stepwise approach was used to screen potential covariates (Supplementary Table S3), and concomitant with fludarabine was incorporated with OFV reduced by −10.8 (p < 0.001).

The AIC value was not reduced; however, the GSH compartment was added to describe the relative amount of GSH available at any time, based on theoretical mechanisms. The S_GSH was fixed at 0.026 h/mg following the study by Langenhorst et al. (2020), which indicates a net relevant GSH reduction of 0.26% per hour for each milligram of busulfan metabolism scaled to a 1L central volume of distribution. The influence of GST activity on S_GSH was also determined exponentially. The S_GSH was increased by 40.6% as GST enzyme activity increased from 0.9 nmol/min/mL to 20.7 nmol/min/mL (Figure 3). The IIV of the GST slope on the S_GSH was 103.9%, possibly owing to the small sample size.

Figure 3

Figure 3. Fraction of the factor S_GSH versus glutathione S-transferase enzyme activity. S_GSH is used to scale the association of busulfan metabolism with relevant glutathione depletion.

The additional estimation of IOV for CL significantly improved model predictions (ΔOFV −43.3, p < 0.001), when considering four distinct sampling occasions. During the backward process, when concomitant with fludarabine was removed from the model, the OFV increased by 5.4, which was <10.83. Therefore, this factor was excluded from the final model.

In the final model, all retained covariates significantly increased the OFV upon removal. Therefore, this model was accepted as the definitive final model. The final model parameter estimates and associated precisions are presented in Table 2. The condition number of the final model was 396.3. Shrinkage analysis for CL revealed a mean η-CL shrinkage and ε-shrinkage of 3.3% and 21.7%, respectively, which was accepted as it was below the critical threshold of 30% (Savic and Karlsson, 2009).

Table 2

Table 2. Parameter estimates for the final model and the bootstrap procedure.

3.3 Model evaluation

Ten patients in the evaluation group were included in the analysis to examine the predictability of the final model. The GOF plots for the final models, as presented in Supplementary Figure S2, show no apparent bias, with over 99.0% of observations falling within the four conditional weighted residuals. The pcVPC results showed good predictability of drug concentrations, as presented in Figure 4. The simulated data closely aligned with the observed data, indicating a lack of significant model misspecifications. The final model parameters were within the 95% confidence intervals of the bootstrap estimates, confirming the model’s stability (Table 2).

Figure 4

Figure 4. Prediction-corrected visual predictive checks (pcVPCs) for the final semi-mechanistic model. The median observed values per bin (red solid line), the 5^th and 95^th percentiles (red dashed lines) of the observations (blue circles) with the 95% confidence interval of the 5^th and 95^th percentiles (blue areas), and the confidence interval of the median (red area), are shown.

3.4 Virtual clinical trial of dosing strategies

The predicted time course of busulfan concentrations in the ten individuals involved in the evaluation dataset, which was simulated based on 200 hypothetical individuals, is presented in Figure 5. All observed concentrations were within the 5th and 95th percentiles of the simulation data, showing no trends or biases. The MDPE, MAPE, F₂₀, and F₃₀ values were −0.44%, 16.6%, 60.0%, and 74.0%, respectively. The relatively low MDPE and MAPE values further confirmed the high prediction accuracy of the final model.

Figure 5

Figure 5. The predicted time course of busulfan concentrations in the ten individuals involved in the evaluation dataset and typical patient. The 5^th - 95^th percentiles (deep pink area) and outside 5^th - 95^th percentiles (light pink area), the median (red solid line) of the simulated data, and the observations (blue circles) are shown.

The results of the probability of target attainment for different dosing strategies based on Monte Carlo simulations are presented in Table 3. Notably, three virtual patients have 8.0%–10.5% probability of exceeding the C_max safety threshold of 1.88 ng/mL when receiving a WT-based dosing strategy. In contrast, this safety risk was eliminated under age-based dosing. The WT of these three patients was within 9–16 kg, indicating the risk of SOS if the patients in this group received a WT-based dosing strategy. Compared with the WT-based dosing strategy, the dosage was relatively low in patients receiving the age-based dosing strategy. No probability of concentrations >1.88 ng/mL was observed according to the Monte Carlo simulation results; however, the infant dosage was relatively high compared with the recommended optimal dosing regimen. Therefore, GST-based dosing strategy targeting cumulative cAUC was suggested, and optimal dosing regimen was recommended for typical patients (Table 3; Figure 5).

Table 3

Table 3. The probability of target attainment of different dosing strategies based on Monte Carlo simulation.

4 Discussion

Over 40 busulfan popPK studies have been reported, of which 68% were developed predominantly in children, from which 69% and 26% of models were developed using first-order elimination and time-varying CL, respectively (Takahashi et al., 2023). Thus far, only one study has been based on a semi-mechanistic enzyme depletion model (Langenhorst et al., 2020). Furthermore, the effects of GST activity remain uncharacterized. Therefore, in this prospective study, factors affecting the PK variability of IV busulfan in pediatric HCT recipients were identified, and the effect of GST activity on the time-varying CL of busulfan was investigated using a semi-mechanistic popPK model.

In the final model, busulfan disposition was best described by including NFM allometrically as a body size metric and age-dependent maturation function. The influence of GST activity on busulfan metabolism was added exponentially based on theoretical mechanisms. The IOV was observed between dosing occasions. The developed model described busulfan PK IIV well, and model-based target attainment of different dosing strategies was evaluated.

Body size scalers and age factors were the most commonly identified covariates impacting busulfan PK in pediatric patients (Du et al., 2022; Takahashi et al., 2023). Given the broad age range of our cohort (0.2–14.1 years), we implemented NFM, a theory-based size descriptor that divides WT into FFM and fat mass (Gonzalez-Sales et al., 2022), to quantify the effect of body size and composition on busulfan PK, consistent with the study by McCune et al. (2014), Van Hoogdalem et al. (2020). In the final model, the fat mass fraction was 90.5% for CL and 68.7% for V, which was reported as 50.9% for CL, 20.3% for V by and 69.2% for CL by Du et al. (2022), McCune et al. (2014). This inconsistency in fat mass may be caused by differences in age distributions between the studies. The mean ages in the study by McCune et al. and Du et al. were 9.8 years (0.1–65.8 years) and 6.1 years (0.6–17 years), respectively. In contrast, our population averaged 2.9 years. The estimated typical CL and V_c standardized to 70-kg adult patient were 9.57 L/h and 28.2 L, respectively, which is consistent with previous studies (Takahashi et al., 2023).

To characterize developmental pharmacology, age-dependent physiological maturation functions were applied to describe changes in CL with age. Consistent with established principles, CL maturation begins before birth, making PMA more physiologically relevant than postnatal age (Anderson and Holford, 2008). According to our results, the maturation of busulfan CL reached 50% of adult values at 45 weeks PMA, which is comparable with McCune et al.’s study (TM50 = 45.7) (McCune et al., 2014). Savic et al. reported that CL increases approximately 1.7-fold between 6 weeks and 2 years by adding a nonlinear function of CL versus age to describe CL maturation (Savic et al., 2013). McCune et al. also reported that size-standardized CL reaches 95% of adult values at 2.5 postnatal years (McCune et al., 2014). Although no information on adults was available in our study, the same tendency was observed (Supplementary Figure S3), suggesting future studies should expand sample sizes to validate these maturation dynamics.

The formation of busulfan-GSH conjugates mechanistically depends on both GSH availability and depletion kinetics. Consistent with this, baseline GSH concentrations and GST polymorphisms are associated with busulfan CL (Almog et al., 2011; Takahashi et al., 2023). Among the five GST classes, the GSTA1 haplotype’ impact on CL exhibits age-dependency, reflecting the protein abundance of GSTAs in the liver increases after birth to reach adult levels during infancy (Strange et al., 1985; Strange et al., 1989). The maturation of liver drug enzymes may partially contribute to the early age-dependent PK of busulfan, indicating that GST enzyme activity may be more suitable than GST polymorphisms when assessing the influence of GST enzymes on busulfan disposition. To test this hypothesis, we incorporated GST activity effects into a semi-mechanistic PK model. Theoretically, increased GST activity would have faster GSH depletion, resulting in higher CL. However, further information on the complete developmental profile of GST enzyme activity is required in future studies.

Apart from the intra-individual variability induced by physiological maturation, busulfan-mediated GSH depletion during the treatment process may result in time-varying CL of busulfan (Langenhorst et al., 2020). The initial metabolism of busulfan occurs primarily through conjugation with endogenous GSH (spontaneously and through GST catalysis) (Myers et al., 2017). Consequently, the depletion of whole-blood GSH may contribute to the observed metabolism-dependent CL reduction (Almog et al., 2011), with clinical studies documenting 17% average CL decline from treatment initiation to day 3 (Schreib et al., 2023). Notably, no specific tendencies were noted in our study; however, IOV was included in the random-effects model to estimate course-to-course variability. The complex relationship between hepatic and blood GSH concentrations during dynamic changes may further contribute to IOV.

Other potential covariates, such as concomitant medications, disease type, and other pathophysiological indicators, were also investigated (Supplementary Table S3); however, they showed no significant effect on the busulfan PK process in this study. Notably, the effect of fludarabine co-administration and disease type has been controversial between studies (Almog et al., 2011; McCune et al., 2014; Takahashi et al., 2023). Furthermore, drug-drug interactions associated with GST depletion, such as N-acetylcysteine, should be used cautiously in the clinic (Palmer et al., 2016; Schreib et al., 2023).

Current busulfan dosing in pediatric HCT patients is based on the recommendations of regulatory agencies, such as the European Medicines Agency and U.S. Food and Drug Administration (Gurlek Gokcebay et al., 2015). Our virtual trial results corroborate previous clinical findings that WT-based dosing strategy has a higher risk of SOS compared with age-based dosing strategy (Gurlek Gokcebay et al., 2015). Specifically, Gurlek et al. previously reported that WT-based dosing was a predictor of SOS (Hazard ratio: 9.46, p = 0.009), with an SOS of 42% compared with the 5% for those receiving age-based dosing (Gurlek Gokcebay et al., 2015). The incidence of SOS is reported at 16% (range: 0%–34%), with risk factors including combining busulfan with cyclophosphamide, GSTA1 genotypes, age, weight <9 kg, weight-based dosing, and the use of once-daily IV busulfan (Lawson et al., 2021). Therefore, model-based precision dosing based on busulfan cAUC should be performed to optimize the dose regimen.

One potential limitation of this study was the lack of assessment of active GSH levels and relationship between plasma and liver GSH levels during GSH resynthesis. Therefore, full GSH dynamics could not be reconstructed in current analysis. Furthermore, comprehensive developmental profiling of GST activity from neonates to adults requires larger sample sizes to quantify its impact on busulfan CL. Additionally, as a single-center study, multicenter validation is required to enhance the model predictability.

In conclusion, we developed a semi-mechanistic popPK model to investigate the PK variability of IV busulfan in pediatric HCT recipients. Our findings demonstrate that physiologically-based descriptions of body composition according to allometry NFM, F_mat, and GST enzyme activity may mediate busulfan PK variability. Virtual clinical trial revealed that WT-based dosing strategy has a higher risk of SOS than age-based dosing strategy. Therefore, model-informed precision dosing targeting cumulative cAUC is essential for optimizing dosing regimens.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Ethics statement

The studies involving humans were approved by the Ethics Committee of the Children’s Hospital of Fudan University. The studies were conducted in accordance with the local legislation and institutional requirements. Written informed consent for participation in this study was provided by the participants’ legal guardians/next of kin.

Author contributions

DC: Data curation, Conceptualization, Writing – review and editing, Validation, Formal Analysis. XQ: Validation, Writing – review and editing, Formal Analysis, Data curation, Conceptualization. PW: Data curation, Software, Writing – review and editing. XYZ: Validation, Methodology, Writing – review and editing, Software. SH: Validation, Methodology, Writing – review and editing, Visualization. ZW: Investigation, Writing – review and editing, Data curation, Visualization. WJ: Formal Analysis, Writing – review and editing, Software, Methodology. LY: Writing – review and editing, Software, Validation, Data curation. XJ: Validation, Data curation, Writing – review and editing. YY: Data curation, Writing – review and editing, Investigation. JM: Writing – original draft, Formal Analysis, Visualization, Resources, Writing – review and editing, Conceptualization, Investigation, Supervision, Methodology, Funding acquisition, Software. XWZ: Project administration, Resources, Writing – review and editing, Supervision, Funding acquisition, Conceptualization.

Funding

The author(s) declare that financial support was received for the research and/or publication of this article. This work was supported in part by grants from the National Key R&D Program of China (No. 2023YFC2706301), National Natural Science Foundation of China (Nos 82141125 and 82204440), AI for Science Foundation of Fudan University (No. FudanX24AI065), National Key R&D Program of Shanghai (No. 2022YFC2705001), and Shanghai Hospital Development (No. SHDC12023109).

Acknowledgments

We would like to thank Editage (www.editage.cn) for English language editing.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Author disclaimer

The authors confirm that the Principal Investigator for this paper is Xiaowen Zhai and that she had direct clinical responsibility for patients.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphar.2025.1632588/full#supplementary-material

References

Almog, S., Kurnik, D., Shimoni, A., Loebstein, R., Hassoun, E., Gopher, A., et al. (2011). Linearity and stability of intravenous busulfan pharmacokinetics and the role of glutathione in busulfan elimination. Biol. blood marrow Transplant. J. Am. Soc. Blood Marrow Transplant. 17, 117–123. doi:10.1016/j.bbmt.2010.06.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Anderson, B. J., and Holford, N. H. (2008). Mechanism-based concepts of size and maturity in pharmacokinetics. Annu. Rev. Pharmacol. Toxicol. 48, 303–332. doi:10.1146/annurev.pharmtox.48.113006.094708

PubMed Abstract | CrossRef Full Text | Google Scholar

Anderson, B. J., and Holford, N. H. (2009). Mechanistic basis of using body size and maturation to predict clearance in humans. Drug metabolism Pharmacokinet. 24, 25–36. doi:10.2133/dmpk.24.25

PubMed Abstract | CrossRef Full Text | Google Scholar

Andersson, B. S., Thall, P. F., Valdez, B. C., Milton, D. R., Al-Atrash, G., Chen, J., et al. (2017). Fludarabine with pharmacokinetically guided IV busulfan is superior to fixed-dose delivery in pretransplant conditioning of AML/MDS patients. Bone marrow Transplant. 52, 580–587. doi:10.1038/bmt.2016.322

PubMed Abstract | CrossRef Full Text | Google Scholar

Bartelink, I. H., Lalmohamed, A., van Reij, E. M., Dvorak, C. C., Savic, R. M., Zwaveling, J., et al. (2016). Association of busulfan exposure with survival and toxicity after haemopoietic cell transplantation in children and young adults: a multicentre, retrospective cohort analysis. Lancet Haematol. 3, e526–e536. doi:10.1016/S2352-3026(16)30114-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Beal, S. L., Sheiner, L. B., Boeckmann, A., and Bauer, R. J. (1989). NONMEM user's guides. Ellicott City, MD, USA: Icon Development Solutions.

Google Scholar

Ben Hassine, K., Nava, T., Théoret, Y., Nath, C., Daali, Y., Kassir, N., et al. (2021). Precision dosing of intravenous busulfan in pediatric hematopoietic stem cell transplantation: results from a multicenter population pharmacokinetic study. CPT pharmacometrics and Syst. Pharmacol. 10, 1043–1056. doi:10.1002/psp4.12683

PubMed Abstract | CrossRef Full Text | Google Scholar

Bergstrand, M., Hooker, A. C., Wallin, J. E., and Karlsson, M. O. (2011). Prediction-corrected visual predictive checks for diagnosing nonlinear mixed-effects models. AAPS J. 13, 143–151. doi:10.1208/s12248-011-9255-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, Y., Xu, L. P., Zhang, X. H., Chen, H., Wang, F. R., Liu, K. Y., et al. (2018). Busulfan, fludarabine, and cyclophosphamide (BFC) conditioning allowed stable engraftment after haplo-identical allogeneic stem cell transplantation in children with adrenoleukodystrophy and mucopolysaccharidosis. Bone marrow Transplant. 53, 770–773. doi:10.1038/s41409-018-0175-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Darwich, A. S., Ogungbenro, K., Vinks, A. A., Powell, J. R., Reny, J. L., Marsousi, N., et al. (2017). Why has model-informed precision dosing not yet become common clinical reality? Lessons from the past and a roadmap for the future. Clin. Pharmacol. Ther. 101, 646–656. doi:10.1002/cpt.659

PubMed Abstract | CrossRef Full Text | Google Scholar

Donohue, M. C., Overholser, R., Xu, R., and Vaida, F. (2011). Conditional akaike information under generalized linear and proportional hazards mixed models. Biometrika 98, 685–700. doi:10.1093/biomet/asr023

PubMed Abstract | CrossRef Full Text | Google Scholar

Du, X., Huang, C., Xue, L., Jiao, Z., Zhu, M., Li, J., et al. (2022). The correlation between busulfan exposure and clinical outcomes in Chinese pediatric patients: a population pharmacokinetic study. Front. Pharmacol. 13, 905879. doi:10.3389/fphar.2022.905879

PubMed Abstract | CrossRef Full Text | Google Scholar

Ette, E. I. (1997). Stability and performance of a population pharmacokinetic model. J. Clin. Pharmacol. 37, 486–495. doi:10.1002/j.1552-4604.1997.tb04326.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Ette, E. I., Williams, P. J., Kim, Y. H., Lane, J. R., Liu, M. J., and Capparelli, E. V. (2003). Model appropriateness and population pharmacokinetic modeling. J. Clin. Pharmacol. 43, 610–623. doi:10.1177/0091270003043006007

PubMed Abstract | CrossRef Full Text | Google Scholar

Gibbs, J. P., Liacouras, C. A., Baldassano, R. N., and Slattery, J. T. (1999). Up-regulation of glutathione S-transferase activity in enterocytes of young children. Drug metabolism Dispos. Biol. fate Chem. 27, 1466–1469. doi:10.1016/s0090-9556(24)14957-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Gonzalez-Sales, M., Holford, N., Bonnefois, G., and Desrochers, J. (2022). Wide size dispersion and use of body composition and maturation improves the reliability of allometric exponent estimates. J. Pharmacokinet. pharmacodynamics 49, 151–165. doi:10.1007/s10928-021-09788-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Gurlek Gokcebay, D., Azik, F., Ozbek, N., Isik, P., Avci, Z., Tavil, B., et al. (2015). Clinical comparison of weight- and age-based strategy of dose administration in children receiving intravenous busulfan for hematopoietic stem cell transplantation. Pediatr. Transplant. 19, 307–315. doi:10.1111/petr.12430

PubMed Abstract | CrossRef Full Text | Google Scholar

Hall, S. D., Thummel, K. E., Watkins, P. B., Lown, K. S., Benet, L. Z., Paine, M. F., et al. (1999). Symposium article: molecular and physical mechanisms of first-pass extraction. Drug metabolism Dispos. Biol. fate Chem. 27, 161–166. doi:10.1016/s0090-9556(24)15271-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Hassan, M., Oberg, G., Ehrsson, H., Ehrnebo, M., Wallin, I., Smedmyr, B., et al. (1989). Pharmacokinetic and metabolic studies of high-dose busulphan in adults. Eur. J. Clin. Pharmacol. 36, 525–530. doi:10.1007/BF00558081

PubMed Abstract | CrossRef Full Text | Google Scholar

Hines, R. N. (2008). The ontogeny of drug metabolism enzymes and implications for adverse drug events. Pharmacol. and Ther. 118, 250–267. doi:10.1016/j.pharmthera.2008.02.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, H., Liu, Q., Zhang, X., Xie, H., Liu, M., Chaphekar, N., et al. (2022). External evaluation of population pharmacokinetic models of busulfan in Chinese adult hematopoietic stem cell transplantation recipients. Front. Pharmacol. 13, 835037. doi:10.3389/fphar.2022.835037

PubMed Abstract | CrossRef Full Text | Google Scholar

Karlsson, M. O., and Sheiner, L. B. (1993). The importance of modeling interoccasion variability in population pharmacokinetic analyses. J. Pharmacokinet. Biopharm. 21, 735–750. doi:10.1007/BF01113502

PubMed Abstract | CrossRef Full Text | Google Scholar

Keizer, R. J., Karlsson, M. O., and Hooker, A. (2013). Modeling and simulation workbench for NONMEM: tutorial on pirana, PsN, and xpose. CPT pharmacometrics and Syst. Pharmacol. 2, e50. doi:10.1038/psp.2013.24

PubMed Abstract | CrossRef Full Text | Google Scholar

Langenhorst, J. B., Boss, J., van Kesteren, C., Lalmohamed, A., Kuball, J., Egberts, A. C. G., et al. (2020). A semi-mechanistic model based on glutathione depletion to describe intra-individual reduction in busulfan clearance. Br. J. Clin. Pharmacol. 86, 1499–1509. doi:10.1111/bcp.14256

PubMed Abstract | CrossRef Full Text | Google Scholar

Lawson, R., Staatz, C. E., Fraser, C. J., and Hennig, S. (2021). Review of the pharmacokinetics and pharmacodynamics of intravenous busulfan in paediatric patients. Clin. Pharmacokinet. 60, 17–51. doi:10.1007/s40262-020-00947-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Long-Boyle, J. R., Savic, R., Yan, S., Bartelink, I., Musick, L., French, D., et al. (2015). Population pharmacokinetics of busulfan in pediatric and young adult patients undergoing hematopoietic cell transplant: a model-based dosing algorithm for personalized therapy and implementation into routine clinical use. Ther. drug Monit. 37, 236–245. doi:10.1097/FTD.0000000000000131

PubMed Abstract | CrossRef Full Text | Google Scholar

Mao, J. J., Jiao, Z., Yun, H. Y., Zhao, C. Y., Chen, H. C., Qiu, X. Y., et al. (2018). External evaluation of population pharmacokinetic models for ciclosporin in adult renal transplant recipients. Br. J. Clin. Pharmacol. 84, 153–171. doi:10.1111/bcp.13431

PubMed Abstract | CrossRef Full Text | Google Scholar

McCune, J. S., Bemer, M. J., Barrett, J. S., Baker, K. S., Gamis, A. S., and Holford, N. H. G. (2014). Busulfan in infant to adult hematopoietic cell transplant recipients: a population pharmacokinetic model for initial and Bayesian dose personalization. Clin. cancer Res. official J. Am. Assoc. Cancer Res. 20, 754–763. doi:10.1158/1078-0432.CCR-13-1960

PubMed Abstract | CrossRef Full Text | Google Scholar

Miyagi, S. J., Brown, I. W., Chock, J. M., and Collier, A. C. (2009). Developmental changes in hepatic antioxidant capacity are age-and sex-dependent. J. Pharmacol. Sci. 111, 440–445. doi:10.1254/jphs.09223sc

PubMed Abstract | CrossRef Full Text | Google Scholar

Myers, A. L., Kawedia, J. D., Champlin, R. E., Kramer, M. A., Nieto, Y., Ghose, R., et al. (2017). Clarifying busulfan metabolism and drug interactions to support new therapeutic drug monitoring strategies: a comprehensive review. Expert Opin. drug metabolism and Toxicol. 13, 901–923. doi:10.1080/17425255.2017.1360277

PubMed Abstract | CrossRef Full Text | Google Scholar

Nguyen, L., Fuller, D., Lennon, S., Leger, F., and Puozzo, C. (2004). I.V. busulfan in pediatrics: a novel dosing to improve safety/efficacy for hematopoietic progenitor cell transplantation recipients. Bone marrow Transplant. 33, 979–987. doi:10.1038/sj.bmt.1704446

PubMed Abstract | CrossRef Full Text | Google Scholar

Owen, J. s., and Fiedler-Kelly, J. (2014). Introduction to population pharmacokinetic/pharmacodynamic analysis with nonlinear mixed effects models.

Google Scholar

Palmer, J., McCune, J. S., Perales, M. A., Marks, D., Bubalo, J., Mohty, M., et al. (2016). Personalizing busulfan-based conditioning: considerations from the American society for blood and marrow transplantation practice guidelines committee. Biol. blood marrow Transplant. J. Am. Soc. Blood Marrow Transplant. 22, 1915–1925. doi:10.1016/j.bbmt.2016.07.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Philippe, M., Neely, M., Rushing, T., Bertrand, Y., Bleyzac, N., and Goutelle, S. (2019). Maximal concentration of intravenous busulfan as a determinant of veno-occlusive disease: a pharmacokinetic-pharmacodynamic analysis in 293 hematopoietic stem cell transplanted children. Bone marrow Transplant. 54, 448–457. doi:10.1038/s41409-018-0281-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Savic, R. M., and Karlsson, M. O. (2009). Importance of shrinkage in empirical bayes estimates for diagnostics: problems and solutions. AAPS J. 11, 558–569. doi:10.1208/s12248-009-9133-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Savic, R. M., Cowan, M. J., Dvorak, C. C., Pai, S. Y., Pereira, L., Bartelink, I. H., et al. (2013). Effect of weight and maturation on busulfan clearance in infants and small children undergoing hematopoietic cell transplantation. Biol. blood marrow Transplant. J. Am. Soc. Blood Marrow Transplant. 19, 1608–1614. doi:10.1016/j.bbmt.2013.08.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Schreib, K. M., Bram, D. S., Zeilhofer, U. B., Muller, D., Gungor, T., Kramer, S. D., et al. (2023). Population pharmacokinetic modeling for twice-daily intravenous busulfan in a large cohort of pediatric patients undergoing hematopoietic stem cell Transplantation-A 10-Year single-center experience. Pharmaceutics 16, 0. doi:10.3390/pharmaceutics16010013

PubMed Abstract | CrossRef Full Text | Google Scholar

Scian, M., Guttman, M., Bouldin, S. D., Outten, C. E., and Atkins, W. M. (2016). The myeloablative drug busulfan converts cysteine to dehydroalanine and lanthionine in redoxins. Biochemistry 55, 4720–4730. doi:10.1021/acs.biochem.6b00622

PubMed Abstract | CrossRef Full Text | Google Scholar

Shukla, P., Goswami, S., Keizer, R. J., Winger, B. A., Kharbanda, S., Dvorak, C. C., et al. (2020). Assessment of a model-informed precision dosing platform use in routine clinical care for personalized busulfan therapy in the pediatric hematopoietic cell transplantation (HCT) population. Front. Pharmacol. 11, 888. doi:10.3389/fphar.2020.00888

PubMed Abstract | CrossRef Full Text | Google Scholar

Strange, R. C., Davis, B. A., Faulder, C. G., Cotton, W., Bain, A. D., Hopkinson, D. A., et al. (1985). The human glutathione S-transferases: developmental aspects of the GST1, GST2, and GST3 loci. Biochem. Genet. 23, 1011–1028. doi:10.1007/BF00499944

PubMed Abstract | CrossRef Full Text | Google Scholar

Strange, R. C., Howie, A. F., Hume, R., Matharoo, B., Bell, J., Hiley, C., et al. (1989). The development expression of alpha-mu- and pi-class glutathione S-transferases in human liver. Biochimica biophysica acta 993, 186–190. doi:10.1016/0304-4165(89)90162-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Takahashi, T., Jaber, M. M., Brown, S. J., and Al-Kofahi, M. (2023). Population pharmacokinetic model of intravenous busulfan in hematopoietic cell transplantation: systematic review and comparative simulations. Clin. Pharmacokinet. 62, 955–968. doi:10.1007/s40262-023-01275-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Ten Brink, M. H., Swen, J. J., Bohringer, S., Wessels, J. A., van der Straaten, T., Marijt, E. W., et al. (2013). Exploratory analysis of 1936 SNPs in ADME genes for association with busulfan clearance in adult hematopoietic stem cell recipients. Pharmacogenetics genomics 23, 675–683. doi:10.1097/FPC.0000000000000007

PubMed Abstract | CrossRef Full Text | Google Scholar

Trame, M. N., Bergstrand, M., Karlsson, M. O., Boos, J., and Hempel, G. (2011). Population pharmacokinetics of busulfan in children: increased evidence for body surface area and allometric body weight dosing of busulfan in children. Clin. cancer Res. official J. Am. Assoc. Cancer Res. 17, 6867–6877. doi:10.1158/1078-0432.CCR-11-0074

PubMed Abstract | CrossRef Full Text | Google Scholar

van Hoogdalem, M. W., Emoto, C., Fukuda, T., Mizuno, T., Mehta, P. A., and Vinks, A. A. (2020). Population pharmacokinetic modelling of busulfan and the influence of body composition in paediatric fanconi anaemia patients. Br. J. Clin. Pharmacol. 86, 933–943. doi:10.1111/bcp.14202

PubMed Abstract | CrossRef Full Text | Google Scholar

West, G. B., Brown, J. H., and Enquist, B. J. (1997). A general model for the origin of allometric scaling laws in biology. Science 276, 122–126. doi:10.1126/science.276.5309.122

PubMed Abstract | CrossRef Full Text | Google Scholar