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.

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.

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 m a t =     1     1 + P M A T M 50 − Hill (1) where TM50 is the PMA at which maturation achieving 50% of the adult value, and Hill defines the steepness of the sigmoid decline.

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).

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). d A G S H d t = S G S H V 1 × A G S H × k 10 × A b u 1 (2) Where A b u 1 represents the amount of busulfan in the central compartment, A G S H 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.

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.

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.

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.

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).

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).

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.

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).