Abstract
Aim : The aim of this study was to develop a population pharmacokinetic (PPK) model in Chinese children for intravenous busulfan, and to develop a novel busulfan dosing regimen to support better area under the concentration-time curve (AUC) targeting.
Methods : We collected busulfan concentration-time samples from 69 children who received intravenous busulfan prior to allogeneic hematopoietic stem cell transplantation (allo-HSCT). A population pharmacokinetic model for busulfan was developed by nonlinear mixed effect modelling and was validated by an external dataset (n=14). A novel busulfan dosing regimen was developed through simulation on 1000 patients. Limited Sampling Strategy (LSS) was established by the Bayesian forecasting. Absolute Prediction Error (APE), Mean Absolute Prediction Error (MAPE) and relative Root Mean Squared Error (rRMSE) were calculated to evaluate predictive accuracy.
Results : A one-compartment model with first-order elimination best described the data. GSTA1 genotypes, BSA and AST were found to be significant covariates of Bu clearance, and BSA had remarkable impact of the volume. Moreover, recommended dose regimens for children with different GSTA1 genotypes and BSA were developed with a perfect AUC targeting. A two-point LSS, two hours and four hours after dosing, behaved well with acceptable prediction precision.
Conclusion : This study developed a PPK model for busulfan that firstly incorporated GSTA1 genotypes in Chinese pediatric population. We recommend a BSA-based dosing for personalizing busulfan therapy in pediatric population. Additionally, an optimal LSS (C2h and C4h) provides convenience for therapeutic drug monitoring (TDM) in the future.
KEYWORDS : busulfan, HSCT, individualized therapy, population pharmacokinetics, pediatric patients
Introduction
Busulfan (Bu) is a common alkylating agent, which can bind to the guanine of intracellular DNA to damage its structure and function (1). Bu-based conditioning schemes are regarded as the cornerstone of allogeneic hematopoietic stem cell transplantation (allo-HSCT) because of their myeloablative activity (1, 2). However, IV Bu usually leads to large variability in pharmacokinetics and the treatment window is narrow. Higher exposure (expressed as area-under-the-curve of 0 - 6 h, AUC0-6h) (>1350 μM·min) or low AUC0-6h (< 900 μM·min) of Bu may lead to a higher probability of sinusoidal obstructive syndrome (SOS) or acute graft-versus-host disease (aGVHD) (3-5). The exposure and clinical outcomes of Bu are largely affected by different conditioning regimen agents, which has been confirmed in children receiving Bu/cyclophosphamide (Bu/CTX), Bu/fludarabine (Bu/FLU) and Bu/ cyclophosphamide/etoposide (Bu/CTX/VP16) (6-8). Thus, therapeutic drug monitoring (TDM) for busulfan is the integral component of HSCT. Personalized Bu dosing via TDM based on the first-dose pharmacokinetics (PK) could make contribution to lower the occurrence of toxicity and to improve rate of engraftment (9, 10).
Bu is metabolized by the formation of a glutathione conjugate in the liver (11, 12). This reaction is primarily catalyzed by glutathione S-transferase (GST) enzymes, such as GSTA1 ,GSTM1 , GSTT1 and GSTP1 (13).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 GSTT1 and GSTP1 is relatively weak (2, 14). Hence, polymorphisms in the GSTA1 or GSTM1 isoenzymes would more likely affect Bu metabolism.
Studies about the relationships between GSTA1 polymorphisms and Bu PK have yielded inconsistent results. Although it has been reported thatGSTA1 haplotypes were not significant influence factors of intravenous (IV) Bu clearance (15), Yin et al. (16) stated adult patients with the GSTA1 *A/*B genotype had a significantly higher AUC, higher peak concentration (Cmax) and lower clearance (CL). Pediatric patients have an increased Bu clearance when compared with adults (17), which is partially caused by higher expression or activity of GST enzymes in children (18, 19). Research on the relationship between GSTA1 polymorphism and PK in children is also controversial. Ansari et al. (20) reported thatGSTA1 *A/*A was associated with lower drug concentration and more extensive metabolism. Conversely, Zwaveling et al. (21) asserted that GSTA1 polymorphisms had no influence on population PK parameters of IV Bu in children undergoing HSCT. In addition to age, the activity of GST enzymes is also affected by race. Minor allele frequencies (MAF) taken from HapMap were reported to be differing between Caucasian and Asian populations. Until now, the correlation between GST polymorphisms and PK in Chinese children has not been reported.
It has been suggested that, apart from GSTs polymorphism, the effect of age, body-weight (WT), primary disease, hepatic function, renal function and drug interactions (fludarabine and phenytoin) may partly explain inter-individual variability on Bu PK (22-25). Little is known about correlates of Bu PK in Chinese children, which may be potential factors to optimize dosage of Bu.
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. Furthermore, Maximum A Posteriori (MAP) Bayesian forecasting made use of a PPK model and limited number of samples to forecast AUC0-6h and to formulate an optimal LSS, which is an alternative monitoring strategy.
Methods
2.1 Patients and treatment regimens
From March 2019 to April 2020, we collected 76 patients received allo-HSCT. This study was approved by the Beijing’s Children Hospital and all patients/parents provided informed consent.
In the pre-treatment regimen, intravenous infusion of busulfan (BUSULFEX, Otsuka Pharmaceutical Co. LTD, Zhejiang, China) typically started eight days prior to transplantation. The infusion frequency was once every six hours within 3-4 days and each infusion process last for two hours. 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-23 kg, 0.95 mg/kg for children between 23-34 kg, and 0.8 mg/kg for children more than 34 kg) in accordance with the European Medicines Agency (EMA) (26). 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 non-malignant tumors. Specific dosage and duration of different conditioning regimens have been exhibited in Supporting Information Table 1 . Phenytoin (PHT) started 30 minutes 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.
2.2 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) (1). The lower limit of quantitation was 10 ng/mL and the range of quantitation was from 10 to 10000 ng/ml. Pre-transplant genomic DNA was isolated and extracted from whole blood or hemocyte prior to the first Bu infusion. GST genotypes of patients,GSTA1 (rs3957356 and rs3957357, which defines haplotype *Aand *B ) and GSTM1 (rs3754446), were detected with the ABI 3730XL DNA Analyzer (ABICo.; BioSune Biotechnology Co., Ltd, Shanghai, China). Supporting Information Table 2 displays the primer sets and Tm used for the genotyping assays.
2.3 PPK analysis
BU plasma concentration-time data was analyzed using a Non-Linear Mixed-Effects (NLME) model implemented in Phoenix 8.0 (Certara USA Inc., Princeton, NJ, USA). Typical PPK parameters and their random inter-individual 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 Pi = P × \(e^{\eta i}\), wherePi is the individual PK parameter of thei th 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 intra-individual (residual) variability. 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 versus covariates and clinical plausibility, we selected candidate covariates for each PK parameters (28). The influence of potential covariates [sex, age, body weight, body surface area (BSA), alanine transaminase (ALT), aspartate aminotransferase (AST), alkaline phosphatase (ALP), total bilirubin (TBIL), creatinine (CRE), creatinine clearance (CLcr), GSTA1genotypes (*A/*A, *A/*B, *B/*B ), GSTM5 genotypes, primary disease (malignant disease, non-malignant 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 (27). 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.
2.4 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(29).
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 value (IPRED) (30). 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% confidence interval (CI) for concentrations by simulating 1,000 individuals. Simulated concentrations were compared with the 5th, 50th, 95th percentiles of the observed concentrations. In the bootstrap analysis (31), 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 (32).
2.5 Dosing Simulations
A new dosage scheme using a simulated dataset of 1000 patients was designed to achieve a targeted AUC of 1125 μM·min. The BSA of simulated patients were distributed from 0.2 m2 to 1.6 m2. Meanwhile, the genotypes of GSTA1 included *A*A and *A*Band the range of AST was coincident with the clinical dataset. Then, 1000 AUCsimulated values and their variabilities were generated in Phoenix-NLME. Model-based doses were recommended according to the following Equation 1 :