Model-Based Simulation of Maintenance Therapy of Childhood Acute Lymphoblastic Leukemia

Acute lymphoblastic leukemia is the most common malignancy in childhood. Successful treatment requires initial high-intensity chemotherapy, followed by low-intensity oral maintenance therapy with oral 6-mercaptopurine (6MP) and methotrexate (MTX) until 2–3 years after disease onset. However, intra- and inter-individual variability in the pharmacokinetics (PK) and pharmacodynamics (PD) of 6MP and MTX make it challenging to balance the desired antileukemic effects with undesired excessive myelosuppression during maintenance therapy. A model to simulate the dynamics of different cell types, especially neutrophils, would be a valuable contribution to improving treatment protocols (6MP and MTX dosing regimens) and a further step to understanding the heterogeneity in treatment efficacy and toxicity. We applied and modified a recently developed semi-mechanistic PK/PD model to neutrophils and analyzed their behavior using a non-linear mixed-effects modeling approach and clinical data obtained from 116 patients. The PK model of 6MP influenced the accuracy of absolute neutrophil count (ANC) predictions, whereas the PD effect of MTX did not. Predictions based on ANC were more accurate than those based on white blood cell counts. Using the new cross-validated mathematical model, simulations of different treatment protocols showed a linear dose-effect relationship and reduced ANC variability for constant dosages. Advanced modeling allows the identification of optimized control criteria and the weighting of specific influencing factors for protocol design and individually adapted therapy to exploit the optimal effect of maintenance therapy on survival.

Acute lymphoblastic leukemia is the most common malignancy in childhood. Successful treatment requires initial high-intensity chemotherapy, followed by low-intensity oral maintenance therapy with oral 6-mercaptopurine (6MP) and methotrexate (MTX) until 2-3 years after disease onset. However, intra-and inter-individual variability in the pharmacokinetics (PK) and pharmacodynamics (PD) of 6MP and MTX make it challenging to balance the desired antileukemic effects with undesired excessive myelosuppression during maintenance therapy. A model to simulate the dynamics of different cell types, especially neutrophils, would be a valuable contribution to improving treatment protocols (6MP and MTX dosing regimens) and a further step to understanding the heterogeneity in treatment efficacy and toxicity. We applied and modified a recently developed semi-mechanistic PK/PD model to neutrophils and analyzed their behavior using a non-linear mixed-effects modeling approach and clinical data obtained from 116 patients. The PK model of 6MP influenced the accuracy of absolute neutrophil count (ANC) predictions, whereas the PD effect of MTX did not. Predictions based on ANC were more accurate than those based on white blood cell counts. Using the new cross-validated mathematical model, simulations of different treatment protocols showed a linear dose-effect relationship and reduced ANC variability for constant dosages. Advanced modeling allows the identification of optimized control criteria and the weighting of specific influencing factors for protocol design and individually adapted therapy to exploit the optimal effect of maintenance therapy on survival.

INTRODUCTION
Acute lymphoblastic leukemia (ALL), characterized by malignant white blood cells (WBCs) and displacement of normal hematopoiesis, is the most common childhood malignancy (Hoffbrand et al., 2016). The treatment of childhood ALL is based on combination chemotherapy and begins with intensive, high-dose treatment for approximately 6 months (the so-called induction and consolidation therapy) followed by less-intensive, low-dose treatment [so-called maintenance therapy (MT)] that lasts for 2-3 years after disease onset. The goal of induction and consolidation therapy is to achieve remission via lymphoblast elimination below the limit of detection, and the high intensity of these therapy elements limits further therapy intensification using conventional chemotherapy. Subsequent MT is essential to prevent disease relapse, and aims to maintain prolonged antileukemic activity against residual lymphoblasts, with minimal adverse events. MT includes daily oral 6-mercaptopurine (6MP) and weekly oral methotrexate (MTX). Both drugs cause myelosuppression through their metabolized active forms (Schmiegelow et al., 2014). Blood count tests are performed regularly to ensure adequate WBC and absolute neutrophil count (ANC) suppression as a surrogate marker for antileukemic activity, without unintended excessive myelotoxicity. However, there exists no international consensus for MT dosing strategies and target levels for WBC and ANC suppression (i.e., what dose to start with, and when and how to increase or decrease chemotherapy). Empirical evaluation of different MT strategies using randomized clinical trials would be extremely challenging, due to the probably moderate effect size, the length of MT, the latency of clinically relevant endpoints, and the risk of compromising the current overall favorable outcome of childhood ALL. However, certain levels of WBC and ANC are established factors for survival, relapse or death, and other adverse events (e.g., infection), respectively. Therefore, a simulation model of childhood ALL MT could support the development of future MT strategies by identifying those strategies that achieve established survival factors best while avoiding established risk factors. Mathematical models describing the pharmacokinetics (PK) of 6MP and MTX and their pharmacodynamic (PD) effects on neutrophils may help clarify the drug-exposure relationship, predict the ANC dynamics, adapt subsequent dosing amounts, and stratify patients into groups with different drug responses. Several PK models for 6MP (Hawwa et al., 2008;Jayachandran et al., 2014Jayachandran et al., , 2015 and MTX (Godfrey et al., 1998;Panetta et al., 2002Panetta et al., , 2010Nagulu et al., 2010;Rühs et al., 2012;Korell et al., 2013;Hui et al., 2019) have been published, but not all have been developed with low-dosage treatments and validated in the pediatric population. To the best of our knowledge, there are only three publications (Jayachandran et al., 2014;Le et al., 2018;Karppinen et al., 2019) in which some of the PK models or their simplifications were linked to transient PD compartment models (Upton and Mould, 2014). The models were individually fitted to WBC counts and different prediction and optimization studies were conducted.
Here, we developed a population PK/PD model for maintenance treatment of ALL in children based on the approach used by Le et al. (2018) with a modified underlying PK model. As ANCs are the best established risk and survival factors, we adapted the model to predict ANCs instead of WBCs. The model was fitted to and validated on a dataset consisting of weekly ANC measurements obtained from 116 patients treated with daily oral 6MP and weekly oral MTX over an average of 459 (range, 200-581) days. We started our investigations with a PK/PD model considering 6MP and MTX but the constant administration ratio hampered the identification of separate PD effects. Further, the PK of MTX had no significant impact on the improvement of the model fitting, similar to the mathematical approach in Karppinen et al. (2019) and the clinical findings of NUDT15 genetics conferring 6MP but no MTX sensitivities (Tsujimoto et al., 2018). Thus, the final model only contains the PK of 6MP. We come back to this issue in the discussion. Then, for each patient, we simulated different therapy protocols (6MP dosing regimens), and compared the resulting predictions.

Data
The data used in this study were obtained retrospectively from 116 children who were diagnosed with de novo ALL at university hospitals in Erlangen and Dresden and treated according to the AIEOP-BFM 2000 and 2009 protocols. A subset of this data set (WBC counts from nine patients) was used and described similarly in a previous study (Le et al., 2018). Patients were eligible if they were diagnosed with precursor B-cell or T-cell ALL, negative for the BCR-ABL-and MLL-AF4 translocations, and started MT (i.e., did not experience relapse before the end of consolidation therapy and did not undergo stem cell transplantation). During MT administered according to the AIEOP-BFM 2000 and 2009 protocols, patients received oral chemotherapy with daily 6MP and once-weekly MTX until 2 years after ALL diagnosis. During MT, chemotherapy was applied to achieve antileukemic activity against lymphoblasts below the limit of detection. As a surrogate for antileukemic activity, WBC and ANC were measured regularly, with ANC <2 G/L, being correlated to a significantly better relapse-free survival (Schmiegelow et al., 2014), and ANC <0.5 G/L being an indicator of excessive myelosuppression. The target range for the WBC count was 1.5-3 G/L. The chemotherapeutic dose was reduced when cell counts fell below the lower limits (WBC count <1.5 G/L, ANC <0.5 G/L, lymphocyte count <0.3 G/L, and platelet count <0.05 G/L) or liver toxicity was suspected. For each patient included in the analysis, data regarding the following variables were recorded: gender, age, weight, height, body surface area (BSA), prescribed 6MP and MTX dosages (absolute and per BSA), WBC count, platelet count, lymphocyte and neutrophil counts, and therapy interruptions. In this study, we focused on 5897 ANCs and 6640 WBC counts, disregarding measurements of other cell types. We used both WBC counts and ANC separately and compared the accuracy of the resulting mathematical models. In all, 1150 ANC and 1289 WBC count measurements were excluded due to concurrent high C-reactive protein (CRP) levels indicating periods in which patients probably suffered from an infection. More precisely, we excluded measurements in the interval from 2 weeks before until 2 weeks after CRP levels of >5 mg/L were recorded. Among the remaining 4747 ANC measurements 56% were below the ANC threshold of 2 G/L, only 2% were below 0.5 G/L, and 54% were in the ANC target range 0.5-2 G/L. The demographic and clinical characteristics of the pediatric ALL population are shown in Table 1. The body surface area was calculated using the Mosteller formula.

Non-linear Mixed-Effects Modeling and Parameter Estimation
The non-linear mixed-effects (NLME) modeling (Bonate and Steimer, 2011) was based on the PK/PD model of Le et al. (2018). It describes the absorption of both drugs, 6MP and MTX, through the gastrointestinal (GI) tract into the plasma after oral administration and their metabolization to their active forms. The MTX metabolites MTXPG 2 to MTXGP 7 inhibit several enzymes responsible for DNA synthesis (Panetta et al., 2002). The active form of 6MP, 6-thioguanine nucleotides (6-TGNs), is incorporated into the DNA (Hawwa et al., 2008). Thus, both drugs negatively affect the hematopoiesis of neutrophils. During the model development, we replaced the 6MP PK model of Jayachandran et al. (2014) with the PK model described by Hawwa et al. (2008) to obtain a better response to 6MP dosage. The PK model of Jayachandran et al. (2014) was validated on concentration data of eight patients (adults) from Hindorf et al. (2006). However, the simulated 6-TGN concentrations coincided with data from pediatric patients reported by Hawwa et al. (2008); hence, it was a priori unclear which would give better results. Both compartment models have a comparable representation of the absorption and metabolic pathway of 6MP but the model of Hawwa et al. (2008) describes the metabolic transformations by first order kinetics instead of Michaelis-Menten kinetics. Further, the clearance is described by a BSA-dependent term, thus providing individualized PK profiles through patient characteristics. We also tested the influence of weekly MTX administration by either ignoring or considering the administrations and their resulting concentrations through the MTX PK model with a second PD parameter during model fitting. Additionally, we tested the myelosuppression model from Jayachandran et al. (2014), which contained a different feedback term for ANC recovery, but the accuracy decreased and this line of research was not further investigated. As a result, we identified one PK/PD model which described the clinical data best. This model was formulated as a system of ordinary differential equations (ODEs): the linear pharmacodynamic effect and the patient-specific bioavailable 6MP amount F u(t) of 6MP (implemented as point administration in NONMEM). The PK of 6MP is described by a three compartment model altered from Hawwa et al. (2008). A fraction of the orally administered 6MP dosage enters the GI tract where bioavailable 6MP is absorbed to the central compartment with the first order rate k a . In the central compartment, 6MP is eliminated by k 20 . The elimination also comprises metabolization of 6MP with the rate k me out of which a fraction FM 3 is metabolized to the active form 6-TGN. 6-TGN is then cleared by the BSA-dependent clearance term CL 6tgn . The hematopoiesis of neutrophils is described via a chain of five compartments with equivalent transition rates k tr representing the mean maturation time of the neutrophils (De Souza et al., 2018). The proliferation rate of hematopoietic stem cells k prol is equivalent to the transition rate k tr guaranteeing homeostasis (De Souza et al., 2018). Deviations from the neutrophil baseline Base are compensated by the feedback regulation (Base/x ma ) γ reflecting the granulocyte colony-stimulating factor (GCSF) controlled proliferation of neutrophils (Friberg et al., 2002;Quartino et al., 2012;Henrich et al., 2017;Jost et al., 2019). As the active forms of both drugs affect the proliferation process, the PD effect is modeled via a linear term with one joint parameter slope multiplied to the feedback-regulated first order proliferation rate constant. Other modeling approaches for the incorporation of the PD effect previously showed worst results in model fitting such that we focused on the described term which is additionally more plausible regarding the PD effect, i.g. an impaired proliferation through the incorporation of the metabolized drug into the DNA (Jost et al., 2019). Matured neutrophils die by the process of apoptosis with the rate k ma . A schematic representation of the model is shown in Figure 1 and model constants are listed in Table 2. As no PK biomarkers were measured in the examined dataset, we relied on published PK models and individualized the PD models with respect to individual sets of PD parameters. In the following, we describe the NLME parameter estimation approach. Therefore, we summarize model (1) for patient i aṡ with u i (t) the individual treatment schedule and θ i = (Base i , k tr,i , γ i , slope i ) T the patient specific parameter values of FIGURE 1 | Visualization of the final compartment model used for the population PK/PD analysis. The underlying mathematical models for the PK of 6MP and the myelosuppression were published by Hawwa et al. (2008), respectively (Le et al., 2018). The PK of orally administered 6MP is described as a three compartment model. A fraction of the 6MP dosage (6MP dose multiplied with the bioavailability factor F) enters the gastrointestinal (GI) tract where bioavailable 6MP is absorbed into the central compartment by the rate k a . In the central compartment, 6MP is eliminated with the first order kinetics k 20 . The elimination rate also comprises metabolization of 6MP to its active form 6-thioguanine nucleotide (6-TGN) with the rate FM 3 k m . The hematopoiesis of neutrophils is described by a chain of five compartments. The first compartment represents the hematopoietic stem cells proliferating with the rate k prol . The maturation process with equivalent transition rates k tr is represented by three intermediate compartments after which matured cells enter the circulating blood (last compartment). Matured cells die by the process of apoptosis with the rate k ma . The neutrophil baseline Base is maintained by the feedback term (B/x ma ) γ . As 6-TGN is incorporated into the DNA leading to cell apoptosis, the proliferation process is negatively affected by a linear PD function E.
the steady state of neutrophils Base, the transition rate k tr , the feedback term γ , and the PD effect slope. The vector θ i contains the fixed effect parameters Base, k tr , γ and slope for all patients and the individual realizations η i ∈ R 4 , i = 1, . . . , 116 of the random variable with the mean 0 ∈ R 4 and the diagonal variance matrix ∈ R 4×4 with the diagonal vector ω 2 = (ω 2 1 , ω 2 2 , ω 2 3 , ω 2 4 ) T . Interindividual variability (IIV) was assumed as log-normally distributed for all four parameters resulting in the following relation between fixed and random effects: summarized as θ i = g(θ , η i ) and for the description of the residual variability a proportional error model was used with a normally distributed measurement error ε ∼ N (0, σ 2 ) and ANC count measurements y ij .
The parameters were estimated using the first order conditional estimation (FOCE) method with η-ε interaction. This approximation method results in the parameter estimation method consisting of two nested optimization problems and t f i being the time point of patient i's last ANC measurement. The two parameter estimation problems (estimating θ , ω 2 , σ 2 with fixed η i and vice versa) are iteratively solved until a convergence criterion is fulfilled (Bae and Yim, 2016). For the detailed derivation of the FOCE method with η-ε interaction and the formulation of the resulting objective functions L FOCEi i,outer and L FOCEi MAP we refer the interested reader to Wang (2007) and Demidenko (2013) as we confine our analysis on the application of the parameter estimation method.

Out-of-Sample Validation
The reliability of the final population PK/PD model was tested via out-of-sample cross-validation. For each patient, the first 70% of ANC measurements were used for parameter estimation and the final 30% were used to evaluate the model predictions. Model accuracy and predictability were evaluated using the root mean squared error (RMSE) and the mean absolute error (MAE).

Simulation Study
We compared individual simulated minimal, median, and maximal ANCs resulting from the application of different dosing regimens (MT dosage over time). The choice of the different doses described in Table 3 was based on ALL treatment protocols (AIEOP-BFM 2009 with EudraCT number 2007-004270-43, NOPHO-ALL 2008-003235-20, andUKALL 2010-020924-22). In particular, we sought to investigate the relationship between an increased total amount of chemotherapy (higher dosage) and plausibly reduced ANC in the in silico simulations. Throughout, we used the fitted models (estimated model parameters) from section 2.2 and only varied the chemotherapy dosage. The simulated ANC values were obtained from the individual actual measurement time points.

Software
The population PK/PD analysis was performed with the NLME modeling program NONMEM 7.4 (ICON Plc., Dublin, Irland) (Beal et al., 2009). There exist several other software packages for parameter estimation of NLME models providing the same or similar algorithms. A variety of algorithms are provided in R Core Team ( Table 4 shows RMSEs, MAEs, and final objective function values for four different parameter estimations. Here, we compared the usage of different PK/PD models and parameter estimations based on either WBC counts or ANCs. First, the explicit consideration of MTX within the PK/PD model of Le et al. (2018) only had a minimal/non-significant effect on the model accuracy, so we fixed it to the ratio 2.5:1 between 6MP and MTX and neglected the PK of MTX in the following. Second, our results showed that the use of the PK model of Hawwa et al. (2008) increased the sensitivity of the PD effect and the model accuracy compared to the 6MP PK model of Jayachandran et al. (2014). Third, ANC measurements resulted in higher accuracy than did WBC measurements.   Figure 3 shows the good agreement of model response and measurements for the median (solid line) and  Table 4). For each patient, the individual 6-mercaptopurine (6MP [mg]) dosing protocol is presented in rows 2, 4 and 6, indicating dose changes for efficacy adjustments. The daily oral 6MP administration, ranging from 10 to 60 mg for the three patients, are presented as filled areas and corresponds to the control function u(t) in model (1).

Parameter Estimation
Frontiers in Physiology | www.frontiersin.org  97.5th percentile (dashed line) with a slight underprediction of the model for low ANC values. The 95% confidence interval of the model simulation median was very thin, indicative of high prediction accuracy. The fixed effect estimate for the ANC steady state was slightly higher than the target range limit of 2 G/L. The estimated transition rate of 0.148 resulted in a mean maturation time (MMT = n tr /k tr ) of 487 h (20.3 days) (De Souza et al., 2018). The interindividual variability and residual error were within reasonable ranges. The goodness-of-fit plot in Supplemental Data shows the results of out-of-sample cross-validation. It reflected reasonable model accuracy for fitted (blue) and predicted (red) ANC measurements with spreading around the line of identity because the model was not able (and not intended) to hit the lower and upper peaks of the measurements. The values of estimated model parameters both for the in-sample and out-of-sample calculations are shown in Table 5. The parameter values for slope and Base were reduced and the value of γ was slightly increased for the estimates based on 70% of the ANC. The interindividual variability (IIV) for the slope was significantly larger whereas the IIV of k tr was smaller. To evaluate the model accuracy, we calculated the median and standard deviation of the individual MAEs and RMSEs, showing the expected decrease in accuracy for out-of-sample predictions.  Table 3. We want to stress three main observations. First, a comparison of the first two entries of the three boxplots confirmed an already known result. The personalized models could reproduce the clinical ANC data on average quite well, with the exception of extreme values quantitatively confirming the observation made in Figure 2. Given the similarity of simulated and observed median values, we continued with an objective comparison only of the simulated results (protocols 2-6).

Simulation
Second, a comparison of the protocols 3-6 (25, 50, 75, and 100 mg/m 2 BSA 6MP) showed a significant and linear dosageeffect relationship with respect to the total amount of 6MP administered, which is, of course, proportional to the daily dose. All (minimal, median and maximal) ANC values decreased linearly, when daily dosing was increased linearly.
Third, a comparison of protocol 2 (the simulation of the real treatment) and protocols 3 and 4 (which gave lower and upper bounds on the total amount of administered 6MP in protocol 2, respectively) showed that the median ANC value of protocol 2 was indeed bounded by the two other values, however, for significantly lower minimal and higher maximal ANC values. Figure 5 shows an exemplary comparison of protocols 2-6 for one patient, highlighting lower peak values and smaller drug-induced steady state values when the dosing is linearly increased from 25 mg/m 2 to 100 mg/m 2 . The actual dosage administered to the patient (blue) ranged between the 25 mg/m 2 and 50 mg/m 2 protocols and resulted in similar ANC dynamics. At approximately day 240, the actual dosing was stopped for a short period, inducing stronger ANC oscillations in the subsequent treatment period and revealing a significant impact of the dosing regimen on the ANCs. This observation is even stronger regarding the proliferating cells as well as cells in the first transit compartment. Similar plots for all 116 patients are provided in the Supplemental Data.

Mathematical Model
We developed and fitted a population PK/PD model to assess the ANC dynamics during 6MP/MTX treatment, get a better understanding of dose adjustments, and identify solutions to the challenges that arise throughout MT. During the model development process we also fitted the model to WBC measurements. The resulting MAEs and RMSEs were worse compared to the values resulting from ANC measurements. This is probably due to the fact that WBCs comprise different cell lineages, with additional physiological effects that are not accounted for in the mathematical model. In future studies, the current model might be extended to further cell lineages. The models brought forth by Quartino et al. (2012) and Fornari et al.
(2019) might serve as a basis and drive the modeling process from a semi-mechanistic approach toward a more mechanistic one. In addition to using a population estimation approach and applying it to ANC instead of WBC, two modifications brought forth by Le et al. (2018) were shown to yield better results. First, the 6MP PK model of Jayachandran et al. (2014) was replaced by that of Hawwa et al. (2008). The first order kinetics in the PK model of Hawwa et al. (2008) compared to the Michaelis-Menten terms in the PK model of Jayachandran et al. (2014) resulted in more significant concentration changes with altered drug amounts consequently in a more sensitive PD effect. Second, the MTX PK model was completely omitted as the constant ratio of administered 6MP and MTX prevents a differentiation of separate PD effects. Further studies with measurements of drug concentrations, metabolites and clinical effects as cell counts would push forward the development of a mathematical model additionally including the PK of MTX to provide two distinct PD effects and to account for varying ratios of 6MP to MTX. For the currently available data, our new model, which indirectly agglomerates the effects of 6MP and MTX, appears to be a good choice (compare for Table 4).

Model Parameter Estimates
Looking at the resulting model parameter estimates listed in Table 5, the question arises as to how these values relate to known biological properties of hematopoiesis and myelosuppression and to other values from the literature. The estimated ANC steady state value Base was below the normal ANC range for children, but still higher than the desired ANC range of 0.5-2 G/L. Without treatment, the model-based ANCs would increase to normal patient-specific steady states. Thus, low ANC values were induced via MT or some of the aforementioned external events.
The estimated fixed-effects parameter value of the transition rate k tr = 0.148 was comparable with the published mean value (k tr = 0.1431) obtained from eight pediatric ALL patients from Riley Hospital for Children in Indianapolis (Jayachandran et al., 2014). For better interpretability, the transition rate parameter k tr can be transformed to the MMT (n tr /k tr ) of the neutrophils. The estimated MMT in our study, as well as the MMT from Jayachandran et al., are extremely high and do not coincide with biological findings of 3.9 days obtained by Hearn et al. (1998). This mismatch is a large disadvantage of the model as it fails to comply with biological properties, leading to falsely characterized physiological mechanisms and thus reduced model reliability. Jayachandran et al. did not discuss this issue, but a similar observation was made by Craig et al. (2016) who determined an estimated proliferation time of 26 days (Craig et al., 2016). In their work, the authors further presented model modifications to obtain a more realistic maturation time of 3.9 days. For this value we performed two parameter estimations with either Base as a parameter or fixed to 5 resulting in promising dynamics but worse RMSEs and MAEs. In future studies, the falsely determined MMT and possible model limitations for continuous low-dose treatments should be further investigated.
The feedback parameter (γ ) is significantly higher compared with published values (Friberg et al., 2002), indicating a stronger feedback mechanism during the daily chemotherapy over a long FIGURE 4 | Boxplots of minimal, median, and maximal (from left to right) individual ANCs for all 116 patients. Shown are values for the 6 different protocols from Table 3, observed for the first column and simulated for protocols 2-6. The target range (0.5-2.0 G/L) of the NOPHO/UK treatment protocol is shown as the gray background. Horizontal lines within the boxes are the medians, the upper and lower box limits are the first and third quartiles of the data, respectively. The whiskers indicate an even larger confidence region of these quartiles plus/minus 1.5-times the interquartile range. Beyond the whiskers, data are considered as outliers and are plotted as individual points. For the columns representing 25 mg/m 2 to 100 mg/m 2 , the total amount of 6MP administered is increasing. The median average individual daily doses actually administered for protocols 1 and 2 were 43.15 ± 10.5 mg/m 2 .
period. This is the first time estimated slope values of the linear PD function from the PK model of Hawwa et al. (2008) are presented; thus there are no available comparisons.

Simulation Results
The newly developed mathematical model enables us to perform a virtual comparison of different treatment protocols. The boxplots in Figure 4 show several interesting results.
First, the median and standard deviation of actual ANC measurements were very accurately matched by the simulation using the estimated parameters (compare the first two entries in the middle boxplot of Figure 4). Concerning the patientwise observed and simulated minimal and maximal ANC values, the model demonstrates a corresponding weakened chemotherapyinduced myelosuppression, respectively overproduction of ANCs compared to the high measured variability. This variability is biologically and clinically very plausible due to the aforementioned external events and uncertainties, although periods of severe infections were already excluded. The reproducibility of the median and avoidance of over-fitting of the extreme values are in our opinion good properties of a mathematical model. Given this good correspondence between cross-validated data and simulations, we felt encouraged to compare simulations of different treatment protocols as specified in Table 3. Note, however, that generalizations of mathematical models personalized for data from one protocol to another have to be considered with extreme care (compare the discussion for acute myeloid leukemia models by Jost et al., 2019). Further, we want to highlight that the current model is not intended to describe the ANC extrema such that the results of the simulation  Table 3 and an exemplary patient. Colors of the trajectories are identical to those used in Figure 4. The linear increase in dosing from 25 to 100 mg/m 2 forces the neutrophils (ANC) to lower peak values and a smaller drug-induced steady state value at the end of treatment. The actual dosage administered to the patient (blue) ranged between the 25 mg/m 2 and 50 mg/m 2 protocols and resulted in similar ANC dynamics. At approximately day 240, the actual dosing was stopped for a short period, inducing stronger ANC oscillations in the subsequent treatment period. This observation is even stronger regarding the proliferating cells as well as cells in the first transit compartment. Interestingly, these oscillations also continued for some time after the end of treatment. Due to the long simulation horizon, the 6-TGN dynamics are squeezed such that the concentrations between two administrations are not visible. study have to be treated with caution. The results shall serve as a preliminary assessment of the dose-effect relationship which has to be confirmed in future studies. The relationship might be stronger compared to the current model predictions and demonstrated by the clinical data in Figure 5. The impact of model variations on the outcome of simulation studies is usually significant. We tested the value of fixing the k tr parameter to represent a biologically plausible MMT of 3.9 days. This decreased the model accuracy (which is why the results are not included here), but still led qualitatively to the same subsequent effects.
Second, an approximately linear decrease in minimal, median and maximal values could be observed as the dosage increased linearly from 25 to 100 mg/m 2 with a slightly reduced decrease of the maximal ANC values. Again, this linear dose-effect relationship seems biologically plausible. For most of the simulations such as those shown in Figure 5, the maximal ANC value decreased. However, for other simulations (see Supplementary Material) stronger myelosuppression led to identical maximal ANC values. This effect is due to a feedback mechanism that may lead to increased proliferation for reduced ANC which leads to larger ANC values after some delay.
Third, a tendency for higher oscillations for treatments with pauses and changes in dosage was seen in a comparison of the simulated actual treatment protocol 2 and the constant administrations of protocols 3 and 4, which used lower/higher total amounts of 6MP. Again, an example of this can be seen in Figure 5. We believe that in the future adapted dosing schedules might take advantage of the chemotherapy-induced oscillations for an optimized dosing regimen. In the consolidation therapy of acute myeloid leukemia it was shown in silico that the timing of the treatment start can have a beneficial influence on the reduction of myelosuppression (Jost et al., 2019). However, high dose chemotherapy administered every 3 to 4 weeks provokes stronger periodic oscillations compared to the daily oral dosing which makes it more challenging to identify and capture the oscillations. For high dosage, previously a multi-compartment hematopoietic model was analayzed regarding Hopf bifurcation and an explicit analytical expression for the bifurcation point was provided depending on model parameters (Knauer et al., 2019). Oscillations of various blood cell populations have been observed in clinical data and partly investigated for different hematological disorders Colijn et al., 2006). The exact mechanisms and interaction between (1) stem cell cycling, (2) hematological disorder, and (3) drug exposure are still not fully understood. In our case, for all 116 patients in silico simulations showed that the oscillations were damped (in 84 cases into a steady state) once the chemotherapy was stopped, albeit with long time ranges of up to one year (see Supplemental Data for examples). Therefore, we assume that the oscillations in the ANCs observed in our simulations could be attributed to the influence of chemotherapy on the nonlinear dynamics of hematopoiesis. The connection between model-intrinsic and chemotherapy-induced oscillations should be assessed in detail in future studies. A stability analysis (Edelstein-Keshet, 2005) of the steady state could be performed (e.g., similar to Stiehl, 2014;Tetschke et al., 2018) to assess the theoretical properties of the model and relate them to the physiological behavior of neutrophils.

CONCLUSION
We presented a novel NLME model describing myelosuppression for ALL MT among children who received 6MP and MTX and was cross-validated on a data set of 4747 ANC measurements obtained from 116 patients. A comparison with alternative modeling approaches and using WBC counts instead of ANCs showed the benefit of this model. We could show a linear doseeffect relationship superimposed with fluctuations of varying magnitude. Mathematical simulations and more mechanistic modeling approaches will allow to improve the understanding of intrinsic and extrinsic influence factors on the aberrant hematopoiesis and chemotherapy-induced myelosuppression of pediatric ALL patients. Therefore, the monitoring of individual PK profiles and a subsequent analysis of the PK/PD relationship are mandatory next steps for a better dose-effect correlation.
In the future, based on the conduction of further PK and PD experiments driving the development of more advanced mathematical models together with the individual determination of response-related genotyping (Tsujimoto et al., 2018), MT protocols might be developed in silico, leading to individualized treatment protocols with better clinical outcomes.

DATA AVAILABILITY STATEMENT
The datasets generated from the parameter estimation process can be found in the Supplemental Material.

AUTHOR CONTRIBUTIONS
FJ developed the population PK/PD models, performed the numerical computations, and wrote the first draft of the manuscript. JZ, TL, TR, MR, MM, and SS contributed to the model development, the study designs, and the interpretation of the results. MSu and MSt provided clinical data. All authors contributed to writing the final manuscript.

FUNDING
This project has received funding from the European Research Council (ERC, grant agreement No 647573) and from the European Regional Development Fund (grants SynMODEST and SynIsItFlutter) under the European Union's Horizon 2020 research and innovation program.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys. 2020.00217/full#supplementary-material Goodness-of-fit plot presenting observed vs. individually calculated absolute neutrophil counts for 116 patients (see Supplemental Data 1).
Simulations of protocols 2-6 for all 116 patients, similar to Simulations of protocols 2 for all 116 patients to analyze the steady state behavior (see Supplemental Data 2).
The individual parameter estimates of the final PK/PD model are recorded in the file finalParameterEstimates191221.csv.
The individual 6-mercaptopurine (6MP) treatment schedules, the observed and calculated absolute neutrophil counts (ANCs) and the patient's weight, height and body surface area are recorded in the file NONMEMresultsALL191221.csv. The dataset contains the columns ID, TIME, DV, CMT, AMT, MDV, EVID, IPRED, WEIGHT, HEIGHT and BSA. ID serves as an identifier for the appropriate patient. TIME [day] either specifies the measurement times of ANCs or the time of oral 6MP administrations. DV [G/L] is the dependent variable, in our case the individual ANC measurements. The column CMT specifies the compartment in which a dosing or observation event occurs. AMT [mg] defines the amount of oral 6MP administration. The column MDV allows the user to inform NONMEM whether or not the value in the DV field is missing, but in our case the datasets do not contain missing measurements. The column EVID explicitly declares to NONMEM the type of the current record. EVID=0 defines the record as an observation event and EVID=1 defines the record as a dose event. The column IPRED contains the calculated ANCs derived by the PK/PD model after parameter estimation.