\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%

\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[greek,ngerman,english]{babel}



\usepackage{float}








\begin{document}

\title{Optimization of busulfan dosing regimen in pediatric patients using
population pharmacokinetic model incorporating GST polymorphisms}



\author[1]{Jin-jie Yuan}%
\author[1]{Ning Sun}%
\author[1]{Xin-ying Feng}%
\author[2]{Huan He}%
\author[1]{Guang-hua Zhu}%
\author[3]{Li-Bo Zhao}%
\affil[1]{Affiliation not available}%
\affil[2]{Beijing Children’s Hospital, Capital Medical University, National Center for Children’s Health, Beijing 10045, China}%
\affil[3]{Peking University}%


\vspace{-1em}



  
  \date{May 16, 2024}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{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 an Asian 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.%
\end{abstract}%



\sloppy


\textbf{Abstract}

\textbf{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.

\textbf{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.

\textbf{Results} : A one-compartment model with first-order elimination
best described the data. \emph{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 \emph{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.

\textbf{Conclusion} : This study developed a PPK model for busulfan that
firstly incorporated \emph{GSTA1} genotypes in Chinese pediatric
population. We recommend a BSA-based dosing for personalizing busulfan
therapy in pediatric population. Additionally, an optimal LSS
(C\textsubscript{2h} and C\textsubscript{4h}) provides convenience for
therapeutic drug monitoring (TDM) in the future.

\textbf{KEYWORDS} : busulfan, HSCT, individualized therapy, population
pharmacokinetics, pediatric patients

\textbf{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,
AUC\textsubscript{0-6h}) (\textgreater{}1350 \selectlanguage{greek}μ\selectlanguage{english}M\selectlanguage{ngerman}·min) or low
AUC\textsubscript{0-6h} (\textless{} 900 \selectlanguage{greek}μ\selectlanguage{english}M\selectlanguage{ngerman}·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 \emph{GSTA1} ,\emph{GSTM1} ,
\emph{GSTT1} and \emph{GSTP1} (13).\emph{GSTA1} is the predominant GST
enzyme involved in Bu metabolism, and the activity of \emph{GSTM1} is
close to the half of \emph{GSTA1} . However, the activity of
\emph{GSTT1} and \emph{GSTP1} is relatively weak (2, 14). Hence,
polymorphisms in the \emph{GSTA1} or \emph{GSTM1} isoenzymes would more
likely affect Bu metabolism.

Studies about the relationships between \emph{GSTA1} polymorphisms and
Bu PK have yielded inconsistent results. Although it has been reported
that\emph{GSTA1} haplotypes were not significant influence factors of
intravenous (IV) Bu clearance (15), Yin \emph{et al.} (16) stated adult
patients with the \emph{GSTA1} \emph{*A/*B} genotype had a significantly
higher AUC, higher peak concentration (C\textsubscript{max}) 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 \emph{GST} enzymes in children (18, 19).
Research on the relationship between \emph{GSTA1} polymorphism and PK in
children is also controversial. Ansari \emph{et al.} (20) reported
that\emph{GSTA1 *A/*A} was associated with lower drug concentration and
more extensive metabolism. Conversely, Zwaveling \emph{et al.} (21)
asserted that \emph{GSTA1} polymorphisms had no influence on population
PK parameters of IV Bu in children undergoing HSCT. In addition to age,
the activity of \emph{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 \emph{GST} polymorphisms and PK in Chinese children has not been
reported.

It has been suggested that, apart from \emph{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
AUC\textsubscript{0-6h} and to formulate an optimal LSS, which is an
alternative monitoring strategy.

\textbf{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
\textbf{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,\emph{GSTA1} (rs3957356 and rs3957357, which defines haplotype
\emph{*A}and \emph{*B} ) and \emph{GSTM1} (rs3754446), were detected
with the ABI 3730XL DNA Analyzer (ABICo.; BioSune Biotechnology Co.,
Ltd, Shanghai, China). \textbf{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 \emph{P\textsubscript{i} = P ×} \(e^{\eta i}\),
where\emph{P\textsubscript{i}} is the individual PK parameter of
the\emph{i} th individual, \emph{P} is the geometric average population
value, and \selectlanguage{greek}\emph{η\textsubscript{i}}\selectlanguage{english} is the subject-specific random
effect value, normally distributed random variable with a mean of 0 and
a variance of \selectlanguage{greek}ω\selectlanguage{english}\textsuperscript{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),
\emph{GSTA1}genotypes (\emph{*A/*A, *A/*B, *B/*B} ), \emph{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 (\emph{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 (\emph{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 \emph{GSTA1} genotypes of patients. Differences of
\textless{} 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 \selectlanguage{greek}μ\selectlanguage{english}M\selectlanguage{ngerman}·min. The BSA of simulated
patients were distributed from 0.2 m\textsuperscript{2} to 1.6
m\textsuperscript{2}. Meanwhile, the genotypes of \emph{GSTA1} included
\emph{*A*A and *A*B}and the range of AST was coincident with the
clinical dataset. Then, 1000 AUC\textsubscript{simulated} values and
their variabilities were generated in Phoenix-NLME. Model-based doses
were recommended according to the following \textbf{Equation 1} :\selectlanguage{english}
\begin{longtable}[]{@{}ll@{}}
\toprule
Dose\textsubscript{recommended} = 1125 / AUC\textsubscript{simulated} \selectlanguage{ngerman}×
Dose\textsubscript{simulated} & (1)\tabularnewline
\bottomrule
\end{longtable}

Then, AUC values were calculated according to the new recommended doses.
Variabilities of AUC 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 \selectlanguage{greek}μ\selectlanguage{english}M\selectlanguage{ngerman}·min) was also listed.

2.6 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 (e.g. BSA, AST, genotypes) influence the PK of BU. After
evaluating final model predictive performance, filtrating optimal LSS
was considered to optimize the sampling scheme for TDM. Candidate
sampling points included: (1) 2, 2.25, 4 and 6 hours; (2) 2, 2.5, 4 and
6 hours; (3) 2.25, 2.5, 4 and 6 hours; (4) 2, 4 and 6 hours; (5) 2.25, 4
and 6 hours; (6) 2.5, 4 and 6 hours; (7) 2, 2.25 and 4 hours; (8) 2,
2.25 and 6 hours; (9) 2, 2.5 and 4 hours; (10) 2, 2.5 and 6 hours; (11)
2.25, 2.5 and 4 hours; (12) 2.25, 2.5 and 6 hours; (13) 2 and 4 hours;
(14) 2 and 6 hours; (15) 2.25 and 4 hours; (16) 2.25 and 6 hours; (17)
2.5 and 4 hours; (18) 2.5 and 6 hours after the beginning of first
infusion. The blood drug concentration at 2.25 hours after
administration was simulated by the final model. Each candidate LSS was
tested by Bayesian method to obtain AUC\textsubscript{0-6h} of each
patient. APE, MAPE and rRMSE were calculated by comparing the predicted
AUC\textsubscript{0-6h}derived from the LSS with the actual
AUC\textsubscript{0-6h} obtained from each patient's full
concentration-time samples, as according to\textbf{Equation 2, 3 and 4}
to evaluate predictive accuracy:\selectlanguage{english}
\begin{longtable}[]{@{}ll@{}}
\toprule
APE (\%) =\(\ \frac{|predicted\ \text{AUC}_{0-6h}-\text{actual}\ \text{AUC}_{0-6h}|\ }{\text{actual}\ \text{AUC}_{0-6h}}\ \times\ 100\)\% & (2)\tabularnewline
\midrule
\endhead
MAPE (\%) =\(\ \frac{1}{N}\sum{(\frac{|predicted\ \text{AUC}_{0-6h}-\text{actual}\ \text{AUC}_{0-6h}\ |}{\text{actual}\ \text{AUC}_{0-6h}}\ \times\ 100\%)}\) & (3)\tabularnewline
rRMSE (\%) = \(\frac{1}{N}\sqrt{\sum{(\frac{\text{predicted\ }\text{AUC}_{0-6h}-\text{actual}\ \text{AUC}_{0-6h}\ }{\text{actual}\ \text{AUC}_{0-6h}}\ \times\ 100\%)}^{2}}\) & (4)\tabularnewline
\bottomrule
\end{longtable}

The lower the MAPE value, the better the LSS. A sampling strategy was
considered to display a good predictive performance when the 95\%
confidence interval around the MAPE was less than 20\% of the reference
Bu AUC\textsubscript{0-6h} values (33). Bland-Altman plots were applied
to assess the agreement between the predicted and the actual
AUC\textsubscript{0-6h}, and the difference was expressed as a mean \selectlanguage{ngerman}±
1.96 SD.

2.7 Definition of Clinical Outcomes

The definition of engraftment was absolute neutrophil count (ANC) \selectlanguage{english}[?] 0.5
x 10\textsuperscript{9} per liter for the first three days in a row and
platelet count (PLT) \selectlanguage{english}[?] 20 x 10\textsuperscript{9} per liter for the
first one week after HSCT (10, 16). Primary engraftment failure was
defined as failing to reach an ANC of 0.5 x 10\textsuperscript{9} /L
within 30 days after HSCT (16). In brief, the diagnosis of SOS was
according to the modified Seattle criteria (34). The Mount Sinai Acute
GVHD International Consortium (MAGIC) was taken as the diagnostic
criterion for aGVHD (35).

2.8 Statistical Analysis

Continuous data, like pharmacokinetics parameters, was compared using a
two-tailed t-test, whereas time-to-event data was compared by the
log-rank test. Univariate analysis was performed with GraphPad Prism
(GraphPad Software, San Diego, CA). In this study, \emph{p} values
\textless{} 0.05 are considered as statistical significance, and\emph{p}
values \textless{} 0.1 were deemed as trends.

\textbf{Results}

3.1 Population characteristics

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

3.2 PPK model

A total of 398 blood samples were collected for the PPK
analysis.\textbf{Supporting Information 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 inter-individual 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
(\selectlanguage{greek}Δ-\selectlanguage{english}2LL \textgreater{} 20, \emph{p} \textless{} 0.001).

In the forward model building step, BSA, AST, types of primary diseases
and \emph{GSTA1} genotypes declined the value of -2LL by more than 3.84
(\emph{p} \textless{} 0.05) at each addition. The \emph{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 \emph{GSTM1} genotypes had no
significant effect on PK parameters.

During backward elimination, the value of -2LL increased significantly
(\selectlanguage{greek}Δ-\selectlanguage{english}2LL \textgreater{}6.63, \emph{p} \textless{} 0.01) when respectively
removing the \emph{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 (\selectlanguage{greek}Δ-\selectlanguage{english}2LL
= 6.333, \emph{p} \textgreater{} 0.01). \textbf{Table 3} presented the
detailed PK parameters of the base and final model. The CL and V of
busulfan were affected by the \emph{GSTA1} genotypes, BSA and AST as
illustrated by the following \textbf{Equation 5 and 6} :\selectlanguage{english}
\begin{longtable}[]{@{}ll@{}}
\toprule
CL (L/h) = CL\textsubscript{pop}\selectlanguage{ngerman}·\selectlanguage{greek}\(({\frac{\text{BSA}_{i}}{\text{BSA}_{\text{median}}})}^{covBSA(CL)}({\frac{\text{AST}_{i}}{\text{AST}_{\text{median}}})}^{covAST(CL)}e^{\text{cov}GSTA1(CL)GATA1}e^{\text{ηCL}}\selectlanguage{english}\) & (5)\tabularnewline
\midrule
\endhead
V (L) = V\textsubscript{pop}\(\ \\selectlanguage{greek} ({\frac{\text{BSA}_{i}}{\text{BSA}_{\text{median}}})}^{\text{covBSA}\left(V\right)}e^{\text{ηV}}\selectlanguage{english}\) & (6)\tabularnewline
\bottomrule
\end{longtable}

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

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

3.3 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 \emph{GSTA1} phenotype. As
demonstrated in the scatter plots of OBS \emph{vs.} PRED (\textbf{Figure
2A and 2B} ), 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 = \selectlanguage{ngerman}± 2) than that in the base model (\textbf{Figure 2C and 2D}
). As a whole, the final model was visually improved in terms of the GOF
and had a more accurate predictive performance.

The estimated values generated by bootstrap analysis were close to the
parameters in the final model, which also fell within 95\% CIs
(\textbf{Table 3} ). Thus, it can be considered that the final model was
accurate and robust. \textbf{Figure 3} 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 (\textbf{Figure 4} ).

3.4 New dosing strategy and simulations

To achieve a target AUC of 1125 \selectlanguage{greek}μ\selectlanguage{english}M\selectlanguage{ngerman}·min, a simple linear relationship
between the total doses and BSA is illustrated in \textbf{Figure 5(A).}
The linear \textbf{Equation 7 and 8} for \emph{GSTA1 *A/*A} and
\emph{*A*B} was respectively:\selectlanguage{english}
\begin{longtable}[]{@{}ll@{}}
\toprule
Dose \textsubscript{recommended} (mg) = 34.14 \selectlanguage{ngerman}× BSA
(m\textsuperscript{2}) + 3.75 & (7)\tabularnewline
\midrule
\endhead
Dose \textsubscript{recommended} (mg) = 30.99 \selectlanguage{ngerman}× BSA
(m\textsuperscript{2}) + 3.21 & (8)\tabularnewline
\bottomrule
\end{longtable}

The dosages were normalized by BSA (mg/m\textsuperscript{2}) due to the
fact that body size-based dosing is more familiar to pediatricians. As
demonstrated in \textbf{Figure 5 (B)} , the doses normalized by BSA were
shown to be a log-linear function of BSA. Dosages according to this
log-linear curve 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 based on
BSA normalization was listed in \textbf{Table 4} .

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

3.5 Limited sampling strategies

Eighteen models with different sampling time points are all listed
in\textbf{Table 5} . The relationship between the predicted and actual
AUC\textsubscript{0-6h} for these models is shown in \textbf{Supporting
Information Figure 2} . From two to four-point models, 83.33\% LSSs
fitted well with the correlation (\emph{r\textsuperscript{2}} ) of more
than 0.85. Prediction precision of LSSs expressed as rRMSE and MAPE is
also given in \textbf{Table 5} . Model 2 (C\textsubscript{2h},
C\textsubscript{2.5h}, C\textsubscript{4h} and C\textsubscript{6h})
showed not only the best fit to the Bu AUC\textsubscript{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\textsubscript{2h},
C\textsubscript{2.5h}and C\textsubscript{4h}) and Model 10
(C\textsubscript{2h}, C\textsubscript{2.5h} and C\textsubscript{6h})
both behaved well, and with similar prediction precision to Model 2.
Within these models, no patients had an AUC\textsubscript{0-6h} lower
than -15\% or higher than +15\% of the reference value. As for these
two-sampling schemes, Model 13 (C\textsubscript{2h} and
C\textsubscript{4h}), Model 14 (C\textsubscript{2h} and
C\textsubscript{6h}), Model 17 (C\textsubscript{2.5h} and
C\textsubscript{4h}) and Model 18 (C\textsubscript{2.5h} and
C\textsubscript{6h}) had relatively low rRMSE and MAPE. The
Bland--Altman plots have verified the excellent capacity of prediction
of the seven models above \textbf{(Supporting Information Figure 3)} .

3.6 Clinical outcomes

Among the clinical outcomes analyzed, 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\textsubscript{0-6h} of 659.9
and 482.7 \selectlanguage{greek}μ\selectlanguage{english}M\selectlanguage{ngerman}·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 \textbf{Figure 6} ,
there were significant differences in ANC recovery and survival rate
between patients with two \emph{GSTA1}genotypes (\emph{p} \textless{}
0.05). Seven patients (10.14\%) developed SOS and aGVHD I--IV was
documented in 57 patients (82.61\%), of whom seventeen had grade III--IV
disease (24.64\%). The clinical outcomes of patients with
two\emph{GSTA1} genotypes were compared in \textbf{Supporting
Information Table 4} .

In a univariate model, type of donor (HLA matched \emph{vs.} mismatched)
was the single variable that had significant association with the
incidence of aGVHD (\emph{p} = 0.022) (\textbf{Figure 7} ). The logistic
regression model was again put into use to investigate the relationship
between the incidence of aGVHD I--IV, aGVHD \selectlanguage{english}III/\selectlanguage{english}IV, SOS and Bu
AUC\textsubscript{0-6h}, while \textbf{Supporting Information Figure
4}revealed that no correlation between Bu AUC\textsubscript{0-6h} and
regimen-related toxicity or mortality was observed.

Discussion

This study developed the first PPK model for busulfan that successfully
incorporated \emph{GSTA1} genotypes in Chinese pediatric population and
partly explained the source of large variability of busulfan exposure,
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 polymorphisms (36). There
have been several researches supporting the positive association
between\emph{GSTs} polymorphisms and Bu clearance since 2016. In a
pediatric multicenter study, Ansari \emph{et al.} (37) reported that the
activity of \emph{GSTA1} promoter was significantly descended in the
case of\emph{*B} haplotype compared with\emph{*A} haplotype, thus the
correlation between \emph{GSTA1} genotypes and clearance was
distinguished. In addition, \emph{GSTA1} diplotypes with slow
metabolizing capacity were associated with higher incidence of SOS,
aGVHD and combined treatment-related toxicity. Another study succeeded
in incorporating the \emph{GST} genetic variants (\emph{GSTA1} ) into a
PPK model for Bu in a Caucasian pediatric population, and then tailored
the dose according to the individual metabolic capacity (38). For
Chinese adults, Yin \emph{et al.} (16) concluded that patients
with\emph{GSTA1} \emph{*A/*B} genotype possessed a higher
C\textsubscript{max}, higher AUC and lower clearance than the group
with\emph{GSTA1 *A/*A} genotype. Furthermore, \emph{GSTA1} expression in
young children has been reported to be higher than adults (19). In a
recent study of Japanese pediatric population (n=20), Nishikawa \emph{et
al.} (39) stated the correlation between \emph{GST} polymorphisms and
clearance was distinguished.

In our two cohorts (n=84), the frequency of \emph{GSTA1}
\emph{*B}haplotype (\emph{GSTA1-52T/-69A} ) was 11.9\% (MAF = 0.119). In
other words, a minority of patients had the \emph{*B} haplotype, of
which 1.19\% were homozygous (\emph{GSTA1 *B/*B} ). Therefore, the MAF
of Chinese patients was lower than that of global population (MAF =
0.306) taken from 1000 Genomes. Due to the significant discrepancy of
the frequency of haplotype \emph{*B} between different populations,
studies in Caucasians may not be able to represent that in the Asian
population. Besides, in a study conducted in patients with acute myeloid
leukemia, Yee \emph{et al.} (40) found that C allele in the \emph{GSTM1}
locus (rs3754446) was associated with decreased Bu AUC of first dose and
lower disease-free survival. Thus,\emph{GSTA1} and \emph{GSTM1}
genotypes were taken into consideration in the present study when
exploring the influence of \emph{GST}polymorphisms 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 of CL was 4.79 L/h and of V was 14.8
L, consistent with those results in previous studies (38). The PPK model
demonstrated that the \emph{GSTA1} polymorphisms were associated with Bu
clearance. Patients carrying the \emph{GSTA1 *A/*B} genotype had a
17.3\% lower clearance than those carrying \emph{GSTA1 *A/*A} , which is
consistent with previous studies reporting that the presence of mutation
allele probably resulted in the decreased activity of GST enzyme (11,
37). However, genetic variation in \emph{GSTM1} showed no significant
impact on Bu CL, likely because the function of the \emph{GSTM1}enzyme
involved in Bu metabolism was less than the \emph{GSTA1} enzyme.

As shown in \textbf{Figure 1} , there was a significant difference in
clearance per body surface area between the two \emph{GSTA1} genotypes,
but there was no significant difference in AUC per body surface area.
These results were indeed confusing. As is well-known, AUC of drug
administered intravenously is affected by the distribution, metabolism
and excretion of drug, while CL is influenced only by metabolism and
excretion. There are many complex physiological factors that affect drug
distribution, especially in children. Thus, AUC may be not completely
inversely proportional to CL. As shown in \textbf{Supporting Information
Figure 1} , the distribution phase (0-2 hours after dosing) shown large
heterogeneity between two genotypes, which may slightly interfere the
statistical~analysis~of AUC. Nevertheless, in the eliminate phase, it
was obvious that patients with\emph{GSTA1 *A*A} appeared faster
clearance than those with \emph{GSTA1 *A*B} . Generally, the non-obvious
significance of AUC\textsubscript{0-6h} seemed reasonable. Moreover, the
number of patients treated with busulfan with an AUC value of more than
900 \selectlanguage{greek}μ\selectlanguage{english}M\selectlanguage{ngerman}·min was eight. In other words, 88.4\% of children with HSCT could
not reach the range of 900-1350 \selectlanguage{greek}μ\selectlanguage{english}M\selectlanguage{ngerman}·min recommended by FDA. The main
factor resulting in the low AUC\textsubscript{0-6h} may be the fact
that, for Chinese patients, the frequency of \emph{GSTA1 *A} with high
metabolizing capacity was much higher when compared to \emph{GSTA1 *B}
with slow metabolizing capacity. Thus, the relatively low level of
AUC\textsubscript{0-6h} in this study population was rational.

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. This study divided the patients into two subpopulations
(malignant and non-malignant diseases), which had no significant effect
on the pharmacokinetic parameters of Bu. In fact, malignant diseases
could be subdivided into seven types according to pathology, and
non-malignant 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 variability of Bu PK. In oral busulfan-based pediatric
research, the influence of underlying diseases on busulfan disposition
was significant and CL/F was significantly lower in group with immune
deficiencies than other groups (metabolic diseases, hemoglobinopathies
and hematological malignancies) (25). 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.

There are theoretical and documented medication interaction with
fludarabine and busulfan (41). Clearance of IV Bu decreased
significantly in patients receiving concomitant fludarabine
administration (\emph{p} = 0.0016) and the average of reduction was
9.7\% (42). However, there was no other studies drawing similar
conclusions (43, 44). Not surprisingly, in this PPK model, the inclusion
of fludarabine as a covariate failed to significantly decline 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. Due to the use of Bu/CTX in all patients without
exception, the influence of CTX has already existed in the model.
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. Additionally, in a retrospective study,
SOS and early infectious complications occurred more frequently in the
weight-based dosing group (45). Besides, a PPK model, developed among
patients of all ages, revealed that the maturation of Bu clearance
reaches half of adult values at 6 weeks after birth (46). Also, 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 age was not a significant factor affecting Bu clearance.
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\emph{GSTA1}
polymorphism into consideration, but could also be applied to patients
with different liver function. More importantly, the simulation analysis
demonstrated that the three fixed doses given on an
mg/m\textsuperscript{2} basis enabled almost all of the young patients
(0.2\textasciitilde{}1.6 m\textsuperscript{2}) to achieve the target
AUC\textsubscript{0-6h}. 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\textsubscript{0-6h} values were performed by using various
combinations of 2-4 sampling times within 6 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 and
lacked professionalism. For example, Vaughan \emph{et al.} (47)
concluded that 4-5 sampling points (3, 4, 5 and 6 hours or 2, 3, 4, 5
and 6 hours after dosing) could predict well in adult patients receiving
IV Bu four times daily. Teitelbaum \emph{et al.} (48) developed an LSS
with four sampling points (0, 2, 3 and 4 hours 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.
In this study population, C\textsubscript{max} generally appeared in
2-2.5 hours after dosing according to the concentration-time curve.
Therefore, peak concentration was picked among C\textsubscript{2h},
C\textsubscript{2.25h} or C\textsubscript{2.5h}, and
C\textsubscript{4h}and C\textsubscript{6h} were regarded as time points
of elimination phase. Since the metabolism of busulfan conformed to the
one-compartment model, selecting C\textsubscript{4h} or
C\textsubscript{6h} on behalf of elimination phase had a similar
predictive function. Model 9 (C\textsubscript{2h}, C\textsubscript{2.5h}
and C\textsubscript{4h}) and Model 10 (C\textsubscript{2h},
C\textsubscript{2.5h} and C\textsubscript{6h}) elevated the accuracy and
precision of prediction while increasing medical costs and pain for
children. Model 13 (C\textsubscript{2h} and C\textsubscript{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 \textbf{Supporting
Information Table 5} . A multicenter retrospective research stated that
graft-failure occurred more frequently in the low AUC group, meanwhile,
acute toxicity and TRM were significantly higher in the high AUC group
(5). Additionally, Andersson \emph{et al.} (49) found that the
probability of developing aGVHD increased with increasing AUC. Similar
correlations between toxicities and AUC have been verified in both
Korean children and adult groups (14, 50). However, no relationship
between busulfan PK and toxicity was observed. Based on the data of 27
children with sickle cell disease undergoing HSCT, hepatic toxicity
(SOS) was not associated with busulfan AUC (51). Moreover,
Jessica\emph{et al} . (52) concluded that there was no significant
association between AUC \textsubscript{dose1} and death, relapse, or a
composite of the two.

In terms of ANC recovery and survival rate, there were significant
differences between patients with two GSTA1 genotypes
(\emph{p}\textless{} 0.05) in our study. Till now, no research has been
reported concerning the relationship between \emph{GSTs} polymorphisms
and engraftment or mortality. However, Ansari \emph{et al.} reported
that GSTA1 *B haploid was associated with higher incidence of
treatment-related toxicity (37). Hence, we deduced that the high
incidence of toxicity was likely connected to the postponement and the
low rate of ANC recovery, which further led to higher chance of
infection, the main cause of death in this study. In the next stage, we
will include more abundant samples to verify the above conjecture.

In this study population, type of donor was the major predictor of the
occurrence of aGVHD, whereas neither \emph{GSTM1} genotypes nor
pharmacokinetic parameters were found to be significantly correlative to
SOS. These results are difficult to interpret because multiple
pre-transplantation and transplantation-related factors have been
implicated in the rate of engraftment and the incidence of toxic effect.
So far, no single factor has been reported to result in aGVHD or SOS
independently. On one hand, it is worth noting that the homogeneity of
patients' disease would be beneficial to the correlation analysis. For
example, Srivastava at al. (53) concluded that GSTM1-null genotype was
relevant to the incidence of SOS in 114 patients undergoing HSCT with
uniform disease, \selectlanguage{greek}β-\selectlanguage{english}thalassemia. On the other hand, the influence of
complex combined medicines could not be ignored. Sixty-eight out of 69
patients adopted a combination of three or more myeloablative drugs and
57.97\% of pediatric patients used four drugs with myeloablative
activity in succession, which led that the AUC of a single drug, Bu, may
be not the only influencing factor in curative effect or toxic reaction.
It should be also recognized that when Bu was combined with different
agents, relationships between Bu PK and outcome or toxicity would be
varied (54), as has been confirmed in BU/CY/TBI regimens (55).
Consequently, in order to clarify the correlation, analysis should be
conducted according to disease types and conditioning regimens in the
further.

\textbf{CONCLUSION}

In summary, this study was the first PPK model for busulfan that
successfully incorporated \emph{GSTA1} genotypes in a Chinese pediatric
population and partly explained the source of large variability of
busulfan exposure, suggesting that TDM would be necessary for
optimization of pediatric Bu treatment. Moreover, a BSA-based dosing
regimen is 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\textsubscript{2h}
and C\textsubscript{4h}) would be convenient for TDM in the future.

\textbf{REFERENCES}

1. Jinjie Y, Sun N, Zhang S, Feng X, Chen X, Zhao D, et al. A rapid
HPLC-MS/MS method for determining busulfan in hemolysis samples from
children with hematopoietic stem cell transplantation. Biomedical
chromatography : BMC. 2020:e4898.

2. ten Brink MH, Zwaveling J, Swen JJ, Bredius RG, Lankester AC,
Guchelaar HJ. Personalized busulfan and treosulfan conditioning for
pediatric stem cell transplantation: the role of pharmacogenetics and
pharmacokinetics. Drug Discov Today. 2014;19(10):1572-86.

3. Dix SP, Wingard JR, Mullins RE, Jerkunica I, Davidson TG, Gilmore CE,
et al. Association of busulfan area under the curve with veno-occlusive
disease following BMT. Bone Marrow Transplant. 1996;17(2):225-30.

4. Geddes M, Kangarloo SB, Naveed F, Quinlan D, Chaudhry MA, Stewart D,
et al. High busulfan exposure is associated with worse outcomes in a
daily i.v. busulfan and fludarabine allogeneic transplant regimen.
Biology of blood and marrow transplantation : journal of the American
Society for Blood and Marrow Transplantation. 2008;14(2):220-8.

5. Bartelink IH LA, Elisabeth ML, \emph{et al} . Association of busulfan
exposure with survival and toxicity after haemopoietic cell
transplantation in children and young adults: a multicentre,
retrospective cohort analysis. \emph{Lancet Haematol} .
2016;3(11):e526--e36.

6. Ayala E, Figueroa J, Perkins J, Kim J, Yue B, Riches M, et al.
Myeloablative intravenous pharmacokinetically targeted busulfan plus
fludarabine as conditioning for allogeneic hematopoietic cell
transplantation in patients with non-Hodgkin lymphoma. Clinical
lymphoma, myeloma \& leukemia. 2015;15(6):335-40.

7. Parmar S, Rondon G, de Lima M, Thall P, Bassett R, Anderlini P, et
al. Dose intensification of busulfan in the preparative regimen is
associated with improved survival: a phase I/II controlled, randomized
study. Biology of blood and marrow transplantation : journal of the
American Society for Blood and Marrow Transplantation.
2013;19(3):474-80.

8. Andersson BS, Valdez BC, de Lima M, Wang X, Thall PF, Worth LL, et
al. Clofarabine \selectlanguage{ngerman}± fludarabine with once daily i.v. busulfan as
pretransplant conditioning therapy for advanced myeloid leukemia and
MDS. Biology of blood and marrow transplantation : journal of the
American Society for Blood and Marrow Transplantation.
2011;17(6):893-900.

9. Malar R, Sjoo F, Rentsch K, Hassan M, Gungor T. Therapeutic drug
monitoring is essential for intravenous busulfan therapy in pediatric
hematopoietic stem cell recipients. Pediatr Transplant.
2011;15(6):580-8.

10. Wall DA, Chan KW, Nieder ML, Hayashi RJ, Yeager AM, Kadota R, et al.
Safety, efficacy, and pharmacokinetics of intravenous busulfan in
children undergoing allogeneic hematopoietic stem cell transplantation.
Pediatr Blood Cancer. 2010;54(2):291-8.

11. Czerwinski M, Gibbs JP, Slattery JT. Busulfan conjugation by
glutathione S-transferases alpha, mu, and pi. Drug metabolism and
disposition: the biological fate of chemicals. 1996;24(9):1015-9.

12. Gibbs JP, Czerwinski M, Slattery JT. Busulfan-glutathione
conjugation catalyzed by human liver cytosolic glutathione
S-transferases. Cancer research. 1996;56(16):3678-81.

13. Huezo-Diaz P, Uppugunduri CR, Tyagi AK, Krajinovic M, Ansari M.
Pharmacogenetic aspects of drug metabolizing enzymes in busulfan based
conditioning prior to allogenic hematopoietic stem cell transplantation
in children. Current drug metabolism. 2014;15(3):251-64.

14. Choi B, Kim MG, Han N, Kim T, Ji E, Park S, et al. Population
pharmacokinetics and pharmacodynamics of busulfan with GSTA1
polymorphisms in patients undergoing allogeneic hematopoietic stem cell
transplantation. Pharmacogenomics. 2015;16(14):1585-94.

15. Abbasi N, Vadnais B, Knutson JA, Blough DK, Kelly EJ, O'Donnell PV,
et al. Pharmacogenetics of intravenous and oral busulfan in
hematopoietic cell transplant recipients. Journal of clinical
pharmacology. 2011;51(10):1429-38.

16. Yin J, Xiao Y, Zheng H, Zhang YC. Once-daily i.v. BU-based
conditioning regimen before allogeneic hematopoietic SCT: a study of
influence of GST gene polymorphisms on BU pharmacokinetics and clinical
outcomes in Chinese patients. Bone Marrow Transplant.
2015;50(5):696-705.

17. Hassan M, Ljungman P, Bolme P, Ringden O, Syruckova Z, Bekassy A, et
al. Busulfan bioavailability. Blood. 1994;84(7):2144-50.

18. Gibbs JP, Murray G, Risler L, Chien JY, Dev R, Slattery JT.
Age-dependent tetrahydrothiophenium ion formation in young children and
adults receiving high-dose busulfan. Cancer research.
1997;57(24):5509-16.

19. Gibbs JP, Liacouras CA, Baldassano RN, Slattery JT. Up-regulation of
glutathione S-transferase activity in enterocytes of young children.
Drug metabolism and disposition: the biological fate of chemicals.
1999;27(12):1466-9.

20. Ansari M, Rezgui MA, Theoret Y, Uppugunduri CR, Mezziani S, Vachon
MF, et al. Glutathione S-transferase gene variations influence BU
pharmacokinetics and outcome of hematopoietic SCT in pediatric patients.
Bone Marrow Transplant. 2013;48(7):939-46.

21. Zwaveling J, Press RR, Bredius RG, van Derstraaten TR, den Hartigh
J, Bartelink IH, et al. Glutathione S-transferase polymorphisms are not
associated with population pharmacokinetic parameters of busulfan in
pediatric patients. Ther Drug Monit. 2008;30(4):504-10.

22. Kusama M, Kubota T, Matsukura Y, Matsuno K, Ogawa S, Kanda Y, et al.
Influence of glutathione S-transferase A1 polymorphism on the
pharmacokinetics of busulfan. Clin Chim Acta. 2006;368(1-2):93-8.

23. Rocha V, Porcher R, Fernandes JF, Filion A, Bittencourt H, Silva W,
Jr., et al. Association of drug metabolism gene polymorphisms with
toxicities, graft-versus-host disease and survival after HLA-identical
sibling hematopoietic stem cell transplantation for patients with
leukemia. Leukemia. 2009;23(3):545-56.

24. Vassal G, Deroussent A, Challine D, Hartmann O, Koscielny S,
Valteau-Couanet D, et al. Is 600 mg/m2 the appropriate dosage of
busulfan in children undergoing bone marrow transplantation? Blood.
1992;79(9):2475-9.

25. Bertholle-Bonnet V, Bleyzac N, Galambrun C, Mialou V, Bertrand Y,
Souillet G, et al. Influence of underlying disease on busulfan
disposition in pediatric bone marrow transplant recipients: a
nonparametric population pharmacokinetic study. Ther Drug Monit.
2007;29(2):177-84.

26. Nguyen L, Fuller D, Lennon S, Leger F, Puozzo C. I.V. busulfan in
pediatrics: a novel dosing to improve safety/efficacy for hematopoietic
progenitor cell transplantation recipients. Bone Marrow Transplant.
2004;33(10):979-87.

27. Aouri M, Barcelo C, Guidi M, Rotger M, Cavassini M, Hizrel C, et al.
Population Pharmacokinetics and Pharmacogenetics Analysis of Rilpivirine
in HIV-1-Infected Individuals. Antimicrobial agents and chemotherapy.
2017;61(1).

28. Ette EI, Ludden TM. Population pharmacokinetic modeling: the
importance of informative graphics. Pharmaceutical research.
1995;12(12):1845-55.

29. Karlsson MO, Savic RM. Diagnosing Model Diagnostics. Clinical
Pharmacology \& Therapeutics. 2007;82(1):17-20.

30. Hooker AC, Staatz CE, Karlsson MO. Conditional weighted residuals
(CWRES): a model diagnostic for the FOCE method. Pharmaceutical
research. 2007;24(12):2187-97.

31. Ette EI. Stability and performance of a population pharmacokinetic
model. Journal of clinical pharmacology. 1997;37(6):486-95.

32. van der Galiën R, Boveneind-Vrubleuskaya NV, Peloquin C, Skrahina A,
Touw DJ, Alffenaar JC. Pharmacokinetic Modeling, Simulation, and
Development of a Limited Sampling Strategy of Cycloserine in Patients
with Multidrug-/Extensively Drug-Resistant Tuberculosis. Clin
Pharmacokinet. 2020;59(7):899-910.

33. Nankivell BJ, Gruenewald SM, Allen RD, Chapman JR. Predicting
glomerular filtration rate after kidney transplantation.
Transplantation. 1995;59(12):1683-9.

34. Jones RJ, Lee KS, Beschorner WE, Vogel VG, Grochow LB, Braine HG, et
al. Venoocclusive disease of the liver following bone marrow
transplantation. Transplantation. 1987;44(6):778-83.

35. Harris AC, Young R, Devine S, Hogan WJ, Ayuk F, Bunworasate U, et
al. International, Multicenter Standardization of Acute
Graft-versus-Host Disease Clinical Data Collection: A Report from the
Mount Sinai Acute GVHD International Consortium. Biology of blood and
marrow transplantation : journal of the American Society for Blood and
Marrow Transplantation. 2016;22(1):4-10.

36. Kim MG, Kwak A, Choi B, Ji E, Oh JM, Kim K. Effect of glutathione
S-transferase genetic polymorphisms on busulfan pharmacokinetics and
veno-occlusive disease in hematopoietic stem cell transplantation: A
meta-analysis. Basic Clin Pharmacol Toxicol. 2019;124(6):691-703.

37. Ansari M, Curtis PH, Uppugunduri CRS, Rezgui MA, Nava T, Mlakar V,
et al. GSTA1 diplotypes affect busulfan clearance and toxicity in
children undergoing allogeneic hematopoietic stem cell transplantation:
a multicenter study. Oncotarget. 2017;8(53):90852-67.

38. Nava T, Kassir N, Rezgui MA, Uppugunduri CRS, Huezo-Diaz Curtis P,
Duval M, et al. Incorporation of GSTA1 genetic variations into a
population pharmacokinetic model for IV busulfan in paediatric
hematopoietic stem cell transplantation. Br J Clin Pharmacol.
2018;84(7):1494-504.

39. Nishikawa T, Yamaguchi H, Ikawa K, Nakayama K, Higashi E, Miyahara
E, et al. Influence of GST polymorphisms on busulfan pharmacokinetics in
Japanese children. Pediatrics international : official journal of the
Japan Pediatric Society. 2019;61(6):558-65.

40. Yee SW, Mefford JA, Singh N, Percival ME, Stecula A, Yang K, et al.
Impact of polymorphisms in drug pathway genes on disease-free survival
in adults with acute myeloid leukemia. Journal of human genetics.
2013;58(6):353-61.

41. Perkins JB, Kim J, Anasetti C, Fernandez HF, Perez LE, Ayala E, et
al. Maximally tolerated busulfan systemic exposure in combination with
fludarabine as conditioning before allogeneic hematopoietic cell
transplantation. Biology of blood and marrow transplantation : journal
of the American Society for Blood and Marrow Transplantation.
2012;18(7):1099-107.

42. Yeh RF, Pawlikowski MA, Blough DK, McDonald GB, O'Donnell PV,
Rezvani A, et al. Accurate targeting of daily intravenous busulfan with
8-hour blood sampling in hospitalized adult hematopoietic cell
transplant recipients. Biology of blood and marrow transplantation :
journal of the American Society for Blood and Marrow Transplantation.
2012;18(2):265-72.

43. de Lima M, Couriel D, Thall PF, Wang X, Madden T, Jones R, et al.
Once-daily intravenous busulfan and fludarabine: clinical and
pharmacokinetic results of a myeloablative, reduced-toxicity
conditioning regimen for allogeneic stem cell transplantation in AML and
MDS. Blood. 2004;104(3):857-64.

44. Russell JA, Tran HT, Quinlan D, Chaudhry A, Duggan P, Brown C, et
al. Once-daily intravenous busulfan given with fludarabine as
conditioning for allogeneic stem cell transplantation: study of
pharmacokinetics and early clinical outcomes. Biology of blood and
marrow transplantation : journal of the American Society for Blood and
Marrow Transplantation. 2002;8(9):468-76.

45. Gürlek Gökçebay D, Azik F, Ozbek N, Isik P, Avci Z, Tavil B, et al.
Clinical comparison of weight- and age-based strategy of dose
administration in children receiving intravenous busulfan for
hematopoietic stem cell transplantation. Pediatr Transplant.
2015;19(3):307-15.

46. McCune JS, Bemer MJ, Barrett JS, Scott Baker K, Gamis AS, Holford
NH. Busulfan in infant to adult hematopoietic cell transplant
recipients: a population pharmacokinetic model for initial and Bayesian
dose personalization. Clinical cancer research : an official journal of
the American Association for Cancer Research. 2014;20(3):754-63.

47. Vaughan WP, Carey D, Perry S, Westfall AO, Salzman DE. A limited
sampling strategy for pharmacokinetic directed therapy with intravenous
busulfan. Biology of blood and marrow transplantation : journal of the
American Society for Blood and Marrow Transplantation.
2002;8(11):619-24.

48. Teitelbaum Z, Nassar L, Scherb I, Fink D, Ring G, Lurie Y, et al.
Limited Sampling Strategies Supporting Individualized Dose Adjustment of
Intravenous Busulfan in Children and Young Adults. Ther Drug Monit.
2020;42(3):427-34.

49. Andersson BS, Thall PF, Madden T, Couriel D, Wang X, Tran HT, et al.
Busulfan systemic exposure relative to regimen-related toxicity and
acute graft-versus-host disease: defining a therapeutic window for i.v.
BuCy2 in chronic myelogenous leukemia. Biology of blood and marrow
transplantation : journal of the American Society for Blood and Marrow
Transplantation. 2002;8(9):477-85.

50. Lee JW, Kang HJ, Lee SH, Yu KS, Kim NH, Yuk YJ, et al. Highly
variable pharmacokinetics of once-daily intravenous busulfan when
combined with fludarabine in pediatric patients: phase I clinical study
for determination of optimal once-daily busulfan dose using
pharmacokinetic modeling. Biology of blood and marrow transplantation :
journal of the American Society for Blood and Marrow Transplantation.
2012;18(6):944-50.

51. McPherson ME, Hutcherson D, Olson E, Haight AE, Horan J, Chiang KY.
Safety and efficacy of targeted busulfan therapy in children undergoing
myeloablative matched sibling donor BMT for sickle cell disease. Bone
Marrow Transplant. 2011;46(1):27-33.

52. Davis JM, Ivanova A, Chung Y, Shaw JR, Rao KV, Ptachcinski JR, et
al. Evaluation of a Test Dose Strategy for Pharmacokinetically-Guided
Busulfan Dosing for Hematopoietic Stem Cell Transplantation. Biology of
blood and marrow transplantation : journal of the American Society for
Blood and Marrow Transplantation. 2019;25(2):391-7.

53. Srivastava A, Poonkuzhali B, Shaji RV, George B, Mathews V, Chandy
M, et al. Glutathione S-transferase M1 polymorphism: a risk factor for
hepatic venoocclusive disease in bone marrow transplantation. Blood.
2004;104(5):1574-7.

54. Krivoy N, Hoffer E, Lurie Y, Bentur Y, Rowe JM. Busulfan use in
hematopoietic stem cell transplantation: pharmacology, dose adjustment,
safety and efficacy in adults and children. Current drug safety.
2008;3(1):60-6.

55. Demirer T, Buckner CD, Appelbaum FR, Bensinger WI, Sanders J,
Lambert K, et al. Busulfan, cyclophosphamide and fractionated total body
irradiation for autologous or syngeneic marrow transplantation for acute
and chronic myelogenous leukemia: phase I dose escalation of busulfan
based on targeted plasma levels. Bone Marrow Transplant.
1996;17(4):491-5.

\textbf{Hosted file}

\verb`Table.docx`  available at \url{https://authorea.com/users/359718/articles/481548-optimization-of-busulfan-dosing-regimen-in-pediatric-patients-using-population-pharmacokinetic-model-incorporating-gst-polymorphisms}

\textbf{Hosted file}

\verb`Figures.docx`  available at \url{https://authorea.com/users/359718/articles/481548-optimization-of-busulfan-dosing-regimen-in-pediatric-patients-using-population-pharmacokinetic-model-incorporating-gst-polymorphisms}

\selectlanguage{english}
\FloatBarrier
\end{document}

