Population pharmacokinetics of treosulfan and development of a limited sampling strategy in children prior to hematopoietic stem cell transplantation

Purpose There is an increasing interest in use of treosulfan (TREO), a structural analogue of busulfan, as an agent in conditioning regimens prior to hematopoietic stem cell transplantation (HSCT), both in pediatric and adult populations. The aim of this study was to develop a population pharmacokinetic model and to establish limited sampling strategies (LSSs) enabling accurate estimation of exposure to this drug. Methods The study included 15 pediatric patients with malignant and non-malignant diseases, undergoing conditioning regimens prior to HSCT including TREO administered as a 1 h or 2 h infusion at daily doses of 10, 12, or 14 g/m2. A population pharmacokinetic model was developed by means of non-linear mixed-effect modeling approach in Monolix® software. Multivariate regression analysis and Bayesian method were used to develop 2- and 3-point strategies for estimation of exposure to TREO. Results Pharmacokinetics of TREO was best described with a two-compartmental linear model with proportional residual error. Following sampling schedules allowed accurate estimation of exposure to TREO: 1 h and 6 h or 1 h, 2 h, and 6 h for a TREO dose 12 g/m2 in a 1 h infusion, or at 2 h and 6 h or 2 h, 4 h, and 8 h for a TREO dose of 12 g/m2 and 14 g/m2 in a 2 h infusion. Conclusions A two-compartmental population pharmacokinetic model of TREO was developed and successfully used to establish 2- and 3-point LSSs for accurate and precise estimation of TREO AUC0→∞.


Introduction
Allogeneic hematopoietic stem cell transplantation (allo-HSCT) is a procedure aimed at reconstituting normal hematopoiesis in malignancies and non-malignant hematopoietic disorders [1]. A conditioning regimen prior to HSCT is required. It should demonstrate myeloablative and immunosuppressive and, in case of malignant disorders, anti-malignancy properties to prevent graft rejection, graft versus host disease and post-transplant relapse of malignancy. Commonly, the myeloablative conditioning procedure consists of fractionated total body irradiation or myeloablative dose of busulfan combined with high-dose cytostatics, such as etoposide, fludarabine, melphalan, thiotepa or cyclophosphamide [2]. Recent clinical studies indicate that also high-dose treosulfan (TREO), which is a structural analogue of busulfan, demonstrates significant myeloablative and immunosuppressive properties as well as anti-malignant activity in case of hematological malignancies and some solid tumors [3][4][5]. According to these studies, TREO has relatively mild toxicity profile [5,6].
TREO is a prodrug and undergoes a non-enzymatic pHd e p e n d e n t r e a c t i o n t o t h e m o n o e p o x y -a n d diepoxytransformers ( Fig. 1) [7]. The products of this reaction alkylate DNA at the N7 position of guanine [8]. The optimal dosing regimen of high-dose TREO is still not established and may vary between medical centers [5]. Most commonly, TREO is administered on three subsequent days prior to the transplant procedure (total dose of TREO 30-42 g/m 2 ) [4,6]. A daily dose of TREO (10-14 g/m 2 ) is administered in a single infusion.
Since the formation of epoxybutane derivatives from TREO is a non-enzymatic process, it might be assumed that measurement of the prodrug is adequate to describe the alkylating activity [9]. Until now only one limited sampling strategy (LSS) for estimating the exposure to TREO was proposed [10]. However the underlying population pharmacokinetic model assumed a one-compartment linear pharmacokinetics, while other authors suggest that two-compartmental model might more accurately describe changes of TREO concentration over time [9,[11][12][13]. Therefore the aim of this study was to develop a population pharmacokinetic model and to establish LSSs enabling accurate estimation of exposure to TREO.

Patients' characteristics
The study included 15 Table 1. Conditioning regimens prior to HSCT included TREO administered as a 1 h or 2 h infusion at daily doses of 10, 12, or 14 g/m 2 . The body surface area was calculated at clinical sites by means of Mosteller method [14]. The study protocol was approved by the local Ethical Committee at the Poznan University of Medical Sciences and is in accordance with the 1964 Declaration of Helsinki and its later amendments. Informed consent was obtained from the parents prior to initiating the study.

Sampling protocol and determination of TREO
The samples were drawn in the first day of the therapy from all of the patients. Two different sampling protocols were applied. From 7 patients the full blood samples were drawn at 0.5, 1, 3, 4, 6 and 8 h after the beginning of infusion, while from the remaining 8 patients a more dense sampling was allowed, at 0.5, 1, 1.5, 2, 2.5, 3, 4, 5, 8 and 12 h after the start of infusion. Immediately after collection 50 μl of 1 M citric acid per 1 ml of full blood was added, to avoid ex vivo transformation of TREO to its epoxides. Subsequently, the samples were centrifuged and the obtained plasma was stored at − 20°C until the analysis.
Concentrations of TREO were determined by a validated HPLC-MS/MS method. The method validation, as well as preliminary pharmacokinetic analysis was published in details elsewhere [12,15]. Briefly, the applied method allowed accurate determination of TREO in the plasma samples in ranges 0.2-5720 μM (56.6 ng/ml-1.59 mg/ml). The lower limit of quantitation was 0.2 μM. The inter-day and intra-day precision and accuracy were calculated according to regulations of the European Medicines Agency for bioanalytical method validation. The precision of the method, described by the coefficient of variation was 1.8-11.5%, while method accuracy described with the relative error was 0.02-11.8%.

Methods and software
The model development was performed in Monolix 2016R1 software (Lixoft SAS, Antony, France, http://lixoft.com/ Fig. 1 Metabolic activation of treosulfan to its active mono-and diepoxide products/monolix/) by means of stochastic approximation of the standard expectation maximization (SAEM) algorithm for non-linear mixed-effects models without approximation. The maximum number of iterations at each stage of population parameter estimation was automatically determined by the algorithm (K 1 = 'auto', K 2 = 'auto'). The quality of SAEM algorithm convergence was inspected at each model estimation step. Minimum 4 Markov chains were set at estimation. Conditional means and standard deviations of individual pharmacokinetic parameters were estimated with a Markov Chain Monte Carlo method (MCMC). Improvement in the model fit was evaluated with the likelihood ratio test. The difference in the minimum objective function value (MOFV) of 10.8 (p < 0. 001) between nested models was considered significant. Also, Akaike information criterion (AIC) and Bayesian Information criterion (BIC) were calculated and models with lower values of AIC and BIC were considered as better fitted to the observed data. Calculation of MOFV was performed by linearization in the initial stages of decision-making, while an importance sampling method was used for the final model selection. Visual examination of goodness-of-fit was based on the following plots: individual (IPRED) and population predicted (PPRED) concentrations versus observed concentrations, individual fits, population weighted (PWRES) and individual weighted residuals (IWRES) versus time and predicted concentrations, normalized prediction distribution errors (NPDE) versus time and predicted concentrations, histograms and quantile-quantile plots.

Structural and error model selection
One-, two-and three-compartmental models for intravenous infusion with first-order elimination were examined. Lognormal distribution of pharmacokinetic parameters was assumed and interindividual variability elements (IIV) were described with an exponential model as follows (Eq. 1): where θ ij is a value of j-th pharmacokinetic parameter for i-th individual, θ j is the population parameter estimate and η ij is a random variable characterizing IIV. Covariance between IIV elements was inspected after building a structural model. Diagonal, full and partial covariance matrices were examined and significant values were retained in the model. Standard errors of model parameters were calculated by means of the linearization algorithm; however, the values obtained from the final model were obtained by stochastic approximation, which is more time consuming, but at the same time a more precise method [16]. Also, ηshrinkage was calculated with a following equation (Eq. 2): whereη is the posterior estimate for individuals based on empirical Bayes estimates andω is the estimated standard deviation for the corresponding random effect [16]. According to Savic and Karlsson [17], high shrinkage (above 20-30%) is associated with insufficient informativeness of diagnostics based on empirical Bayes estimates, which include IPRED and IWRES. Additive, proportional and combined (additive and proportional) error models describing residual unexplained variability (RV) were examined and following equation was applied (Eq. 3): where C obs and C pred are observed and predicted concentrations of treosulfan, ε 1 is a variable associated with proportional RV and ε 2 defines additive portion of RV.

Covariate selection
Covariate model was established in a forward inclusionbackward elimination manner. Weight and sex were examined as potentially influential covariates. Creatinine clearance was not considered in the covariate analysis due to lack of the data for 7 patients out of total 15 subjects. Weight, a continuous covariate, was standardized to adult bodyweight (70 kg). Also, allometric scaling factor was included for clearance parameters, as suggested by Anderson and Holford [18]. Wald test implemented in Monolix was used to evaluate the covariates and p value < 0.05 was considered significant. After including all the significant covariates, a full model was created. Next, the covariates were eliminated from the model in a step-by-step manner. A covariate was retained if the MOFV increased by more than 6.67 (p < 0.01). As a result, the final model was obtained and applied in the further analysis.

Model evaluation
Due to low number of patients included in the study, internal methods of model validation were applied. Predictioncorrected visual predictive check (pcVPC) was performed for 1000 simulated observations. Contrary to a standard VPC, pcVPC allows accurate visual presentation of the simulation results without losing power due to data stratification [19]. In this approach observed and simulated dependent variables are normalized basing on a typical population prediction for the median independent variable in each bin [20]. 5th, median and 95th percentiles of the simulated data were plotted and compared to the corresponding percentiles of the observed data. The other recommended validation procedure is bootstrapping. Due to the fact that bootstrapping algorithms are not implemented in Monolix, Wings for NONMEM (WFN, version 742, http://wfn.sourceforge.net/, Nick Holford, University of Auckland, New Zealand) with the mlxbs script was applied. 1000 bootstrapped datasets were used to determine medians of each estimated parameter, as well as 5th and 95th confidence intervals (CI). Calculated values were compared with those obtained from the final model.

Data simulation
A simulation-based approach was applied to develop LSSs for estimation of exposure to TREO in pediatric patients. Firstly, the distribution of significant covariates included in the final model was visually inspected.
Deviations from normal distribution were tested with Shapiro-Wilk's test. Next, a group of 100 virtual patients, with distribution of significant covariates resembling the In the first step, areas under time-concentration curves from time 0 to infinity (AUC 0 → ∞ ) were calculated by means of a non-compartmental method analysis (NCA). The AUC 0 → ∞ was calculated with a linear-up log-down interpolation method by means of PKNCA package (version 0.8.1, https://cran.rproject.org/package=PKNCA) run through the R software. Next, each of the three simulated patient groups were randomly divided into two approximately equal subgroups. First subgroup was a Blearning^one, while the second subgroup, which was labeled Bvalidation,^was used to evaluate the predictive performance of the strategies. Multiple regression equations were calculated for models in which two or three samples were required to predict AUC 0 → ∞ . The calculations were performed in R software by means of leaps package (version 3.0, https://cran.r-project. org/package=leaps). The models were stratified upon the adjusted coefficient of determination value (R 2 ). Models which assumed drawing during the infusion or sampling at12 h after the beginning of infusions were discarded. Finally, best two and three-point models for each type of TREO administration were chosen for further validation.

Bayesian estimation of exposure to TREO
The procedure was performed as described by Alsultan et al. [21]. In the analysis, the simulated datasets described in chapter 2.4.1 were used. Also, the best sampling strategies selected by linear regression fitting were applied. First, the estimates from the final model were fixed at the obtained values. Second, the datasets were created in which only the sampling times from the chosen LSSs were left. Then, the individual pharmacokinetic parameters were calculated for each subject. The exposure to TREO for each individual was calculated as following (Eq. 4): where D is the amount of TREO administered and Cl is total clearance.

Evaluation of the predictive performance of LSSs
The prediction performance of LSSs was evaluated as suggested by Sheiner and Beal [22]. Following parameters were calculated for each strategy: relative prediction error (PE), mean relative prediction error (MPE), mean absolute relative prediction error (MAPE) and root mean squared relative prediction error (RMSE). These parameters were estimated using following equations (Eq. 5-8): For each LSS, distribution symmetry and range of relative errors was verified [23].
For the linear regression method, the best strategies were used to predict AUC 0 → ∞ of patients from the Bvalidationŝ ubgroup. Next, also, the LSS equations were used to predict AUC 0 → ∞ calculated from the acquired experimental data, and the cumulative predictive performance was calculated for the best 2-and 3-point strategies as presented above. For the Bayesian method, the chosen strategies were used to predict AUC 0 → ∞ of all patients from a given dosing subgroup (100 patients in each).

Population pharmacokinetic analysis
A total of 110 samples were obtained from the patients and analyzed in the study. All concentrations of TREO were above the LOQ of the HPLC-MS/MS method. A spaghetti plot with individual time-concentration curves is presented in Fig. 2. It was found that the experimental data were best described by a linear two-compartmental model with a proportional error (63.76 decrease of MOFV), where V 1 and V 2 describe central and peripheral compartment volumes, respectively, Cl is clearance and Q is an intercompartmental clearance. Analysis of covariance matrix showed that there is a significant covariance between random variability of Cl and V 1 . The only covariate included in the final model was patient's weight, as described by following equation (Eq. 9): where θ ij is a value of j-th pharmacokinetic parameter for i-th individual, θ j is the population parameter estimate, BW i is a bodyweight of the i-th individual centered on a typical weight of 70 kg, while β is a scaling exponent and η ij is a random variable characterizing IIV. The values of the exponents were first estimated and equaled 0.804 for Cl, 0.959 for V 1 and 0.925 for V 2 . For subsequent analysis, these values were fixed to most widely used values of 1 for volume parameters or to 0.75 for clearance. Also, the IIV parameters were evaluated for CL, V 1 and Q only.
In the modeling process, the standard error for the estimated V 2 IIV was over 100%. Hence it was decided to remove the IIV on V 2 from the model. It was found that addition of bodyweight as a covariate decreased the IIV of CL (65.3 to 25.5%) and V 1 (84.5 to 51.4%). Interestingly, addition of bodyweight as a covariate for Q increased the MOFV and worsened the model fit. Therefore this relationship was not included in the model. The estimates derived from the final pharmacokinetic model are presented in Table 2. Calculated η-shrinkage was 1% for Cl, 19% for V 1 and 18% for Q. Figure 3 presents basic goodness-of-fit diagnostic plots of the final model. Observed concentrations (OBS) plotted vs. population predicted (PPRED) and individual predicted concentrations (IPRED) are scattered randomly around the line of identity and the spline is close to the identity line (Fig. 3A). Plots of IWRES and PWRES vs. time and PPRED are presented in Fig. 3B. For both types of residuals, the points are randomly scattered along the y = 0 line and no significant trends are visible. Moreover, most points fall within ±2 standard deviation with only a few points with deviations larger than 3. Simulation-based plots of NPDE vs. time and PPRED (Fig. 3C) show that the 5th, median and 95th percentiles of empirical data (solid lines) fall within the corresponding percentiles of the predictions (light and dark gray areas).
Performed pcVPC is presented in Fig. 4. The solid black lines, which represent 5th, median and 95th percentile of observed data fall within the areas representing respective prediction intervals of the simulated data. There are only 2 points which fall outside the plotted prediction intervals.
The results of bootstrapping are presented in Table 2. Calculated means as well as the 5th-95th confidence intervals are very similar to the estimates obtained with the SAEM algorithm for the primary dataset. The largest discrepancy was noted for the IIV of Q. However, the calculated standard error for this particular parameter was large (52.3%) and also small size of study group might account for the difference in the results.

Performance of selected LSSs
Covariate analysis showed that the patient's bodyweight is an important cofactor. Visual inspection of the data and results of Shapiro-Wilk's test (p = 0.182) it was assumed that the bodyweight was normally distributed. 100 patients with mean where BSA is body surface area in cm 2 and BW is bodyweight in g. Obtained body surface areas were used to calculate the exact dose of TREO (12 g/m 2 or 14 g/m 2 ) for the simulation analysis. For all three types of TREO administration the population pharmacokinetic estimates were not significantly different from the values obtained from the final model.
The best strategies and their performance for estimating exposure to TREO, based on the simulation study, are presented in Table 3. For 1-h infusion of 12 g/m 2 of TREO, the best strategies assumed sampling at 1 h and 6 h or at 1.5 h, 2 h and 6 h after the beginning of the infusion. For 2-h infusion of 12 g/m 2 an accurate prediction of AUC 0 → ∞ required determination of TREO concentration in samples drawn 2 h and 6 h or 2 h, 3 h and 8 h after the beginning of infusion. While for 2h infusion of 14 g/m 2 of the drug, sampling at 2 h and 6 h, or at 2 h, 4 h and 8 h was needed for the most accurate estimation of exposure. The R 2 is very close to the unity value and the prediction errors are overall very small. Noteworthy, the variability observed in the Bayesian prediction method is higher, as well as the prediction errors. Also, approximately 2% of the  with dots as observed treosulfan concentrations, bold lines as 5th, median and 95th percentile of observed concentrations, light gray area as 50% interval of simulated data and dark gray areas as 95% intervals of simulated data predictions based on this method were outside the 20% error boundary.
The performance of the proposed LSSs was also evaluated by comparison of observed and predicted AUC 0 → ∞ in the experimental group (Table 4). Bias in prediction of this parameter was within ±15% for all patients. Only 3-point strategies in which the exposure was estimated with Bayesian methods had RMSE slightly above 15%. Unfortunately, due to heterogeneity in the sampling designs and lack of samples in some time points, prediction of AUC 0 → ∞ was not possible for all of the patients, while for some individuals prediction was calculated only by means of 2-point strategies.

Discussion
There is an increasing interest in use of TREO as a basic agent in conditioning regimens prior to hematopoietic stem cell transplantation, both in pediatric and adult populations [3,25,26]. The aim of this study was to determine a population pharmacokinetic model for TREO basing on the data acquired for pediatric population and to develop LSSs for estimation of AUC of this drug.
According to Scheulen et al. [27], the pharmacokinetics of TREO is linear in a wide range of doses (20-56 g/m 2 ) therefore it was possible to pool patients with different dose levels and it might be assumed that the dose would not have an impact on the estimations. It was found that the pharmacokinetics of TREO was best described with a linear two-compartmental model with a proportional residual error. The covariate analysis showed that bodyweight was significantly associated with the estimated values of Cl, V 1 and V 2 . Interestingly, visual inspection of data and analysis of MOFV did not support inclusion of bodyweight as an important covariate of Q. Although the explanation of this phenomenon remains unknown, one of the reasons might be unique metabolic activation of TREO. Scaling of parameters, especially allometric scaling of clearance parameters, describes an increase in the metabolic rate [28]. At the same time activation of TREO is a pH-and temperature-dependent process and does not require enzymatic conversion. Therefore the lack of relationship between Q and bodyweight might, to some extent, reflect this pathway. The conversion of TREO to epoxides might occur throughout the whole body both in the central and peripheral compartments. However, more studies are needed, especially the combined parent-metabolite modeling approach, to sufficiently explore this observation. The other tested covariate, patient's sex, was not found significant. One explanation of these observations might be the fact that the study group included three girls only and this number might be insufficient to observe potentially existing relationships between estimated parameters. Also, some authors note that patient's sex might be an important factor in children older than 12 years of age LSS limited sampling strategy, R 2 adjusted coefficient of determination, PE relative prediction error, MPE mean relative prediction error, MAPE mean absolute relative prediction error, RMSE root mean squared relative prediction error, AUC pred predicted area under time-concentration curve, C nh concentration of treosulfan measured n hours after the beginning of infusion [29]. In the study group only four children (three boys and one girl) could be qualified into this category and therefore a more thorough investigation of this potential covariance could not be performed. Previous studies on pharmacokinetics of TREO indicate that up to 39% of total dose of this drug is excreted renally [12]. As a consequence, parameters which describe renal function such as creatinine clearance might be important covariates for Cl. However, in the present study, this factor was not included in the analysis. First of all, data on creatinine clearance was available for eight children (53% of total). Several methods to overcome the problem of missing data were considered in this analysis-exclusion of patients with lack of information on creatinine clearance or imputation with a median value [30]. Although removing individuals with missing covariate data gives relatively unbiased results, the estimates are less precise. Imputation of a median value might significantly bias the results of the study. Noteworthy, no significant trends were observed in the plots of neither Cl vs. creatinine clearance nor IIV Cl vs. creatinine clearance. Also, the available data on creatinine clearance indicate that the renal function in patients could be described as normal. Therefore, to investigate the influence of this parameter on pharmacokinetics of TREO a larger study group would be necessary. The parameter estimations obtained in the present study are similar to the results reported by other authors. In adult populations [9,13,26] estimated total clearance ranged from 8.7 l/ h to 13.5 l/h, while mean steady-state volume of distribution (V ss ), which represents both central and peripheral compartments, was ranging from 19.5 l to 34 l. The obtained values are also similar to the ones presented by ten Brink et al. (V c / 70 kg = 43.05 l, Cl/70 kg = 16.12 l/h) [10]; however, it has to be noted that the model assumed by those authors is onecompartmental. The IIVof calculated parameters was relatively large, reaching up to 51.4%. According to the presented goodness-of-fit plots (Fig. 3), pcVPC (Fig. 4) and bootstrap analysis (Table 2), proposed model adequately describes changes of TREO concentration over time.
Overall, the results are in some contrast with the results of ten Brink et al. [10] who developed a one-compartment model. These authors note that using a two-compartmental model did not sufficiently lower the MOFV (< 10.8 decrease in the MOFV). However, a systematic bias is noticeable in the conditional-weighted residuals vs. time graphs presented by these authors. It indicates that the addition of peripheral compartment might had improved the model. Therefore, the model developed in the present study might be superior to the previous one and allows better prediction of changes in TREO concentration over time.
The best LSSs developed on a basis of the proposed population pharmacokinetic model were shown to have good predictive properties (Table 3). Interestingly, the equations proposed for estimation of AUC 0 → ∞ on a 2-sample basis were very similar for 2-h infusions of 12 g/m 2 and 14 g/m 2 . This is an understandable consequence of the assumed linearity of TREO pharmacokinetics. Moreover, the strategies were applied to predict the AUC 0 → ∞ in the experimental group and the predictive performance was acceptable.
As administration of TREO in HSCT in pediatric populations is still a relatively new issue, because standard procedures rather include busulfan, a uniform system for estimation of exposure is yet to be evaluated. Also, the exact and detailed protocols for monitoring of TREO have also not been published yet, beside the work by ten Brink et al. [10]. TREO itself is a prodrug, and its pharmacological effect is dependent on the epoxytransformers. Since the concentrations of TREO are much higher than the concentrations of its Bmetabolites,t hey are much easier to monitor. Therefore the monitoring procedure based on the prodrug concentration was proposed. To our knowledge, no studies were published which would bind directly the efficacy and safety of TREO with its pharmacokinetics and several authors point out the need of further studies in this particular area [4,5]. It might be possible that the levels of TREO epoxides have higher predictive value for safety, toxicity and efficacy of this treatment. However, exact data which combine follow-up data and the epoxides concentrations would be required for this assessment. LSS limited sampling strategy, PE relative prediction error, MPE mean relative prediction error, MAPE mean absolute relative prediction error, RMSE root mean squared relative prediction error The major limitation of the study is small sample size. Therefore, an external validation procedure is required prior to employing the proposed LSS in the clinical practice.

Conclusions
In conclusion, in the present study a two-compartmental population pharmacokinetic model of TREO was developed and successfully used to establish 2-and 3-point LSSs for accurate and precise estimation of TREO AUC 0 → ∞ .