Optimization of Busulfan Dosing Regimen in Pediatric Patients Using a Population Pharmacokinetic Model Incorporating GST Mutations

Introduction

Busulfan (Bu) is a common alkylating agent, which can bind to the guanine of intracellular DNA to damage its structure and function.¹ Bu-based conditioning schemes are regarded as the cornerstone of allogeneic hematopoietic stem cell transplantation (allo-HSCT) and autologous transplantation because of their myeloablative activity.^1,2 However, the IV treatment window of IV Bu recommended by FDA³ is narrow (900–1350 μM·min). Higher exposure (expressed as area-under-the-curve of 0–6 h, AUC_0-6h) (>1350 μM·min) may lead to a higher probability of sinusoidal obstructive syndrome (SOS) and acute graft-vs-host disease (aGVHD), while low AUC_0-6h (<900 μM·min) of Bu might increase the risk of graft failure.^4–6 Therapeutic drug monitoring (TDM) for busulfan based on the first-dose pharmacokinetics (PK) has been confirmed to decrease the occurrence of toxicity and to improve efficacy.^7,8

Bu is metabolized by the formation of a glutathione conjugate in the liver.^9,10 This reaction is primarily catalyzed by glutathione S-transferase (GST) enzymes, such as GSTA1, GSTM1 and GSTP1.¹¹ GSTA1 is the predominant GST enzyme involved in Bu metabolism, and the activity of GSTM1 is close to the half of GSTA1.¹² However, the activity of GSTP1 is relatively weak.^2,12 Hence, polymorphisms in the GSTA1 or GSTM1 isoenzymes would most likely affect Bu metabolism. It has been reported that nucleotide polymorphisms (SNPs) for GSTA1, which are linked and located within the proximal promoter of GSTA1 gene, with the *B haplotype resulting in the lower gene transcription.¹³

Studies about the relationships between GSTA1 mutations and Bu PK have yielded inconsistent results in adults. For example, Abbasi et al¹⁴ stated that GSTA1 was not significant influence factor of intravenous (IV) Bu clearance. However, Michaud et al¹⁵ reported that adult patients carrying GSTA1 loss of function *B allele were at increased risk of overdosing on their initial busulfan oral dose. Yin et al¹⁶ stated adult patients with the GSTA1 *A/*B genotype had a significantly higher AUC, higher peak concentration (C_max) and lower clearance (CL).

Relationships between GSTA1 mutations and Bu PK in children has not been elucidated as well. According to the basic research, the expression of GST enzymes in children has been reported to be higher than that in adults.¹⁷ Furthermore, the activity of GST enzymes is also affected by race. Minor allele frequencies (MAF) taken from 1000 Genomes were reported to differ between Caucasian and Asian populations (https://www.ncbi.nlm.nih.gov/snp/rs3957356#frequency_tab). Until now, there has been only one population PK model for busulfan that successfully incorporated GSTA1 genotype in a Caucasian pediatric population.¹⁸ Little study has quantificationally analyzed GST mutations potentially affecting Bu PK in an Asian pediatric population.

In addition to GSTA1 genotypes, the influence of other factors on Bu PK could not be ignored. PPK models developed by Trame et al¹⁹ and Booth et al²⁰ emphasized the non-negligible role of BSA and body weight (BW) on individualizing busulfan therapy in children. It has been also suggested that the effect of age,²¹ primary disease, hepatic function and drug interactions (fludarabine and phenytoin) may partly explain interindividual variability on Bu PK.^22–24

The aim of this study was thus to build a population pharmacokinetic (PPK) model with data from pediatric patients, so as to present the PK feature of Bu, to reveal the variability of PK parameters, and to identify the potential contribution of covariates on the disposition of Bu. Then, we tried to develop a novel busulfan dosing regimen, based on the PPK model, to achieve better AUC_0-6h targeting. Furthermore, maximum a posteriori (MAP) Bayesian forecasting made use of a PPK model and limited number of samples to forecast AUC_0-6h and to formulate an optimal LSS, which is an alternative monitoring strategy.

Patients and Methods

Patients and Treatment Regimens

From March 2019 to April 2020, we collected 76 patients received allo-HSCT. This study followed the tenets of the Declaration of Helsinki promulgated in 1964. The research was approved by the Beijing’s Children Hospital and all patients/parents provided informed consent.

In the pretreatment regimen, intravenous infusion of busulfan (Busulfex, Otsuka Pharmaceutical Co. Ltd, Zhejiang, China) was once every six hours for three or four days, for a total of 12 or 16 doses, and infusion for each dose lasted for two hours. The first dose started seven-to-nine days before HSCT (ie day −7, −8, or −9). The dosage of Bu was based on weight (1 mg/kg for children less than 9 kg, 1.2 mg/kg for children within a 9–16 kg range, 1.1 mg/kg for children between 16 and 23 kg, 0.95 mg/kg for children between 23 and 34 kg, and 0.8 mg/kg for children more than 34 kg) in accordance with the European Medicines Agency (EMA).²⁵ Different conditioning regimens were used dependent on the primary diseases and the types of donor. Briefly, Bu combined with cyclophosphamide (CTX) were the basic components of conditioning chemotherapy. For patients with malignant diseases, regimens containing cytarabine (Ara-C) were commonly applied, while fludarabine (FLU) was administrated for myeloablative treatment of patients with nonmalignant tumors. Specific dosage and duration of different conditioning regimens are shown in Supporting Information Table 1. Phenytoin (PHT) started 30 min before the initiation of Bu therapy in order to prevent central nervous system toxicity caused by Bu. Moreover, cyclosporine (CsA) with or without methotrexate (MTX)/mycophenolate mofetil (MMF) was given as GVHD prophylaxis.

Bu Determination and Genotyping

Blood samples for PK analysis and genotyping were withdrawn from central venous lines, in heparinized glass tubes, pre-infusion, 0.5, 1, 2, 2.5, 4 and 6 hours after the first infusion. Plasma concentrations of Bu were determined using a high performance liquid chromatography-tandem mass spectrometry (HPLC-MS/MS).¹ The lower limit of quantitation was 10 ng/mL and the range of quantitation was from 10 to 10,000 ng/mL.

Pretransplant genomic DNA was isolated and extracted from whole blood prior to the first Bu infusion. Improve potassium iodide methods was applied for DNA extraction from whole blood. Ammonium chloride was used to destroy red blood cells, and potassium iodide was used to destroy white blood cells and their nuclear membranes for a short time. Then, proteins, lipids and residual cell debris are precipitated by chloroform/isopropanol. Finally, DNA is precipitated by isopropanol, and was washed by ethanol. The extracted genomic DNA was dissolved with TE. After confirming the DNA concentration and purity, PCR amplification and purification of PCR product were performed. GST genotypes of patients, GSTA1 (rs3957356 and rs3957357, which defines haplotype *A and *B) and GSTM1 (rs3754446), were detected with the ABI 3730XL DNA Analyzer (Applied Biosystem).²⁶ Supporting Information Table 2 displays the primer sets and Tm used for the genotyping assays.

PPK Analysis

BU plasma concentration-time data was analyzed using a nonlinear mixed-effects (NLME) model implemented in Phoenix 8.0 (Certara USA Inc., Princeton, NJ, USA). Typical PPK parameters and their random interindividual variability (IIV) were estimated using a first-order conditional estimation method with extended least squares method (FOCE ELS). One- or two-compartment distribution with first-order elimination were tested as structural model. Exponential errors follows a log-normal distribution, assumed to describe the IIV in PK parameters by the equation P_i=P× , where P_i is the individual PK parameter of the ith individual, P is the geometric average population value, and η_i is the subject-specific random effect value, normally distributed random variable with a mean of 0 and a variance of ω^2.27 Additive, proportional, combined additive, and proportional models were evaluated to account for the intraindividual variability, which were shown by Equation 1~3:

(1)

(2)

(3)

Where and represented the observed and predicted ith concentration in the jth patient respectively. was additive error, and was proportional error. The selection of base model was based on changes in −2 log likelihood (−2LL) and on graphic analyses of the goodness-of-fit (GOF).

Based on scatter-plots of individual parameters estimated by base model vs covariates and clinical plausibility, we selected candidate covariates for each PK parameters.²⁸ The influence of potential covariates ([sex, age, body weight, BSA, alanine transaminase (ALT), aspartate aminotransferase, AST; alkaline phosphatase, ALP; total bilirubin, TBIL; creatinine, CRE; creatinine clearance, CLcr; GSTA1 genotypes, *A/*A, *A/*B, *B/*B; GSTM1 genotypes; primary disease, malignant disease, nonmalignant disease; and drug interactions, FLU) on clearance (CL) and volume (V) were further investigated using a stepwise procedure. Categorical covariates were coded as number, and continuous covariates were centered on their median value.²⁷ During forward selection, covariates were defined as significance if the −2LL decreased by at least 3.84 (p≤0.05) following their inclusion in the model. During backward elimination, one covariate could remain in the final model if the −2LL increased by at least 6.63 (p≤0.01) when removed at a time from the full model.

Model Validation

Shrinkage in individual random effects were evaluated in order to assess whether the final model could be capable to estimate individual PK parameters by taking advantage of population typical values and sparse PK data. Shrinkage values of less than 20% indicate that the individual data are rich enough to compute the PK parameters, whereas larger shrinkage values generally mean that individual Bayesian estimates are biased towards the population mean values.²⁹

Graphical observation of the final model adopted GOF, including conditional weighted residuals (CWRES) over population predicted concentrations (PRED) or time after dose (TAD) and the relationship between observed (OBS) and PRED or individual predicted concentrations (IPRED).³⁰ In addition, the final model was evaluated using visual predictive check (VPC) and bootstrap analysis. VPC was based on the final pharmacokinetic estimates and then calculated the 95%CI for concentrations by simulating 1000 individuals. And, VPC compared 5th, 50th, 95th percentiles of the simulated concentrations with those of the observed concentrations. In the bootstrap analysis,³¹ the 95%CI of the parameter estimates were derived from 1000 datasets of 69 subjects generated by random sampling using the Phoenix NLME.

External validation of the model was performed using an external dataset to evaluate the predictive performance of the final PPK model. The external dataset consisted of 81 busulfan concentrations from 14 children undergoing allo-HSCT. Busulfan plasma concentrations were predicted by fixing the population PK parameters to the final estimates of the previously established model and setting maximum evaluations to 0. From this study, measured concentrations of individual patients assigned to busulfan were compared with calculated concentrations of these individual patients at the same time with our PPK model using their BSA, ALT and GSTA1 genotypes of patients. Differences of <20% between calculated and measured concentrations were allowed.³²

Dosing Regimen Optimization

A new dosage scheme using a simulated dataset was designed to achieve a targeted AUC_0-6h of 1125 μM·min.²⁵ The distributions of BSA and AST level in the simulated dataset were tried to be coincident with the clinical dataset. The BSA of actual patients were distributed from 0.2 m² to 1.6 m², which indicated that BSA of children almost fall within this range. So, the simulated BSA values were set as 0.2, 0.3, 0.4, …, 1.5, and 1.6. As to the simulation of AST, the first consideration is clinical significance. We assumed that every simulated patient (every BSA value) has the possibility of different AST levels, and the setting of simulated AST referred to the opinion of the hematologists. Based on this assumption that the ranges of BSA and AST level in simulated patients were similar between different genotypes, the number of patients carrying GSTA1 *A*A and *A*B was equivalent. A total of 480 simulated patients were included in simulation dataset. The Dose_simulated was designed according to EMA regimen and the simulated dosing frequency was once every six hours. AUC_0-6h was the integral of concentration over time (0–6 h). One thousand simulations were conducted in per patient to take into account intraindividual variability. Then, 1000 simulated AUC_0-6h values per patient and mean value of each patient were generated in Phoenix-NLME. Model-based doses were recommended according to the following Equation 4:

(4)

Then, AUC_{0-6h, recommended} values were simulated according to the new recommended doses (Dose_recommended). Variabilities of AUC_{0-6h, recommended} were shown by the range between minimum and maximum values and the coefficient of variation (CV). The rate of achieving the targeted AUC window (900–1350 μM·min) was also listed. Finally, we plotted the relationship between recommended dose and BSA, and derived the regression equations. Meanwhile, in order to facilitate clinical application, we present a piecewise dosing nomogram based on BSA categories for each genotype.

In order to verify the application of new regimen, we explored how it worked for the real patients whose data was used in the development of PPK model (n=69) and its external validation (n=14). We calculated the recommended doses for the real patients based on regression equations. And, the relevant predicted AUC_0-6h was calculated as follows: AUC_predicted=(AUC_observed/Dose_{administrated})×Dose_designed. Then, we checked whether the expected AUC_0-6h fell within the therapeutic window (900–1350 µM∙min).

Moreover, we performed a comparison with other five analogous pediatric PPK models for busulfan. Doses were designed based on the individual Bu CL obtained from each publication^{18,20,21,33,34} and target AUC_0-6h (1125 µM∙min). The detail dosing regimens have been demonstrated in Supporting Information Table 3. The predicted AUC_0-6h was calculated as mentioned above and was checked.

Optimized Sampling Scheme

The final PPK model provides information on typical PK parameters and variabilities associated with these values as well as how specific covariates (eg BSA, AST, genotypes) influence the PK of Bu. After evaluating final model predictive performance, optimal LSS was considered to optimize the sampling scheme for TDM. Candidate sampling points included: (1) 2, 2.25, 4, and 6 h; (2) 2, 2.5, 4, and 6 h; (3) 2.25, 2.5, 4, and 6 h; (4) 2, 4, and 6 h; (5) 2.25, 4, and 6 h; (6) 2.5, 4, and 6 h; (7) 2, 2.25, and 4 h; (8) 2, 2.25, and 6 h; (9) 2, 2.5, and 4 h; (10) 2, 2.5, and 6 h; (11) 2.25, 2.5, and 4 h; (12) 2.25, 2.5, and 6 h; (13) 2, and 4 h; (14) 2, and 6 h; (15) 2.25, and 4 h; (16) 2.25, and 6 h; (17) 2.5, and 4 h; (18) 2.5, and 6 h after the beginning of the first infusion. The plasma drug concentration at 2.25 hs after administration was simulated by the final model. Each candidate LSS was tested by Bayesian method to obtain AUC_0-6h of each patient. APE, MAPE and rRMSE were calculated by comparing the predicted AUC_0-6h derived from the LSS with the actual AUC_0-6h obtained from each patient’s full concentration-time samples, as according to Equation 5~7 to evaluate predictive accuracy:

(5)

(6)

(7)

The lower the MAPE value, the better the LSS. A sampling strategy was considered to display a good predictive performance when the 95%CI around the MAPE was less than 20% of the reference Bu AUC_0-6h values.³⁵ Bland–Altman plots were applied to assess the agreement between the predicted and the actual AUC_0-6h, and the difference was expressed as a mean ±1.96 SD.

Definition of Clinical Outcomes

The definition of engraftment was absolute neutrophil count (ANC) ≥0.5×10⁹ per liter for the first three days in a row and platelet count (PLT) ≥20×10⁹ per liter for the first one week after HSCT.^8,16 Primary engraftment failure was defined as failing to reach an ANC of 0.5×10⁹/L within 30 days after HSCT.¹⁶ In brief, the diagnosis of SOS was according to the modified Seattle criteria.³⁶ The Mount Sinai Acute GVHD International Consortium (MAGIC) was taken as the diagnostic criterion for aGVHD.³⁷

Statistical Analysis

Associations with pharmacokinetics parameters (CL and AUC_0-6h) and GSTA1 mutations were tested with t-test,³⁸ whereas time-to-event data was compared by the log rank test. All statistical analyses were performed with SPSS version 21.0 (IBM Inc., IL, USA). In this study, P-values <0.05 are considered as statistical significance, and P-values <0.1 were deemed as trends.

Results

Population Characteristics

From March 2019 to April 2020, 76 patients underwent allo-HSCT and received IV busulfan. The GSTA1 genotypes for four patients and GSTM1 genotypes for six patients were missing because of samples being unavailable or detection failure. GSTA1-52C and T defined haplotype *A and *B, respectively. In relation to GSTA1-52T/C and −69A/G SNPs, a higher percentage of CC/GG (78.08%) vs TT/AA (1.37%) and CT/GA (20.55%) was found in our model-building patient cohort. There was only one patient with GSTA1 *B*B. Thus, 69 patients were included in our final PPK analysis. Table 1 summarizes characteristics of patients including demographics, primary diseases, donor types and conditioning regimens. Supporting Information Table 4 further elaborates patients’ diagnoses. The genotyping results are shown in Table 2.

Table 1 Patients’ Demographic and Clinical Characteristics of Model-building

Table 2 Genetic Characteristics of Study Population

PPK Model

A total of 398 blood samples were collected for the PPK analysis. Figure 1 plotted concentration–time profiles of Bu. A one-compartment model with first-order elimination was suitable to describe the profiles of busulfan pharmacokinetics, and the exponential model was feasible to estimate the interindividual variability. Neither proportional error model nor additive error model could perform well, while a combined proportional and additive residual error model provided an adequate fit (Δ-2LL >20, p<0.001).

Figure 1 Observed busulfan concentrations plotted vs time. The filled dark dots represent the observed concentrations of patients with GSTA1 *A/*A. Gray circles represent the observed concentrations of patients with GSTA1 *A/*B. The solid black lines and gray dashed lines represent the mean (A) and “typical subjects” (B) of the observed concentrations of patients with GSTA1 *A/*A and GSTA1 *A/*B, respectively. “Typical subjects” of two genotypes were chosen based on the PPK model profile. Abbreviation: GST, glutathione S-transferase enzyme.

In the forward model building step, BSA, AST, types of primary diseases and GSTA1 genotypes declined the value of −2LL by more than 3.84 (p<0.05) at each addition. The GSTA1 mutation, BSA and the level of AST significantly affected the clearance (CL) of busulfan. Meanwhile, the volume (V) of busulfan was influenced by BSA and types of primary diseases, while the GSTM1 genotypes had no significant effect on PK parameters.

During backward elimination, the value of −2LL increased significantly (Δ-2LL >6.63, p<0.01) when respectively removing the GSTA1 mutation, BSA and the level of AST out of the model. Thus, the three covariates mentioned above could remain in the final model. Types of primary diseases were eliminated in the final model because it failed to significantly increase the −2LL value (Δ-2LL=6.333, p>0.01). Table 3 presented the detailed PK parameters of the base and final model, which were standardized to the median BSA, median AST and GSTA1 *A/*A. The CL and V of busulfan were affected by the GSTA1 genotypes, BSA and AST as illustrated by the following Equations 8 and 9:

(8)

(9)

Table 3 Final PPK Model: Pharmacokinetic Parameter Estimates and Bootstrap Results

When patients carried the GSTA1 *A/*B mutation, GSTA1 = 1, and when patients carried the *A/*A genotype, GSTA1 = 0.

Population estimates for a patient with median BSA with GSTA1 *A/*A were CL 4.79 L/h and V 14.8 L. The values of interindividual variability for CL and V were 18.65% and 23.63%, respectively. CL of patients carrying the GSTA1 *A/*B genotype was predicted to decline by 17.3% compared with those carrying GSTA1 *A/*A. As a result, CL/BSA between two genotypes had a significant difference (p=0.0103), in spite of AUC_0-6h/BSA without significant difference (Figure 2).

Figure 2 The effect of GSTA1 genotypes on CL/BSA (A) and AUC0-6h/BSA (B). AUC0-6h, area under the concentration-time curve of 0–6 h. Abbreviations: BSA, body surface area; CL, clearance; GST, glutathione S-transferase enzyme.

Model Evaluation and Validation

Shrinkage values of clearance and volume were respectively 0.207 and 0.130, which manifested that the model-building dataset was rich enough to compute the PK parameters of busulfan in children.

Compared with the base model, the −2LL value of final model declined by 175.26, which indicated that the model improved substantially after incorporating covariates, BSA, AST and GSTA1 phenotype. As demonstrated in the scatter plots of OBS vs PRED (Figure 3A and B), the PRED strongly deviated from the OBS in the base model, but PRED agreed with the OBS in the final model. In addition, the CWRES in the final model were more uniformly distributed within the accepted range (y=±2) than that in the base model (Figure 3C and D). Supporting Information Figure 1 shown the OBS vs IPRED and CWRES vs TAD of final model. As a whole, the final model was visually improved in terms of the GOF and had a more accurate predictive performance.

Figure 3 Goodness-of-fit plots of the base model and final model. (A) OBS vs PRED of base model; (B) OBS vs PRED of final model; (C) CWRES vs PRED of base model; (D) CWRES vs PRED of final model. The black solid lines are the reference lines, and blue solid lines are the linear regression lines. Abbreviations: CWRES, conditional weighted residuals; OBS, observed concentrations; PRED, population predicted concentrations.

The estimated values generated by bootstrap analysis were close to the parameters in the final model, which also fell within 95%CIs (Table 3). Thus, it can be considered that the final model was accurate and robust. Figure 4 clearly displayed that 5th, 50th, and 95th percentiles of simulation values almost coincided with that of observed concentrations. The results of VPC combined with bootstrap analysis confirmed the exactitude of the parameter estimates and demonstrated the reliability of the final model. For external validation, a mean difference of 18.48% was observed in the simulated Bu concentrations generated using the established PPK model when compared with the observed Bu concentrations (Figure 5).

Figure 4 Visual prediction check of final model. The dots represent observed concentrations and dashed lines represent 5th, 50th, and 95th percentiles of observed concentrations; solid lines represent 5th, 50th, and 95th percentiles of simulated concentrations; shaded areas represent 95%CIof the 5th, 50th, and 95th prediction intervals. Abbreviation: TAD, time after closing.

Figure 5 Comparison between the predictive and the observed busulfan concentrations by time from external population. Black circles represent the observed Bu concentrations; red dots represent the predictive Bu concentrations. Abbreviation: TAD, time after dosing.

New Dosing Strategy and Simulations

Based on target AUC, recommended doses (Q6 h) of per BSA were obtained. There was a linear relationship between dose and BSA (Figure 6A). The linear Equation 10 and 11 for GSTA1 *A/*A and *A*B was respectively, in which the range of BSA was 0.2 m² to 1.6 m²:

(10)

(11)

Figure 6 Dosing curves: total doses in mg (A) and doses in mg/m2 (B). The solid line and dashed line respectively represent dosage for patients with GSTA1 *A/*A and GSTA1 *A/*B. Abbreviations: BSA, body surface area; GST, glutathione S-transferase enzyme.

The dosages were normalized by BSA (mg/m²) due to the fact that body size-based dosing is more familiar to pediatricians. As demonstrated in Figure 6B, the doses normalized by BSA (Dose/BSA) were shown to be a homographic function of BSA. Dosages according to this hyperbola would be inconvenient for application in clinical, so a nomogram of BSA categories was derived. Three dosages for each genotype were then defined according to selected BSA categories in order for convenient clinical application. Finally, a new dosing strategy (once every six hours) based on BSA normalization was listed in Table 4.

Table 4 New Dosing Strategy (Q6h) Based on BSA Normalization

BSA-based dosing recommended by the final PPK model was helpful for targeting the patient AUC. The success rate of achieving the targeted AUC_0-6h window (900–1350 μM·min) was 99.58% in simulated patients. The new dosage yielded homogeneous AUC_0-6h values in different BSA categories, and the CV of 7.57% in AUC was low.

Approximate 80% of the predicted AUC_0-6h of real patients fell within the therapeutic window, which verified the applicability and feasibility of the new regimen once again. As for the other five dosage schemes, their success rates were all below 50%.

Limited Sampling Strategies

Eighteen models with different sampling time points are all listed in Table 5. The relationship between the predicted and actual AUC_0-6h for these models is shown in Supporting Information Figure 2. From two to four-point models, 83.33% LSSs fitted well with the correlation (r²) of more than 0.85. Prediction precision of LSSs expressed as rRMSE and MAPE is also given in Table 5. Model 2 (C_2h, C_2.5h, C_4h and C_6h) showed not only the best fit to the Bu AUC_0-6h, but also better prediction precision (rRMSE=0.72% and MAPE=4.55%) than other LSSs. Among the three-point models, Model 9 (C_2h, C_2.5h, and C_4h) and Model 10 (C_2h, C_2.5h, and C_6h) both behaved well, and with similar prediction precision to Model 2. Within these models, no patients had an AUC_0-6h lower than −15% or higher than +15% of the reference value. As for these two-sampling schemes, Model 13 (C_2h and C_4h), Model 14 (C_2h and C_6h), Model 17 (C_2.5h and C_4h), and Model 18 (C_2.5h and C_6h) had relatively low rRMSE and MAPE. As the Bland–Altman plots shown, the means difference between actual and predicted AUC_0-6h were close to 0, and dots almost fell within the lines of mean ±1.96 SD, which indicated that the actual and predicted AUC_0-6h agreed well in the seven models above (Supporting Information Figure 3).

Table 5 LSSs for IV BU Determined by Bayesian Estimate After the First Dose

Clinical Outcomes

Among the clinical outcomes analyzed (n=69), five patients died after HSCT, one of whom died of SOS and the other four died of severe infection. Graft failure occurred in two patients with AUC_0-6h of 659.9 and 482.7 μM·min. Engraftment was achieved for 97.10% of patients (median time: 12 days, range: 10–19 days) for neutrophils within 30 days after transplantation and 59 patients achieved engraftment for platelets within a median of 15 days (range: 7–30). As shown in Figure 7, there were significant differences in ANC recovery and survival rate between patients with two GSTA1 genotypes (p<0.05). Seven patients (10.14%) developed SOS and aGVHD I–IV was documented in 57 patients (82.61%). The clinical outcomes of patients with two GSTA1 genotypes were compared in Supporting Information Table 5.

Figure 7 GSTA1 polymorphism and clinical outcomes of HSCT. ANC recovery (A) and survival rate (B) plotted against two GSTA1 genotypes. The solid line and dashed line respectively represent clinical outcomes of patients with GSTA1 *A/*A and GSTA1 *A/*B. P values are shown on the plots. Abbreviations: ANC, absolute neutrophil count; GST, glutathione S-transferase enzyme.

The logistic regression model was put into use to investigate whether Bu AUC_0-6h had an impact on the incidence of aGVHD I–IV, SOS or death, while Supporting Information Figure 4 revealed that no correlation between Bu AUC_0-6h and regimen-related toxicity or mortality was observed.

Discussion

This study developed the first PPK model for busulfan that successfully incorporated GSTA1 genotypes in a Chinese pediatric population and partly explained the source of variability of busulfan clearance, suggesting that GST genotyping would be necessary for optimization of pediatric Bu treatment.

The literature until now has been controversial regarding the correlation between Bu clearance with genetic mutations.³⁹ There have been several research studies supporting the positive association between GST mutations and Bu clearance. In a pediatric multicenter study, Ansari et al¹³ reported that the activity of GSTA1 promoter was significantly descended in the case of *B haplotype compared with *A haplotype. Another study succeeded in incorporating the GSTA1 genetic variants into a PPK model for Bu in a Caucasian pediatric population, and then tailored the dose according to the individual metabolic capacity.¹⁸ For Chinese adults, Yin et al¹⁶ concluded that patients with GSTA1 *A/*B genotype possessed a higher C_max, higher AUC and lower clearance than the group with GSTA1 *A/*A genotype. Furthermore, in a recent study of Japanese pediatric population (n=20), Nishikawa et al⁴⁰ stated the correlation between GST mutations and clearance was distinguished.

In our two cohorts (n=84), the frequency of GSTA1 *B haplotype was 11.9% (MAF=0.119). In other words, a minority of patients had the *B haplotype, of which 1.19% were homozygous (GSTA1 *B/*B). Therefore, the MAF of Chinese patients was lower than that of global population (MAF=0.306) taken from 1000 Genomes (https://www.ncbi.nlm.nih.gov/snp/rs3957356#frequency_tab). Apart from the significant discrepancy of the MAF between different populations, the distribution of haplotypes differences existed.¹³ Therefore, studies in Caucasians could not completely represent that in the Asian population. Besides, in a study conducted in patients with acute myeloid leukemia, Yee et al⁴¹ found that C allele in the GSTM1 locus (rs3754446) was associated with decreased Bu AUC of first dose and lower disease-free survival. Thus, GSTA1 and GSTM1 genotypes were taken into consideration in the present study when exploring the influence of GST mutations on Bu PK in Chinese children.

A population pharmacokinetic model of Bu was developed to test the influence of gene mutation on the pharmacokinetic characteristics. In the final model, the estimate value (per BSA of 0.67 m²) of CL was 4.79 L/h and of V was 14.8 L. Patients carrying the GSTA1 *A/*B genotype had a 17.3% lower clearance than those carrying GSTA1 *A/*A, which is consistent with previous studies reporting that the presence of the mutation allele probably resulted in the decreased activity of GST enzyme.¹² It can be seen that GSTA1 mutations were associated with Bu CL either in Chinese children or in Caucasian children. According to the Stevenson formula [BSA (m²)=0.0061×height (cm)+0.0128×weight (kg)–0.1529], the BSA of normal adult of 70 kg and 180 cm was calculated as 1.8411 m². The results of CL in the present study were normalized to 11.08 L/h (per BW of 70 kg) for GSTA1*A/*A. Nava et al¹⁸ reported that CL in Caucasian children was 14.66 L/h for rapid metabolizers, which was a little faster than that in Chinese children. However, genetic variation in GSTM1 showed no significant impact on Bu CL, likely because the function of the GSTM1 enzyme involved in Bu metabolism was less than the GSTA1 enzyme. In the further research, merging data from multiple studies may be helpful to more precisely estimate the busulfan PK.

As shown in Figure 1, the distribution phase (0–2 hs after dosing) shown large heterogeneity between two genotypes, which may slightly interfere the statistical analysis of AUC_0-6h. Nevertheless, in the eliminate phase, it was obvious that patients with GSTA1 *A/*A appeared higher CL than those with GSTA1 *A/*B. Thus, as Figure 2 reflected, there was a significant difference in clearance per body surface area between the two GSTA1 genotypes, but there was no significant difference in AUC_0-6h per body surface area. Theoretically, AUC depends on the dose as well as CL. The dose administered in this study was adjusted according to body weight and the actual dosage was usually based on rounding, which may contribute to some variability observed in the infusion phase. Generally, no significant difference of AUC_0-6h seemed reasonable.

In the final PPK model constructed in this study, AST levels were negatively correlated with CL of Bu in pediatrics and CL declined 38.34% when AST increased from 12.7 to 127.4 U/L. This is unsurprising because busulfan is mainly eliminated by the liver as previously described. AST marks hepatocyte lesion, meaning that the higher the AST level, the more severe the liver damage. Consequently, patients with higher AST have a decreased ability to metabolize busulfan.

This study divided the patients into two subpopulations (malignant and nonmalignant diseases), which had no significant effect on the pharmacokinetic parameters of Bu. The conclusion was consistent with previous studies.^42,43 In fact, malignant diseases could be subdivided into seven types according to pathology, and nonmalignant diseases included six types. Nevertheless, some of the pathologies, such as osteopetrosis, were rare and present in only a small number of patients, resulting in an impossible comprehensive evaluation of the impact of diseases on Bu PK. Hence, the number of different types of cases needs to be expanded in order to analyze the specific impact of every primary disease on Bu PK.

The medication interaction with fludarabine and busulfan has not been widely confirmed. Although Yeh et al⁴⁴ reported that CL of IV Bu decreased significantly in patients receiving concomitant fludarabine administration, majority of studies did not observe changes in busulfan CL between patients with or without fludarabine.^45,46 Not surprisingly, in this PPK model, the inclusion of fludarabine as a covariate failed to significantly reduce the value of −2LL during the forward selection. In fact, patients in our cohort received one of the two regimens, Bu/CTX/FLU and Bu/CTX/Ara-C, as basic conditioning therapies. Although we only regarded FLU as a candidate covariate, in fact the influence of Ara-C has been also reflected in the model. When FLU=1, the impact of FLU was added to the model, and when FLU=0, the effect of Ara-C was included.

BSA was the most predictive covariate for CL and V, explaining 25.50% and 24.17% of the observed IIV, respectively. Weight-based dosing schedules were calculated with five fixed doses, however the model established by this study showed that body weight is not a significant predictor of Bu PK in children. In addition, in a retrospective study, SOS and early infectious complications occurred more frequently in the weight-based dosing group.⁴⁷ Besides, a PPK model, developed among patients of all ages, revealed that the maturation of Bu clearance reaches half of adult values at six weeks after birth.⁴² Furthermore, in children, Bu concentration did not show an obvious trend of change with postnatal age, which strongly supported the conclusion drawn from our cohort that BSA explained interindividual difference better than age. According to the PPK model established with Chinese pediatric patients, BSA-based dosing scheme was recommended. When compared with weight-based dosage, the new dosing scheme not only took the influence of GSTA1 mutation into consideration, but could also be applied to patients with varying degrees of liver damage. More importantly, in weight-based dosing regimens, only eight patients reached the window of 900–1350 μM·min. The validation analysis demonstrated that the new BSA-based dosage enabled almost all of the simulated patients to achieve the target AUC_0-6h, and greatly improved the success rate of obtaining the therapeutic window in real patients. Since this new dosing regimen was based on a retrospective analysis, a prospective study is necessary to confirm the benefits in terms of efficacy and safety.

Based on the final PPK model, Bu Bayesian estimation of individual AUC_0-6h values were performed by using various combinations of two-to-four sampling times within six hours following busulfan administration. Until this point, several LSS strategies have been established to predict BU exposure. Most LSS strategies published were estimated by using the trapezoidal rule (TR) or multiple linear regression (MLR), which usually reduced the accuracy of estimation. For example, Vaughan et al⁴⁸ concluded that four-to-five sampling points (3, 4, 5, and 6 h or 2, 3, 4, 5, and 6 h after dosing) could predict well in adult patients receiving IV Bu four times daily. Teitelbaum et al⁴⁹ developed an LSS with four sampling points (0, 2, 3, and 4 h after the start of the second infusion). To achieve a more accurate estimate of AUC of every LSS, PPK models, considered to be the gold standard, should be applied.

Theoretically, there is no absorption phase of IV Bu and C_max usually occurs at the infusion end. During the actual clinical nursing process, the infusion pump completes the titration in two hours, but there is still a little of the drug solution in the infusion tube, which may lead to the delay of T_max (the time to peak concentration). Assuming the peak concentration fell within 2~2.5 h, we tried to choose optimal T_max among 2, 2.25, and 2.5 h. C_4h and C_6h were regarded as time points of elimination phase. Since the metabolism of busulfan conformed to the one-compartment model, selecting C_4h or C_6h on behalf of elimination phase had a similar predictive function. Model 9 (C_2h, C_2.5h, and C_4h) and Model 10 (C_2h, C_2.5h, and C_6h) elevated the accuracy and precision of prediction while increasing medical costs and pain for children. Model 13 (C_2h and C_4h) not only had a better predictive performance by rRMSE, MAPE and Bland–Altman analysis, but was also more in line with the clinical requirement of reducing sampling points for TDM. Finally, considering the accuracy of prediction and the feasibility of pediatric clinical practice synthetically, Model 13 was selected as the optimal LSS.

The relationship between pharmacokinetics and outcome or toxicity of busulfan reported previously are summarized in Supporting Information Table 6. It can be seen that conclusions drawn out from different research studies were discrepant. In the current study, no relationship between busulfan PK and toxicity was observed (Supporting Information Figure 4).

In terms of ANC recovery and survival rate, there were significant differences between patients with two GSTA1 genotypes (p<0.05) in our study (Figure 7). Until now, no research has reported the relationship between GST mutations and engraftment. Ansari et al reported that individuals with poorer metabolism of Bu had also lower overall survival (OS) and GSTA1 *B haploid was associated with higher incidence of treatment-related toxicity.¹³ To confirm the long-term efficacy and toxicity of busulfan and explore their influencing factors, research with a larger sample size would be necessary.

Conclusion

In summary, this study was the first PPK model for busulfan that successfully incorporated GSTA1 genotypes in a Chinese pediatric population, suggesting that genotypes would be necessary for optimization of pediatric Bu treatment. Moreover, a BSA-based dosing regimen has been recommended for individual busulfan therapy in order to weaken variability in busulfan exposure and to enhance the safety and efficacy of Bu treatment. Finally, an optimal LSS (C_2h and C_4h) would be convenient for TDM in the future.