A mathematical model of white blood cell dynamics during maintenance therapy of childhood acute lymphoblastic leukemia

Abstract

Acute lymphoblastic leukemia is the most common malignancy in childhood and requires prolonged oral maintenance chemotherapy to prevent disease relapse after remission induction with intensive intravenous chemotherapy. In maintenance therapy, drug doses of 6-mercaptopurine (6-MP) and methotrexate (MTX) are adjusted to achieve sustained antileukemic activity without excessive myelosuppression. However, uncertainty exists regarding timing and extent of drug dose responses and optimal dose adaptation strategies. We propose a novel comprehensive mathematical model for 6-MP and MTX pharmacokinetics, pharmacodynamics and myelosuppression in acute lymphoblastic maintenance therapy. We personalize and cross-validate the mathematical model using clinical data and propose a real-time algorithm to predict chemotherapy responses with a clinical decision support system as a potential future application.

childhood acute lymphoblastic leukemia, chemotherapy, maintenance therapy, ordinary differential equations, parameter estimation, real-time algorithm

1. Introduction

Acute lymphoblastic leukemia (ALL) is the most common cancer in children, comprising approximately |$25\%$| of all childhood malignancies. ALL is characterized by the overproduction and accumulation of immature, abnormal white blood cells (WBCs; lymphoblasts) and consecutive displacement of normal haematopoiesis. Current treatment schedules for childhood ALL are based on combination chemotherapy and achieve long-term survival in >90% of children (Hunger et al., 2015). With some international variation, all major treatment protocols start with intensive, high-dose treatment for approximately 6 months (so-called induction and consolidation therapy) followed by less-intensive, low-dose treatment (so-called maintenance therapy) until 2–3 years after disease onset. While severe therapy-induced myelosuppression and frequent associated hospitalizations are acceptable up to a certain level during lymphoblast elimination in intensive treatment periods, maintenance therapy aims to achieve sustained antileukemic activity against lymphoblasts below the limit of detection, with minimal impact on quality of life due to adverse effects. The importance of maintenance therapy for patient survival has been proven and non-adherence and reduced cumulative drug exposure during maintenance therapy are associated with decreased survival (Schmiegelow et al., 2009).

The backbone of maintenance therapy includes daily oral |$ 6 $|-mercaptopurine (⁠|$ 6 $|-MP) and weekly oral methotrexate (MTX) administration. Drug doses are adjusted at the discretion of the treating physicians to achieve WBC suppression without unintended myelotoxicity according to treatment protocol-specific target ranges (i.e. the treatment protocol AIEOP-BFM 2009 specifies a WBC target range of 1.5–3.0 |$\times \ 10^3/\mu l $| and recommends dose reduction for WBC counts <|$1.5 \times 10^3/\mu l $|⁠, neutrophils <|$500/\mu l $|⁠, lymphocytes <|$300/\mu l $| and platelets <|$50\times 10^3/\mu l $|⁠) (Balis et al., 1998; Schmiegelow et al., 2014). To this end, peripheral blood sampling is performed regularly (most often weekly in the first months of maintenance therapy) and drug doses are changed accordingly. However, considerable uncertainty exists with respect to the exact drug dose response on myelosuppression and the timing of dose adjustments and resulting changes in peripheral blood counts. Both |$ 6 $|-MPand MTX are regarded as prodrugs which undergo intracellular metabolism after oral application (Schmiegelow et al., 1994, 2014). MTX is catalysed by the enzyme folylpolyglytamyl transferase to cytotoxic MTX polyglutamates (MTXPGs). Likewise, |$6$|-MP is converted to several metabolites before interfering with DNA synthesis (Lennard et al., 1992). The metabolic pathway leading to active cytotoxic |$6$|-thioguanine nucleotides (⁠|$6$|-TGNs) by the enzyme hypoxanthine-guanine phosphoribosyltransferase (HGPRT) results in the antileukaemic effect desired in maintenance therapy (Schmiegelow et al., 1995). However, large interpatient variability in bioavailability of |$6$|-MP and MTX has been reported (Kearney et al., 1979; Zimm et al., 1983; Arndt et al., 1988; Balis et al., 1998) and their major cytotoxic metabolite concentrations are associated with treatment outcomes (Lilleyman et al., 1994; Schmiegelow et al., 1995; Ebbesen et al., 2017).

In clinical practice, adjustment of maintenance therapy drug doses to desired WBC suppression levels (i.e. 1.5–3.0 |$\times 10^3 /\mu l $| WBCs) is challenging. The complex pharmacokinetics (PKs) and pharmacodynamics (PDs) of 6-MP and MTX and their interplay with a child’s bone marrow recovering from ALL induction and consolidation therapy are insufficiently represented with current treatment recommendations (i.e. reduce dose if WBCs <|$1.5 \times 10^3 /\mu l $| and increase dose if WBCs >|$3 \times 10^3 /\mu l $|⁠). We believe that personalized mathematical modelling has huge potential benefits in this context. This may lead to a better understanding of drug response (i.e. to the nonintuitive dynamics of children’s individual haematopoiesis and immune systems). Eventually, this may lead to improved therapy schedules and clinical decision support systems. However, cross-validated mathematical models with a high predictive accuracy have not yet been implemented.

Various aspects of ALL treatment during maintenance therapy have been studied in the literature. PKs and PDs of |$6$|-MP, MTX and their interaction were investigated in Balis et al. (1987, 1998); Innocenti et al. (1996); Giverhaug et al. (1999); Dervieux et al. (2002). Physiologically-based pharmacokinetic modelling of |$6$|-MP and MTX was introduced in detail in Ogungbenro et al. (2014a,b). However, this approach leads to high-dimensional and over-parameterized models. Multi-compartment models describing the PK of |$6$|-MP, MTX and their metabolism were presented in Jayachandran et al. (2014) and Panetta et al. (2002, 2010). The gold standard compartment model of Friberg et al. (2002) has been used widely to describe chemotherapy-induced leukopenia for different types of diseases (Sandstrom et al., 2005; Jayachandran et al., 2014; Rinke et al., 2016). To our knowledge, the only attempt to model the dynamics of WBCs during ALL maintenance therapy was conducted in Jayachandran et al. (2014). They investigated chemotherapy toxicity due to |$ 6 $|-MP in a virtual dataset.

In this paper, we introduce a PK/PD model of chemotherapy-induced leukopenia during maintenance therapy of childhood ALL. In particular, we consider the PK effect of the standard treatment with combined |$6$|-MP and MTX. Clinical data to personalize and cross-validate mathematical models were provided by the Department of Pediatrics and Adolescent Medicine, University Hospital Erlangen, Germany and the Department of Pediatric Hematology and Oncology, University Hospital Dresden, Germany. We show that the proposed model is able to predict the time course of WBC counts following chemotherapy using a moving horizon/iterative parameter estimation (PE) algorithm.

2. Methods

2.1 Data description

We used retrospective patient data from children diagnosed with de novo ALL at the University Hospitals in Erlangen and Dresden and treated according to the AIEOP-BFM |$ 2000 $| and |$ 2009 $| protocols. 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 did start maintenance therapy (i.e. did not relapse before the end of consolidation therapy and did not undergo stem cell transplantation). During AIEOP-BFM |$ 2000 $| and |$ 2009 $| maintenance therapy, patients received oral chemotherapy with daily |$ 6 $|-MP and once-weekly MTX until |$ 2 $| years after ALL diagnosis. Timing and dosage of chemotherapy were adjusted by the treating physicians to keep WBC counts within a protocol-specified target range of 1.5–3.0 |$\times \ 10^3/\mu l $|⁠. Dosage was reduced when cell counts fell below lower limits (WBC count <|$1.5 \times 10^3/\mu l $|⁠, neutrophils <|$500/\mu l $|⁠, lymphocytes <|$300/\mu l $| and platelets <|$50\times 10^3/\mu l $|⁠) or liver toxicity was suspected. For each patient included in the analysis, the following variables were recorded: prescribed |$6$|-MP and MTX dosage (absolute and per body surface area), WBC count, platelet count, lymphocyte and neutrophil counts and therapy interruptions.

After exclusion of patients with complications due to infections and irregular and/or missing measurements, a subset of nine patients (five boys, age range 2–9) was selected using visual inspection (Table 1). Figure 1 shows the corresponding WBC counts, the target range and chemotherapy schedules.

Fig. 1.

Underlying clinical data with treatment days on the x-axis [days]. Rows 1, 3 and 5 show measured WBC counts in [|$ \times 10^3/\mu l $|] with the therapeutic target range of 1.5–3.0 |$\times 10^3/\mu l $| in grey. Rows 2, 4 and 6 show the corresponding absolute drug doses of 6MP (⁠|$\cdot $|⁠) and MTX (⁠|$+$|⁠) in [mg]. A high variability of WBC counts and a considerable proportion of WBC counts outside the target range are observed.

2.2 Development of mathematical models

In this section, the PK models of |$ 6 $|-MP, MTX and their corresponding cytotoxic metabolites |$6$|-TGNs and MTXPGs and the resulting PK/PD model of chemotherapy-induced leukopenia are introduced. The models comprise systems of ordinary differential equations (ODEs) together with initial conditions. The concentrations of |$6$|-TGNs and MTXPGs used to represent their cytotoxic effect are fully incorporated into the PK/PD model. We tested several different mathematical models, varying initial values and penalizing deviations of parameter values from reference values. Numerical results were evaluated on small-time horizons with respect to visual assessment of solutions. The chosen model, which we describe in the following, performed best and was eventually cross-validated with a moving horizon approach over the full-time horizon.

Fig. 2.

Schematic |$ 6 $|-MP and |$6$|-TGN model.

2.2.1 |$ 6 $|-MP and |$6$|-TGN model

A simplified version of the compartment model in Jayachandran et al. (2014) is used to describe PK of |$6$|-MP and |$6$|-TGNs. A diagram of the model is depicted in Fig. 2. In particular, |$6$|-MP is absorbed at the rate |$k_a$| into plasma after oral intake to the gastrointestinal (GI) tract. From plasma, it is partly excreted at the rate |$k_e$| and gets into red blood cells (RBCs). Here, it undergoes intracellular metabolism. |$6$|-MP is metabolized by HGPRT and other enzymes at the rates |$k_{pt}$| and |$k_{pm}$|⁠, respectively. Finally, |$6$|-TGNs is assumed to be eliminated from RBC at the rate |$k_{te}$|⁠. Since the |$6$|-MP concentration in RBCs is negligible (Hawwa et al., 2008), we assumed that it is metabolized as soon as it enters RBC. Moreover, only the metabolic pathway leading to |$6$|-TGNs by HGPRT is modelled in detail since |$6$|-TGNs are the primary mediators of the cytotoxic effect of |$6$|-MP through their incorporation as false nucleotides into DNA (Schmiegelow et al., 1995). The mathematical model is

$$\begin{equation} \begin{aligned} \dot x_1 &= -k_a x_1 + \frac{\alpha F D_{6MP}(t)}{T_{dur}} ,\\ \dot x_2 &= k_ax_1-k_ex_{2}-\frac{k_{pt}(1-e_{rel})x_2}{K_t+x_2}-\frac{k_{pm}e_{rel}x_2}{K_m+x_2},\\ \dot x_3 &=\frac{ v_{pt}k_{pt}(1-e_{rel})x_2}{K_t+x_2}-k_{te}x_3 \end{aligned} \end{equation}$$

(2.1)

with initial values

$$\begin{equation} x_1(0) = x_2(0)= x_3(0)=0. \end{equation}$$

(2.2)

The |$6$|-MP control function |$D_{6MP}(t)$| is defined as

$$\begin{equation*} D_{6MP}(t) = \begin{cases} D_i, t\in [t_i,t_i+ T_{dur}] \text{ if an amount of {$6$}-MP dose {$D_i$} was taken at time } t_i,\\ \textrm{otherwise } 0. \end{cases} \end{equation*}$$

Notice that |$v_{pt}$| is used only for unit consistency. Since our clinical data does not contain concentration measurements of |$6$|-MP and |$6$|-TGNs, most values of parameters appearing in the |$ 6 $|-MP and |$6$|-TGNs model (2.1) are taken from Jayachandran et al. (2014) and are used for all patients. Simulations show that |$T_{dur}$| does not have a strong effect on the concentration of |$6$|-TGNs in RBCs. Moreover, due to large interpatient variability in bioavailability of |$6$|-MP (see Zimm et al., 1983; Arndt et al., 1988; Balis et al., 1998), |$F$| and |$T_{dur}$| are assigned to values of |$0.45$| and |$1/24$|⁠, respectively. These values will be used for MTX as well for the same reason. All state variables and parameters of (2.1) are described and summarized in Tables 2 and 3.

2.2.2 MTX and MTXPG model

The mathematical model describing PK of MTX and its metabolites |$ \mathrm{MTXPG}_{i} $| (⁠|$i=1,\ldots ,7$| is the number of glutamates attached to each MTX molecule) is based on the previous work of Panetta et al. (2002, 2010), where its detailed description and assumption can be found. A schematic illustration of the model is displayed in Fig. 3. The extracellular PK of MTX after oral intake to the GI tract were presented by a two-compartment model in Panetta et al. (2002, 2010), but were not fully described. In this work, it is simplified and encompassed by a one-compartment model with first-order absorption, see the upper part of Fig. 3 or the first two equations in (2.3). The simplification is explained in Remark 2.1. The mathematical model is the system of ODEs

$$\begin{equation} \begin{aligned} \dot x_4 &= -k_ax_4 + \frac{\beta FD_{MTX}(t)}{T_{dur}\,\mathrm{BSA}},\\ \dot x_5 &=k_ax_4-k_e x_5,\\ \dot x_6 &=\frac{V_{mI}x_5/V}{K_{mI} + x_5/V}+\frac{k_px_5}{V}-K_{eff}x_6-\frac{V_{m-fpgs}x_6}{K_{m-fpgs} + x_6}+K_{ggh}x_7,\\ \dot x_7 &=\frac{V_{m-fpgs}x_6}{K_{m-fpgs} + x_6}-K_{ggh}x_7 \end{aligned} \end{equation}$$

(2.3)

Fig. 3.

Schematic MTX and MTXPG model.

with initial values

$$\begin{equation} x_4(0) = x_5(0)= x_6(0)= x_7(0)=0. \end{equation}$$

(2.4)

The definition of the drug control function |$D_{MTX}(t)$| is similar to that of |$ D_{6MP}(t)$| and following (Du Bois et al., 1989), we calculate the body surface area as |$\mathrm{BSA} = W^{0.425} H^{0.725} 71.84$| with a patient’s weight |$W$| and height |$H$|⁠. The last term in the first equation of (2.3) is divided by |$\mathrm{BSA}$| for unit consistency.

As in the |$6$|-MP case, measurements of MTX and |$\mathrm{MTXPG}_{1-7}$| are not available. We set most model parameters in (2.3) to values from the literature (Panetta et al., 2010, Table |$2$|⁠). The values of |$k_a,k_e$| are obtained via PE based on measurements of MTX concentration reported in Pinkerton et al. (1980). All state variables and parameters of (2.3) are summarized in Tables 4 and 5.

 

Remark 2.1

In this work, we investigated the two-comparment model for MTX from Panetta et al. (2010) and its simplified (one-compartment) version. We concluded that the obtained results of the PK/PD model did not differ significantly between both PK models since MTX in plasma (only this compartment has a direct effect on the intracellular compartments) was manipulated by the changes of |$k_a$| and |$k_e$|⁠, accordingly. The drug concentrations for the considered patients are not available. Moreover, the focus of this paper is not on estimating MTX parameters, but on encompassing several models to describe the treatment outcomes during ALL maintenance therapy. Therefore, we decided to use a one-compartment model for simplicity. Regarding the estimation of |$k_e$| and |$k_a$| (based on the measurements in Pinkerton et al., 1980), first the elimination half-life was fixed to a value from a range provided in the literature, e.g. Bleyer (1978). Then, the area under curve (AUC) was computed by the trapezoid rule (based on the provided measurements) and extrapolation to obtain AUC from the last known concentration to infinity.

2.2.3 Model of chemotherapy-induced leukopenia

We present a mathematical model that is based on the gold-standard myelosuppression model (Friberg et al., 2002) that has previously been successfully applied in the context of chemotherapy treatment (Kloft et al., 2006; Netterberg et al., 2012; Quartino et al., 2012; Pefani et al., 2013, 2014; Rinke et al., 2016). We use five compartments, see Fig. 4. The proliferating compartment describes stem cells and progenitor cells which are sensitive to chemotherapeutic agents and are able to self-renew at the rate |$k_{prol}$|⁠. The next three transit compartments describe the maturation process of cells and are used to account for the delay time between observed effect and drug administration. The transition rate between compartments is set to |$k_{tr}$|⁠. The last compartment illustrates circulating or mature WBC, which die at the rate |$k_{circ}$| due to apoptosis. Cells in the last four compartments are assumed to be drug insensitive. The feedback mechanism |$\displaystyle \left(\frac{\textrm{Base}}{x_{12}}\right)^{\gamma }$| manipulates the mitosis rate by taking the information about the number of circulating cells into account. During ALL maintenance therapy, blood samples are commonly taken once per week or less frequently, resulting in a sparse dataset. To improve the identifiability of model parameters, |$k_{circ}$| is assigned to |$k_{circ}=0.5346/\mathrm{day}$| as in Jayachandran et al. (2014) since it is less sensitive in comparison with the other parameters in (2.5). Furthermore, |$k_{prol} = k_{tr}$| at steady state, therefore |$k_{prol}$| is set to |$k_{tr}$|⁠. The loss of proliferating cells due to chemotherapy is included by a log-linear function |$E_{drug}$|⁠. Since the active intracellular metabolite concentrations of MTX (MTXPGs) and 6MP (6-TGNs) are associated with the treatment outcome of ALL (Lilleyman et al., 1994; Schmiegelow et al., 1995; Ebbesen et al., 2017), they are incorporated into the WBCs’ dynamics via |$E_{drug}$| instead of their direct drug concentrations in plasma. Additionally, long-chain |$ \mathrm{MTXPGs} $| are known to be more likely to exert antileukemic effects than MTX and short-chain |$ \mathrm{MTXPGs} $| (Masson et al., 1996). For a two-drug schedule, the contribution of the component drugs to the total effects remains unclear. Thus, it is assumed to contribute equally in this work. Therefore, |$ E_{drug} = \ln \ (\textrm{slope} \ (\,p_{6MP}x_3 + p_{MTX}x_7 )/1000 + 1) $| is chosen with the constants |$p_{6MP} = 1$| and |$p_{MTX}=1$| for unit consistency. The model can be written as (2.5). The detailed description of the state variables and model parameters are listed in Tables 6 and 7.

$$\begin{equation} \begin{aligned} \dot{x_8}& = k_{prol}\,x_8 (1-E_{drug})\left(\frac{\textrm{Base}}{x_{12}}\right)^{\gamma}-k_{tr}\,x_8,\\ \dot{x_9}& = k_{tr} (x_8- x_9),\\ \dot{x_{10}} &= k_{tr} (x_9- x_{10}),\\ \dot{x_{11}}& = k_{tr} (x_{10}- x_{11}),\\ \dot{x_{12}}& = k_{tr} x_{11}- k_{circ} x_{12}.\\ \end{aligned} \end{equation}$$

(2.5)

Fig. 4.

Schematic model of chemotherapy-induced leukopenia.

We assume that the WBCs oscillate around the steady-state value. Therefore, the initial values of (2.5) are commonly chosen as follows (Quartino et al., 2012):

$$\begin{equation} \begin{aligned} x_8(0)&=x_9(0)=x_{10}(0)=x_{11}(0)=\frac{\textrm{Base}*k_{circ}}{k_{tr}}, \quad x_{12}(0)& = \textrm{Base}. \end{aligned} \end{equation}$$

(2.6)

2.3 Parameter estimation

The dataset described in Section 2.1 contains patients’ physiological data and their WBC counts over time. This allows us to estimate patient-specific parameters of (2.5) by fitting predicted WBC counts to real-world measurements. A general mathematical formulation of a PE problem is typically stated as

$$\begin{equation} \begin{aligned} \min_{\theta } F(x,\theta) \quad \quad \mathrm{s.t.}\quad \dot{x}(t) = f(x(t),\theta), \quad x(0) = x_0, \end{aligned} \end{equation}$$

(2.7)

where |$\theta $| is a vector of unknown model parameters which need to be estimated. The objective functional |$F(x,\theta )$| models deviations between model response and real-world measurements. Optimal model parameters |$\theta ^{\ast }$| are obtained as a solution of (2.7).

2.4 Adaptive PE approach

The clinical dataset contains weekly or bi-weekly WBC measurements. As soon as new measurements arrive, they can be taken into consideration to update model parameters |$\theta $|⁠. We propose to use an adaptive PE procedure to allow model parameters |$\theta $| to slightly vary over time. This is meant to account for systematic changes over time due to children’s physiological development with age and dynamics after long-term chemotherapy. Practically, we only consider the most recent |$N_m$| measurements. Alternatively, monotonically increasing scaling factors (in time) could be used. In every iteration |$k = 1,2,\ldots $|⁠, we solve a PE problem over the time horizon |$[t^{k}_0,t^{k}_f]$|

$$\begin{equation} \min_{\theta^{k}} \,\,\sum_{i=1}^{N_m} \, \left(\frac{x_{12}(t^{k}_i; \theta^{k})-m^{k}_i}{\sigma_i}\right)^2 + \sum_{i = 1}^{4} \varepsilon_i^k \left(\theta^{k-1}_i - \theta^{k}_i\right)^2 \quad\mathrm{s.t.} \quad (2.1, 2.3, 2.5), \quad x\big(t^{k}_0\big) = x^{k}_0. \end{equation}$$

(2.8)

The first term in the objective models is the least squares error between |$N_m$| measured WBC counts |$m_i^{k}$| at times |$t^{k}_i$| and the corresponding calculated model response |$x_{12}(t^{k}_i; \theta ^{k})$|⁠. This model response depends on the unknown personal model parameters |$\theta ^k = (\textrm{slope}^k,k_{tr}^{k},\textrm{Base}^k,\gamma ^k)$| in iteration |$k$|⁠. The variances are chosen as |$\sigma _i=\sigma = 0.04\, M$| with |$i=1,\ldots ,N_m$| and mean value |$M$|⁠. To avoid unrealistically large changes in the model parameters from one iteration to the next, we penalize deviations between |$\theta ^{k-1}$| and |$\theta ^k$| with a second term in the objective. We use the zero vector as |$\theta ^0$| to favor parameters with small values. Once the methodology has been established and applied to several patients, population-based parameter values can be used to initialize |$\theta ^0$|⁠. The penalized vector is chosen as |$\varepsilon ^1 = (0.1,0.1,0.1,0.1)$| and |$\varepsilon ^k = (20,0.5,0.5,0.5)$| with |$k>1$| to produce a stable evolution of the estimated parameters. The initial values of the state variables |$x$| in (2.1, 2.3, 2.5) are set to fixed values |$x^k_0$|⁠. For |$k=1$|⁠, the initial values specified in (2.2, 2.4, 2.6) are used. For |$k>1$|⁠, they are set to the corresponding state values |$x_i(t^{k-1}_1;\theta ^{k-1})$| of the previous iteration |$k-1$|⁠. Note that |$t^{k}_0=t^{k-1}_1$|⁠.

The ODE system (2.1), (2.3) and (2.5) is discretized with a fourth-order Runge–Kutta scheme. The resulting PE problem (2.8) should be solved with a Generalized Gauß–Newton method (Bock, 1987). In our prototype implementation, we used the sequential quadratic programming method in Matlab’s fmincon. The main steps of the adaptive PE algorithm are summarized in Algorithm 2.1.

 

Algorithm 2.1

FOR |$ k=1,2,\ldots $|

• (a) Set up (2.8) with the measurement set |$M_{N_m} = \{m^k_1,\ldots ,m^k_{N_m}\}$| and calculate initial values |$x(t^k_0)$|⁠.

• (b) Solve (2.8) and obtain updated parameters |$\theta ^k$|⁠.

• (c) Predict future WBC dynamics.

There is a vast existing literature on the adaptive PE approach or in other words the moving horizon approach. The readers can refer to Kühl et al. (2011) and references therein for more details.

2.5 Statistical methods

To assess how well the model (2.5) is able to capture the time course of WBC counts with the estimated values of the patient-specific parameters obtained as a solution of (2.8), simulations over a prediction horizon are conducted to produce the predicted WBC counts over time. The predictability of (2.5) is evaluated by comparing model response with real-world measurements at the corresponding time points. In particular, the standard deviation (Sd) of the differences between the real-world measurements and the predicted WBC counts is computed for each patient. Moreover, the paired two-sample |$t$|-test is carried out to calculate |$\mathrm{p}$|-values (see e.g. Jeremy et al., 2014). Model response and real-world measurements are two paired samples. In this scenario, we want to test the assumption that the chemotherapy-induced leukopenia model (2.5) is able to predict WBC count over a certain time horizon. Therefore, the null hypothesis is that the differences between the paired samples were drawn from a normal distribution with zero mean.

Fig. 5.

Algorithm 2.1 applied to virtual data. Virtual WBC counts [|$ \times 10^3/\mu l $|] and estimated state trajectories |$x_{12}(t)$| for measurement days on the x-axis (daily monitoring) are shown. Estimated (red line) and predicted (green line) trajectories almost coincide with the true profile (black line). The left and right (zoomed-in) figures display the results after one and five iterations, respectively. The adaptive PE approach remains stable over iterations and the uncertainty (shaded) region covers most measurements. The x- and y-axes represent measurement schedule [days] and WBC count [|$ \times 10^3/\mu l\, $|], respectively.

Fig. 6.

Residual errors (RSME) of PE (left) and prediction (right) for different measurement densities and errors that were added to the simulated data. Model–data mismatch increases from retrospective estimation (left) toward predictive simulation (right) and from daily monitoring toward weekly monitoring of WBC counts.

3. Numerical results

There are several types of uncertainty involved, when Algorithm 2.1 is applied to clinical data including measurement errors, measurement frequency and modelling and algorithmic (numerical) errors. We consider both virtual and clinical data to get a better understanding of their influence.

3.1 Virtual data

We generate a virtual dataset with the following procedure:

• (i) Take a random patient from our real-world dataset.

• (ii) Estimate parameters using |$N_m$| measurements within the first |$ 100 $| days to obtain |$\theta _{\textrm{true}}$|⁠.

• (iii) Simulate (2.1, 2.3 and 2.5) over |$160$| days with |$\theta _{\textrm{true}}$| as the true parameters to generate the true data.

• (iv) Add normally distributed noise to this true data to obtain the virtual data.

With this virtual patient, it is expected that the estimated parameters remain almost unchanged against time. Moreover, the more often measurements are drawn, the more robust the model parameters are and the better the model (2.5) can predict the future. Figure 5 displays the performance of the adaptive PE approach applied to the virtual data with daily monitoring frequency and the measurement error |$\varepsilon \in N(0,0.25)$|⁠. As illustrated, the estimated and predicted trajectories almost coincide with the true profile (the solid black line). The approach remains stable over iterations and the uncertainty region covers most of the measurements. The root mean square error (RMSE) behaviour of the adaptive PE approach under different uncertainty levels and monitoring frequencies is depicted in Fig. 6. Observe that the errors of PE and prediction increase with data uncertainty and less frequent monitoring.

Fig. 7.

Values of model parameters |$\theta ^k$| during iterations |$k$| of Algorithm 2.1. Most of the parameters are either constant or show a linear change over time. The steady-state baseline WBC value |$Base$| varies stronger than the other parameters. This may be related to recovery of the bone marrow with time after the completion of induction therapy and/or age-dependent changes of normal WBC ranges. The x- and y-axes represent time [days] and values of the estimated parameters |$k_{tr}\,\,[1/\mathrm{day}]$|⁠, |$\gamma $|⁠, |$\mathrm{Base}\,\,[\times 10^9/l]$|⁠, |$\mathrm{slope}\,\,[1/\mathrm{pmol}],$| respectively.

Fig. 8.

Model predictions (red) and real-world measurements of WBC counts (black). The uncertainty tubes (the shaded regions) are obtained from |$ 100 $| simulations using parameters sampled from |$95\%$| of the normal distribution with the estimated parameters and uncertainties as its means and Sds. The predictions were carried out over a prediction horizon of |$10$| days. The x- and y-axes represent measurement schedule [days] and WBC count [|$ \times 10^3/\mu l $|] respectively.

3.2 Clinical data

It is known that intraindividual variation in bioavailability or stress is significant, which results in measurement data with uncontrollable noise. Therefore, the application of the adaptive PE technique to our real-world data becomes more demanding than the ideal case considered in Section 3.1.

The time profiles of the model parameters obtained via the approach in Section 2.4 for all patients in our real-world dataset are displayed in Fig. 7. Most of the parameters are either constant or show a linear change over time. The steady-state baseline WBC value, |$Base$|⁠, varies stronger than the other parameters. This may be related to not yet fully understood processes including long-term recovery of the bone marrow after induction therapy on the one hand and a normal evolution of |$Base$| due to aging on the other. With the parameter values as a solution of (2.8), we carry out simulations over a prediction horizon (i.e. |$10$| days) to produce the predicted WBC counts. The results are shown in Fig. 8. The uncertainty regions cover most of the measurements excluding outliers. The Sd of the differences between the real-world measurements and the predicted WBC counts and |$\mathrm{p}$|-value derived from the paired two sample t-test (as described in Section 2.5) for each patient are depicted in Table 8. The |$p$|-values obtained from the considered data did not support rejection of the null hypothesis, the predictability of the model (2.5), at the significant level of |$5$|%. More evidence and a more detailed study with additional clinical data is needed to clarify the issue (see e.g. Jeremy et al., 2014). Since blood samples were drawn at most weekly, the number of data points is not large enough to produce more intuitive scatter plots and rigorous justification, see e.g. Patients 7 and 8. Figure 9, which was obtained by putting together the results of the nine patients, shows that the WBC counts in the therapeutic target (shaded region) concentrate nicely around the diagonal in most cases.

Fig. 9.

Real-world measurements versus predicted WBC counts (⁠|$ \times 10^3/\mu l $|⁠) of all patients. The shaded area shows the WBC target range.

4. Discussion and summary

In this paper, we present a PK/PD model for chemotherapy-induced leukopenia of pediatric ALL patients during maintenance therapy. The model comprises the PK of |$6$|-MP, MTX and of their cytotoxic metabolites |$6$|-TGNs and MTXPGs, respectively. The drug effects of |$6$|-MP and MTX are described as a log-linear PD function which links the PK model with the myelosuppression model. To our knowledge, this is the first time that both drugs are taken into consideration in a data-driven modelling approach. Due to the lack of drug-related measurements, most parameters of (2.1) and (2.3) were fixed to literature values. This may contribute to the variations of the estimated parameters in the model (2.5) since inter-individual variation in drug dose-response is well established. Therefore, patient-specific measurements are indispensable to develop and validate mathematical models for further purposes, although the PD effects of orally taken drugs are more difficult to determine compared with those of intravenous administrations.

The parameters of (2.5) were obtained by fitting the model response to real-world measurements of each patient available in our dataset, which are therefore patient-specific. By means of the adaptive PE approach, each iteration produced a set of model parameters, which remained stable over time (equivalently updating of the measurement subset used for estimating the parameters) as expected. Furthermore, the predictability of our model over a desired time duration was verified using a real dataset by the suitable statistical methods in Section 2.5. This is especially important if the proposed model is to be refined and extended to support clinical decisions regarding dose adjustments in childhood ALL maintenance therapy. It is also important to point out that some parameters of (2.5) are kept fixed to improve their identifiability due to the sparsity of our data. Therefore, optimal experiment design (OED) should be considered to overcome this issue (Jost et al., 2017). OED is able to suggest the most informative measurement time points. Performing OED leads to reduction of parameter uncertainties and optimization of sampling times, which results in better monitoring of patients.

In summary, we provide a comprehensive model-based procedure to describe chemotherapy-induced leukopenia during childhood ALL maintenance therapy, which includes a combination of |$6$|-MP and MTX. The model predicts WBC counts surprisingly well given the large variation of individual response patterns in the clinical data. We see the proposed mathematical model and algorithmic procedure as important first steps on the way toward personalized clinical decision support in childhood ALL maintenance therapy.

Acknowledgements

Financial support from the European Research Council via the Consolidator Grant MODEST-647573 is gratefully acknowledged. Thanks to Rebecca Terry for proofreading the article.

Appendix

Table 1

Demographic data of nine patients

Patient	Gender	Age	Risk group	Outcome	Treatment protocol
|$ 1 $|	M	|$ 9 $|	MR	Remission	AIEOP BFM |$ 2000 $|
|$ 2 $|	F	|$ 2 $|	SR	Remission	AIEOP BFM |$ 2000 $|
|$ 3 $|	F	|$ 3 $|	MR	Remission	AIEOP BFM |$ 2000 $|
|$ 4 $|	F	|$ 5 $|	HR	Remission	AIEOP BFM |$2009$|
|$ 5 $|	M	|$ 2 $|	MR	Remission	AIEOP BFM |$2000$|
|$ 6 $|	F	|$ 2 $|	MR	Remission	AIEOP BFM |$2000$|
|$ 7 $|	M	|$ 6 $|	MR	Remission	AIEOP BFM |$ 2009 $|
|$ 8 $|	M	|$ 3 $|	SR	Remission	AIEOP BFM |$ 2009 $|
|$ 9 $|	M	|$ 2 $|	SR	Remission	AIEOP BFM |$2000$|

SR: standard risk, MR: medium risk, HR: high risk

Table 2

State variables of the |$ 6 $|-MP and |$6$|-TGN model

Variables	Units	Description
|$ x_1 $|	|$ \mathrm{pmol} $|	Amount of |$6 $|-MP in GI tract
|$ x_2 $|	|$ \mathrm{pmol} $|	Amount of |$ 6 $|-MP in plasma
|$ x_3 $|	|$\mathrm{pmol}/8\times 10^8\, \textrm{RBCs}$|	Concentration of |$ 6 $|-TGN in RBCs

Table 3

Parameters of the |$ 6 $|-MP and |$6$|-TGN model

Parameters	Values	Units	Description
|$ k_a $|	|$ 4.8 $|	|$ 1/\mathrm{day} $|	|$ 6 $|-MP absorption rate from GI tract
|$ k_e $|	|$ 5.0 $|	|$ 1/\mathrm{day} $|	|$ 6 $|-MP elimination rate from plasma
|$ k_{pt} $|	|$ 29.8 $|	|$ \mathrm{pmol} $||$ 6 $|-MP/day	|$ 6 $|-MP to |$ 6 $|-TGN conversion rate
|$ k_{pm} $|	|$ 655.8 $|	|$ \mathrm{pmol} $||$ 6 $|-MP/day	|$6$|-MP to MeMP Conversion Rate
|$ K_{t} $|	|$ 4.04 \times 10^5 $|	|$ \mathrm{pmol} $|	M-M constant for |$ 6 $|-TGN
|$ k_{te} $|	|$ 0.0714 $|	|$ 1/\mathrm{day} $|	|$ 6 $|-TGN elimination rate from RBCs
|$ e_{rel} $|	|$ 0.5 $|		TPMT enzyme activity constant
|$ v_{pt} $|	|$ 1 $|	|$\displaystyle \frac{\mathrm{pmol}\,\,6 \mathrm{-TGN}}{\mathrm{pmol}\,\,6\mathrm{-MP}/ 8 \times 10^8\, \mathrm{RBCs} } $|	|$ 6 $|-TGN elimination rate from RBCs
|$F$|	|$0.45$|		Bioavailability factor
|$T_{dur}$|	|$1/24$|	|$\mathrm{day}$|	Time duration for drug absorption
|$D_{6MP}(t)$|		|$\mathrm{mg}$|	|$6$|-MP control function
|$\alpha $|	|$10^{12}/152177 $|	|$\mathrm{pmol}/\mathrm{mg}$|	Unit consistency constant

MeMP: methyl-mercaptopurines, TPMT: thiopurine methyl-transferase

Table 4

State variables of the MTX and MTXPG model

Variables	Units	Description
|$ x_4 $|	|$ \mu \mathrm{mol}/\mathrm{m}^2 $|	Amount of MTX in GI tract
|$ x_5 $|	|$ \mu \mathrm{mol}/\mathrm{m}^2 $|	Amount of MTX in plasma
|$ x_6 $|	|$\mathrm{pmol}/10^9$| cells	Intracellular concentration of |$ \mathrm{MTXPG}_1 $|
|$ x_7 $|	|$\mathrm{pmol}/10^9$| cells	Intracellular concentration of |$ \mathrm{MTXPG}_{2-7} $|

Table 5

Parameters of the MTX and MTXPG model

Parameters	Values	Units	Description
|$ k_a $|	|$ 26.64 $|	|$ 1/\mathrm{day} $|	MTX absorption rate from GI tract
|$ k_e $|	|$ 5.76 $|	|$ 1/\mathrm{day} $|	MTX elimination rate from plasma
|$ k_p $|	|$ 9.6 $|	|$ 1/\mathrm{day} $|	Passive influx rate
|$ V $|	|$ 11.606 $|	|$ \mathrm{L}/\mathrm{m}^2$|	Systematic volume
|$V_{mI} $|	|$ 2.3895\times 10^4 $|	|$ \mathrm{pmol}/10^9 \mathrm{cells}/\mathrm{day}$|	M-M parameter for active (RFC) influx
|$ K_{mI} $|	|$2.898 $|	|$ \mu \mathrm{M} $|	M-M parameter for active (RFC) influx
|$ K_{eff} $|	|$ 179.76 $|	|$ 1/\mathrm{day} $|	first-order efflux parameter
|$ V_{m-fpgs} $|	|$ 7.0119\times 10^3 $|	|$ \mathrm{pmol}/10^9 \mathrm{cells}/\mathrm{day}$|	M-M parameter for FPGS activity
|$ K_{m-fpgs}$|	|$35.262$|	|$ \mathrm{pmol}/10^9 \mathrm{cells}$|	M-M parameter for FPGS activity
|$ K_{ggh}$|	|$4.992$|	|$ 1/\mathrm{day} $|	First-order GGH activity
|$F$|	|$0.45$|		Bioavailability factor
|$T_{dur}$|	|$1/24$|	|$\mathrm{day}$|	Time duration for drug absorption
|$ \mathrm{BSA} $|		|$ \mathrm{m}^2 $|	Body surface area
|$D_{MTX}(t)$|		|$\mathrm{mg}$|	MTX control function
|$\beta $|	|$10^{6}/454440 $|	|$\mu \mathrm{mol}/\mathrm{mg}$|	Unit consistency constant

RFC: reduced folate carrier, M–M: Michaelis–Menten

Table 6

State variables of chemotherapy-induced leukopenia model

Variables	Units	Description
|$ x_8 $|	|$\sharp \times 10^9/\mathrm{l}$|	Number of proliferating cells
|$ x_9,x_{10},x_{11} $|	|$\sharp \times 10^9/\mathrm{l}$|	Number of maturating cells
|$ x_{12} $|	|$\sharp \times 10^9/\mathrm{l}$|	Number of circulating cells

Table 7

Parameters of the chemotherapy-induced leukopenia model

Parameters	Values	Units	Description
|$ k_{prol} $|		|$ 1/\mathrm{day} $|	Proliferation rate
|$ k_{tr} $|		|$ 1/\mathrm{day} $|	Transition rate
|$ k_{circ} $|	|$0.5346$|	|$ 1/\mathrm{day} $|	Apoptosis rate
|$ \gamma $|			Feedback factor
|$ \textrm{Base} $|		|$\sharp \times 10^9/\mathrm{l}$|	Baseline value of WBC count
|$ \textrm{slope} $|		|$\mathrm{l}/\mathrm{pmol}$|	Drug effect coefficient
|$ p_{6MP}$|	|$1$|	|$8\times 10^8 \mathrm{RBCs}/l$|	Constant for unit consistency
|$ p_{MTX} $|	|$1$|	|$10^9 \mathrm{cells}/l$|	Constant for unit consistency

Table 8

The Sd of the differences between the real-world measurements and the predicted WBC counts and |$\mathrm{p}$|-value derived from the paired two sample t-test (as described in Section 2.5) for each patient. The p-values show that the null hypothesis, the predictability of the model (2.5), cannot be rejected at the |$ 5\% $| significance level for all considered patients

Patient	|$ 1 $|	|$ 2 $|	|$ 3 $|	|$ 4 $|	|$ 5 $|	|$ 6 $|	|$ 7 $|	|$ 8 $|	|$ 9 $|
Sd	|$0.797$|	|$0.744$|	|$0.691$|	|$1.187$|	|$0.61$|	|$0.605$|	|$1.369$|	|$0.487$|	|$0.479$|
p-value	|$0.121$|	|$0.669$|	|$0.223$|	|$0.022$|	|$0.172$|	|$0.178$|	|$0.301$|	|$0.059$|	|$0.798$|