http://orcid.org/0000-0002-0316-7055
Figures
Abstract
We investigate the personalisation and prediction accuracy of mathematical models for white blood cell (WBC) count dynamics during consolidation treatment using intermediate or high-dose cytarabine (Ara-C) in acute myeloid leukaemia (AML). Ara-C is the clinically most relevant cytotoxic agent for AML treatment. We extend a mathematical model of myelosuppression and a pharmacokinetic model of Ara-C with different hypotheses of Ara-C’s pharmacodynamic effects. We cross-validate the 12 model variations using dense WBC count measurements from 23 AML patients. Surprisingly, the prediction accuracy remains satisfactory in each of the models despite different modelling hypotheses. Therefore, we compare average clinical and calculated WBC recovery times for different Ara-C schedules as a successful methodology for model discrimination. As a result, a new hypothesis of a secondary pharmacodynamic effect on the proliferation rate seems plausible. Furthermore, we demonstrate the impact of treatment timing on subsequent nadir values based on personalised predictions as a possibility for influencing/controlling myelosuppression.
Citation: Jost F, Schalk E, Rinke K, Fischer T, Sager S (2019) Mathematical models for cytarabine-derived myelosuppression in acute myeloid leukaemia. PLoS ONE 14(7): e0204540. https://doi.org/10.1371/journal.pone.0204540
Editor: Obul Reddy Bandapalli, German Cancer Research Center (DKFZ), GERMANY
Received: September 5, 2018; Accepted: May 30, 2019; Published: July 1, 2019
Copyright: © 2019 Jost et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: The data which were used to perform the parameter estimations are deposited publicly in form of NONMEM datasets and a Python file. The data sets are available at https://bitbucket.org/fjost21/jost2017data/src/default/. The data sets are also available in the Supporting Information (S1, S2 and S3 Files).
Funding: This project has received funding from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation programme (grant agreement No 647573) and from the “International Max Planck Research School (IMPRS) for Advanced Methods in Process and System Engineering” in Magdeburg.
Competing interests: The authors have declared that no competing interests exist.
Introduction
Acute myeloid leukaemia (AML) is a malignant clonal disorder of myeloid stem and progenitor cells. In untreated AML, immature neoplastic myeloid blasts rapidly proliferate and suppress the generation and maturation of blood cells in the bone marrow. While being a curable disease using chemotherapy including anthracyclines and/or cytarabine (Ara-C), this approach leads to prolonged myelosuppression with extremely low white blood cell (WBC) counts (leukopenia), i.e. values below 1 G/L, associated with a high risk of infection and treatment-related mortality [1].
Consolidation treatment, repetitive (up to 4) cycles of intermediate-/high-dose Ara-C (1 − 3 g/m2) [2], is given once patients achieve complete remission (CR) and is considered the most important part of chemotherapy in preventing relapses. The treatment of 3 g/m2 Ara-C infusion lasting 3 hours every 12 hours on days 1, 3 and 5 for patients aged 60 years and younger was established by Mayer et al. [3] and remains the current standard of treatment to this day.
If predictions from personalised mathematical models including all relevant biomarkers were reliable and accurate, they could be used for providing better care to AML patients receiving Ara-C consolidation treatment, e.g. in an automatised measurement–decision support loop [4, 5]. Precisely identifying the period of Ara-C-induced profound leukopenia and modification of treatment schedules based on such predictions might enable prevention of severe infectious complications, sepsis, and thus delay to undergo subsequent treatment cycles. Therefore, by realising timely adherence to consolidation therapy cycles and by avoiding delays in treatment schedule, the density of chemotherapy cycles may be increased and thus deeper remissions and lower relapse rates may be achieved. This may ultimately translate into improved overall survival rates.
There are many different levels on which haematopoiesis [6–8], granulopoiesis [9–11], myelosuppression [12] and dynamics of leukemic cells [8] can be modelled [13, 14]. A comprehensive overview and summary of the various models is given in the recently published reviews [15, 16]. We analysed models that capture only the most important dynamics for non-leukemic cells and “agglomerate” different physiological effects into simplified expressions. Higher levels of detail in more sophisticated models, covering many physiologial properties and thus providing a deeper understanding of biological phenomena, come at the price of needing more observed biomarkers and model parameters. Depending on the medical question, the required outcome and the available biomarker information, more complex models can be reformulated into minimalistic models that concentrate on the fundamental physiological mechanisms without a qualitative loss of the outcome [17, 18]. In our setting, the current lack of clinical measurements of endogenous granulocyte-colony stimulating factor (G-CSF) concentrations and leukemic cell counts (as no relapse events occurred) leads to identifiability issues with the related dynamics. Due to these issues and our main focus on myelosuppression and WBC recovery, we concentrated on agglomerating effects of Ara-C on proliferation and maturation rates and did not consider models including G-CSF or leukemic cell dynamics.
Mathematical models for myelosuppression due to various chemotherapy agents have been proposed [12, 19–23] and applied successfully to predict the dynamics of neutrophils [4, 24]. However, this is not the fact for high-dose Ara-C, the most important component in consolidation therapy [2, 25]. Pharmacology of Ara-C is particularly difficult, as its exact mechanisms of action both on normal and leukemic cells are not fully understood. The main effect of Ara-C on normal and leukemic proliferating cells is the inclusion of intracellular Ara-C triphosphate (Ara-CTP) into DNA and RNA, which impairs cell replication [26]. Yet, the synthesis of intracellular Ara-CTP is saturable such that the clinical success of intermediate-/high-dose Ara-C is not well explained [3, 27, 28]. Additional effects are the subject of ongoing research [28, 29].
Here, we surveyed different published and new hypotheses of the pharmacodynamic (PD) effects of Ara-C on WBC dynamics during AML consolidation therapy. We used models for myelosuppression and Ara-C pharmacokinetics (PK) from the literature to quantify prediction accuracies. All model variations are based on the myelosuppression model developed by Friberg et al. [12] and are tailored to the special case of Ara-C via a parameterised two-compartment PK model. The general modelling goals were to include possible secondary effects of Ara-C and to obtain a good balance between modelling detail, prediction accuracy, and the number of patient-specific parameters. As a successful methodology, we considered predictions of WBC recovery times trec (defined as the time when the WBC count recovers above 1 G/L) for different Ara-C schedules and compared them to published average WBC recovery times.
Our work is another contribution towards the ultimate goal of mathematically optimising and individualising consolidation therapy for AML patients. Here, we are focusing on one important aspect of AML therapy: Ara-C-derived myelosuppression.
Remark on terminology and potentially confusing synonyms
The manuscript is categorised in the intersection of mathematics, control theory, systems biology, pharmacology, and medicine. Words like “model” or “parameter” have different meanings in these scientific communities, and similar concepts have different names like “calibration”, “estimation”, or “personalisation”. For convenience, we list some synonyms that we did (not) use in S1 Table.
Methods
The section starts with a detailed description of the different model variations. Next, three initial value approaches are introduced which are used during parameter estimation. Afterwards we describe our clinical data and specify how we personalised the mathematical models. The sections Prediction & Cross-Validation and Schedule Timing contain the description of several simulation studies which were performed to discriminate between the different model variations and to analyse the treatment timing in consecutive CCs on the nadir value. The mathematical approaches to parameter estimation, uncertainty quantification, statistical analysis and the nonlinear mixed-effects modelling approach can be found in the S1 Appendix.
Mathematical Models
Fig 1 illustrates the basic assumptions from which we derived twelve model variations of the original Friberg model which we denote by M1–M12 from now on. They differ concerning the number of transition compartments (M1–M3), initial conditions for the differential equations (M3–M5), and model assumptions for the possible effects of Ara-C on proliferation and maturation rates (M5–M12). In this designation, the original Friberg model is denoted by M2. The ordering of the first three models was chosen with respect to the decending number of transition compartments. Quartino et al. [30] proposed a model with six instead of three transition compartments and is thus denoted as M1. After intermediate evaluations of accuracies we concentrated on the most promising choice of scaling, transition compartments, and initial conditions, and included different modelling assumptions in the models M6–M12 which are alternatives to M5, our reference myelosuppression model extended to Ara-C. Most models refer to previous approaches in the literature and are included for a comprehensive comparison and evaluation of our new hypotheses.
We assumed clustering of cells and cytarabine (Ara-C) concentrations in compartments with identical properties. White blood cell (WBC) differentiation is represented by a proliferating compartment xpr, a number ntr of transit compartments xtr with different levels of maturation, and a compartment xma with mature, circulating WBCs. Cells mature with a maturation rate G. Mature cells xma die by apoptosis with a death rate of kma. The pharmacodynamic effect of Ara-C is described as a log-linear function E targeting the proliferating cells in the bone marrow. It depends on the concentration x1 of Ara-C in an assumed central compartment including the circulating blood. The proliferation rate F of xpr models the replication speed of proliferating progenitor cells. Modelling assumptions were incorporated by choosing different functions F and G (compare Table 1). The estimated model parameters used for personalisation were B, slope, ktr, γ, and initial conditions.
M1–M5: The basic PK/PD model, number of compartments and initial conditions.
In 2002 Friberg et al. published a PK/PD model describing myelosuppression induced by different chemotherapeutic agents (docetaxel, paclitaxel, and etoposide) [12]. The well studied model showed a good trade-off between capturing the important aspects of the dynamics, containing a moderate number of identifiable model parameters, and being applicable for different cytostatic drugs. It has become the gold-standard model in pharmacometrics [16] with different PK and population-based modifications to topotecan [31], to daunorubicin [32], to a combination therapy of Ara-C (low-dose), etoposide and daunorubicin in the induction treatment for AML [20], to a physiologically based PK model for the induction therapy of AML patients with daunorubicin and Ara-C (low-dose) [33], to a combination therapy of carboplatin, etoposide and thiotepa [34], to paclitaxel [22], to an individual-based approach [35], and to drug specific optimisations [21].
WBCs derive from differentiated, matured haematopoietic stem cells that have passed through several intermediate stages during maturation. The chain of maturation is reflected in the mathematical model as a clustering of cells in several consecutive compartments with identical properties. Each compartment is described as a differential state. The 4 + ntr differential states of our mathematical models M1–M12 are the amounts x1 and x2 of Ara-C in two PK compartments, respectively, the amounts xpr of proliferating cells, of differentiating cells in ntr transient compartments, and xma of mature, circulating WBCs.
The differential equations that correspond to Fig 1 are
(1)
(2)
(3)
(4)
(5)
(6)
(7)
with constants k10, k12, k21, BSA, duration, cV, kma, Ara-C administration u(t) and parameters ktr, γ, B and slope specified in the following model description and listed in S3 Table. Functions F and G are chosen differently in the models M1–M12, compare Table 1.
For fixed transition rate ktr, the number of compartments can be used to model the delay between the proliferating and circulating cells (mean maturation time [36]). As there is no common consensus on the precise number of differentiation stages [7, 37] we compared ntr = 6, as proposed by Nock [34], ntr = 3, as proposed by Friberg et al. [12] and ntr = 1, proposed by us in which we comprise the whole maturation process into one transition compartment for models M5-M12. The fusion of the differentiation steps into one compartment is justified by the mean maturation time (MMT) from proliferating stem cells to circulating mature WBCs, which we compared with published values and presented in the section Results. Cells mature with a maturation rate constant G = ktr summarising the fraction of cells performing self-renewal and differentiation into one parameter. This is a simplified assumption made by Friberg et al. guaranteeing homeostasis [17] and identifiability of the estimated parameters. During the modelling process we analysed a model considering separate parameters for the fraction of self-renewal and for differentiation in each compartment. The model has a similarly high accuracy but more challenging identifiability properties. Our findings are summarised in S5 Table.
Mature cells xma die by apoptosis with a death rate constant kma. As Monte Carlo simulations were not very sensitive, we fixed kma to a constant value as previously proposed [19]. The time dependent dosage of Ara-C is denoted by u(t) and is specified by the individual treatment plan and the body surface area (BSA).
Each compartment represents the number of cells per liter. Liter is referred to the peripheral blood, so that the WBC counts from the last compartment coincide with the measured WBC counts from the clinical data. The model also determines the cell numbers of the upstream compartments. These numbers are provided per liter (peripheral blood) as well. If we are interested in the correct number of cells per liter in one of the compartments located in the bone marrow, we need a conversion factor from blood volume to bone marrow volume. We will not discuss this issue further and refer to [38]. Therein, the authors state that the functional haematopoietic marrow volume of about 1.75 L can increase sixfold depending on infection or haemorrhage.
We used a standard two-compartment PK model of high-dose Ara-C, which is administered in the consolidation phase, with zero-order input and linear elimination based on published drug concentration-time data [39]. The elimination rate constant k10, the transfer rate constants k12, k21 and the distribution volume of the central compartment VC were estimated in a previous step and defined as constants for all further computations. A detailed discussion of the PK model is presented in the next section.
The PD effect, e.g. the negative effect of Ara-C on the proliferating cells, linking the PK model, especially the Ara-C concentration, to the myelosuppression model was modelled by a log-linear function E = slope ln(1+ cV x1), using the parameter slope for patient-specific calibration and chemotherapeutical effects and the constant cV for unit consistency (see S3 Table). We also implemented a linear PD function with discouraging results. Additionally, we tested a (sigmoid) Emax model without achieving better model accuracies. The three different PD functions are the commonly used mechanistic models describing pharmacodynamic effects in PK/PD modelling [40].
The function F(x, ktr, γ, B, slope) is a general description of the proliferation rate of xpr and incorporates the PD effect E on the proliferating cells, as discussed in Minami et al. [41], Derendorf et al. [42] and applied, e.g. in Hing et al. [43]. The basic structure of the function F derived in [12], is (1 − E)ktr(B/xma)γ in which the mature cells influence the proliferation rate ktr of xpr with a feedback term (B/xma)γ that leads to higher rates if the number of circulating cells xma is below the baseline WBC count B, and vice versa. It is motivated from studies showing that the proliferation rate can be affected by endogenous growth factors and cytokines [44] and that circulating neutrophil counts and the growth factor G-CSF levels are inversely related [45]. Including this term allows a temporary overshoot of WBC compared with the baseline value B. The proliferation exponent γ indicates the strength or speed of this feedback. The estimation parameters were B, slope, ktr, and γ plus a varying number of additional parameters depending on the initial condition approach, see Table 1. The different initial condition approaches are introduced in the next sections. The WBCs completely recover after cytotoxic therapy (see Fig 2a in [46]) and each cycle is scheduled such that the WBCs should be in homeostasis before treatment start. However, our clinical data indicate that not for all patients the WBCs are yet completely recovered from myelosuppression or they already recovered and overshoot their steady state value before the start of the next treatment cycle due to an overproduction of WBCs. This carry over effect was already mentioned by Kloft et al. [31] but was not considered by them. Thus, we implemented three different strategies to treat initial conditions: I1 assumes a steady state, I2 assumes a steady state only for the proliferating and transient compartments, and I3 penalises deviations from the steady state. I2 and I3 are considered as alternative initial conditions as our clinical data indicate that the steady state assumption after induction phase and between the consolidation cycles may be violated.
We used I1 and ntr = 6 for M1, I1 and ntr = 3 for M2, I1 and ntr = 1 for M3, I2 and ntr = 1 for M4, and I3 and ntr = 1 for M5–M12. Model constants, patient-specific constants and units of model parameters are summarised in S3 Table.
Apart from different PK models which were linked to the myelosuppression model, modifications of the structural model were also proposed [22, 23]. Both models have a more detailed description of the stem cell compartment. The model from Henrich et al. covers a consecutive decrease of the leukocyte’s nadir in the treatment cycles achieved by a prior additional compartment mimicking the slow replication of pluripotent stem cells in the bone marrow. Mangas-Sanjuan et al. models a cell-cycle occurring in the bone marrow compartment covering quiescent cells which do not enter the proliferation process and are not sensitive to the pharmacodynamic effect of the treatment.
Development of PK Model.
We developed a two-compartment PK model for intraveneous high-dose Ara-C infusions and compared it with two previously published PK models [20, 47]. 86 Ara-C concentration measurements (μg/mL = mg/L) from 11 patients were collected from Fig 2 of Kern et al. [39] and presented in Fig 2.
As Fig 2 in [39] was our only source of data, no inter–individual variability analysis could be performed. The patients received high-dose (3 g/m2) Ara-C infusions over 3 hours every 12 hours on days 1, 2, 8 and 9. The measurements we are using were collected at day 1 and 8 and we assume BSA = 1.78 m2. The resulting model with unknown parameters k10, k12, k21 and Vc denoting the elimination rate, distribution rates and the volume of the central compartment, is formulated as
(8)
(9)
We estimated the unknown parameters using a naive pooling approach with exponential error model where
. The naive pooling approach was used as the collected Ara-C measurements could not be assigned to the corresponding patients. The error model is transformed to
and the following parameter estimation problem
(10)
(11)
is solved with a Gauss-Newton algorithm implemented in CasADi [48] with single shooting (CVODES). The estimated parameters and their relative standard deviations are presented in Table 2 together with the parameter values from Solans et al. [47] and Krogh et al. [20].
The two-compartment PK model representing a central and peripheral compartment, see Fig 1, adequately described the concentration-time data and coincides with the derived values for clearance and the elimination rate constant k10 from Table 6 in [39]. We used our derived two-compartment model with the given estimated parameter values in all calculations. We did not use the PK models from Krogh et al. and Solans et al., because they were fitted to low-dose treatment schedules. Although Ara-C is reported to have a linear pharmacokinetics [39], simulations with these PK models did not coincide well with the concentration measurements from high-dose treatment schedules published in [39] as shown in Fig 2.
During the model development process we also tested a K-PD approach [49] which resulted in non-identifiability of the elimination rate constant of the virtual compartment during the individual parameter estimations.
As stated above the PK model was fixed in the modelling process. In one of the simulation studies we analysed the effect of the PK variability on the different modelling hypotheses. We present and discuss the results in S8 Fig.
M6: Modelling a direct killing effect of Ara-C on the proliferating cells.
In the model M6, we chose the proliferation rate as discussed in previous works [21, 33, 42] as F = ktr(B/xma)γ − E. The main difference to all other models is that the PD effect E is directly multiplied with xpr and not with ktr(B/xma)γxpr. Multiplying with xpr can be seen as a direct (killing) impact of Ara-C on the amount of proliferating cells, whereas the more plausible mechanism-based rationale is the induced reduction of the proliferation rate constant ktr used in all models except in M6.
M7–M12: Extending the effects of Ara-C.
The root mean squared error (RMSE) values indicate that model M5 with one transition compartment and initial condition approach I3 (described in the next section) provides the highest accuracy after model personalisation compared to M1–M4.
The indirect effect of Ara-C with an impaired proliferation (M5) is more plausible than a direct killing effect (M6), because Ara-CTP is incorporated into DNA and RNA and impairs cell replication [26]. Therefore, M5 became the reference model for all further analysis. We extended the proliferation rate F(⋅) and/or the transition rate G(⋅) in M5 to capture potential secondary effects of Ara-C. To understand the implications of the extensions, we observe that the proliferation rate F = (1 − E) ktr(B/xma)γ is negative when 1 < E. This is the case for
(12)
This corresponds to more proliferating cells being in the process of apoptosis than being in the process of cell division. It is important that the feedback term (B/xma)γ increases the absolute value of F for B > xma, and decreases it for B < xma. Therefore, an analysis of F always has to consider all four cases related to the signs of 1 − E and of B − xma. Inspired by the log-linear behaviour of the PD effect E, we chose
This monotonously increasing function is applied to different expressions in M5.
In M7 we replaced the transition rate ktr by ktr/S(x1) throughout M5. This results in an Ara-C induced reduction of the transition rate.
In M8 we replaced the complete feedback function F in M5 by F/S(x1). This models an Ara-C induced decreased auto-feedback of the proliferating cells. For high values of x1, i.e. when (12) holds, this results in a decreased killing of proliferative cells. For values x1 > 0 below that boundary, we get a decreased positive proliferation rate.
In M9 we replaced both the complete feedback function F by F/S(x1) and the proliferation exponent γ by γ S(x1). Again, depending on x1 either the killing or the proliferation rate of xpr are decreased by F/S(x1). In addition, the impact depends on whether the WBC count is below or above the baseline: for xma < B we have an increased killing/proliferation rate (B/xma)γ S(x1) > (B/xma)γ and vice versa.
In M10 we replaced the proliferation exponent γ in M5 by γS(x1). This is motivated by the observation that the feedback term with exponent γ is related to the endogenous G-CSF [12]. In contrast to M9 the function F itself is not scaled. Like in M9, the γS(x1) scaling results in an increase of killing/proliferation rates for WBC counts below the baseline, and a decrease else.
In M11 we replaced the quotient B/xma by a comparison between cells in the bone marrow and their baseline value. Based on the statement in [7] and the references therein, we assumed that about 1% of the WBC precursor cells in the bone marrow are in the proliferating compartment xpr, and 99% in the transition compartment xtr. In M12 we combined the extensions from M10 and M11.
The parameter vector contains also initial values, which we discuss next.
Initial conditions of the differential states
The initial values of the PK were chosen to be zero, x1(t0) = x2(t0) = 0, due to the fact that the considered time horizons start before administration of chemotherapy. Further, it is known that previous Ara-C treatments have no impact on the pharmacokinetics of subsequent treatments [39] which is supported by simulation studies showing that the values of x1(t) and x2(t) are below 10−6 after 16.35 days of the 1 Ara-C infusion. The remaining initial conditions were chosen using one of the following three strategies.
Initial condition approach I1.
The WBC count xma(t0) was set to the long term WBC baseline (steady state) count B. With this particular choice all feedback terms simplify to
Assuming x1 = 0, also E = 0, S(x1) = 1, and hence F = G = ktr for all models at time t0, which simplifies the analysis. As initial conditions for the cell counts in the bone marrow we chose the bone marrow baseline cell count
,
(13)
(14)
which guarantees that inserting (13) into (1) leads to a right hand side of zero.
The advantage of this approach is the identifiability of the estimation problem, as no additional degrees of freedom in the estimation problem need to be introduced for the initial conditions. However, simulations showed that both the assumption of xma(t0) = B and the steady state assumption were typically violated at the beginning of a new consolidation cycle.
Initial condition approach I2.
One additional parameter B0 was introduced and estimated, as suggested by Nock [34]. It was used for the initialisation as
(15)
together with (14). The time derivative (1) at time t0 is given by
(16)
(17)
(18)
which is not zero for B0 ≠ B. The advantage of this approach is that also increasing or decreasing WBC counts at t0 can be captured, depending on the sign of (18).
Initial condition approach I3.
The initial conditions , and, as in I2, also xma(t0) = B0 were introduced as additional estimation parameters. As this leads to unidentifiability of the estimation problem, a term penalising deviations from (13),
with α = 1/2500 was added to the objective function of the least squares estimation problem. The regularisation parameter α was chosen with respect to the tradeoff between identifiability of xpr(t0) and xtr,i(t0) and the violation degree of the steady state assumption. Larger values of α resulted in similar parameter estimates compared to the initial condition approach I1, as more attention was drawn to the penalising terms guaranteeing a solution close to the steady state assumption after parameter estimation. Smaller values of α weakened the steady state assumption but also increased the uncertainty of parameter estimates for xpr(t0) and xtr,i(t0) resulting in large standard deviations. A good tradeoff was achieved with α = 1/2500. For more information about the regularisation approach see [50] and the references therein. Note that B−B0 was not penalised. This approach is the most flexible with respect to the possibly transient initial dynamics resulting, e.g. from previous treatments.
Clinical data (high density WBC counts) & personalisation
AML patients who had received induction therapy (commonly defined as anthracycline- and Ara-C-based 7+3 regimen [2]) resulting in complete remission and who did not receive granulocyte-colony stimulating factors (G-CSF) during the post-remission consolidation therapy were eligible for data analysis. We focused on patients who did not receive growth factor support, as such effects were not yet accounted for in our mathematical models. Almost daily WBC counts from 42 consolidation Ara-C cycles (CCs) of 23 AML patients (median 62 years, 14 male, mostly de novo AML (19/23), mostly AML FAB-M2 (9/19), mostly intermediate cytogenetic risk (12/20)) from 2008 to 2015 were analysed from clinical charts provided by the Department of Hematology and Oncology, Magdeburg University Hospital, Magdeburg, Germany. The data were retrospectively collected and pseudonymised from records of the clinical routine, and interventions were not performed for this work. For this reason no patients’ agreements were required. This study was approved by the Ethics Committee of the Medical Faculty of the University of Magdeburg (Optimal control of clinically relevant cancer chemotherapy schedules in patients with acute leukaemia—with special emphasis on neutropenia, MARTINA; approval number 124/15). All clinical procedures were performed in accordance with the general ethical principles outlined in the Declaration of Helsinki. The CCs were partitioned in one, two, and three consecutive CCs from nine, nine, and five patients, respectively. Four different schedules D135, d135, D123, or D12, in which the numbers correspond to treatment days 1, 2, 3, and 5, respectively, d to intermediate-dose Ara-C (i.e. 1 g/m2 per BSA twice a day over three hours) and D to high-dose Ara-C (i.e. 3 g/m2 twice a day), were administered 23, 15, two, and two times. Patient PD123 (62 years, male) received two cycles of D123. Patient PD12 (64 years, female) received two cycles of D12. The 21 other patients received 1-3 D135 cycles (median 57 years, 8 male, 4 female) or d135 cycles (median 68 years, 5 male, 4 female).
The cycle- and patientwise longitudinal WBC count measurements are published as NONMEM compatible comma separated values files in the supporting information (S1 and S2 Files). Both datasets contain the columns ID TIME DV CMT AMT RATE DUR MDV and EVID. ID serves as an identifier for the appropriate cycle, respectively patient. In the cyclewise dataset each of the 42 CCs has its individual ID meaning that each cycle is treated independently although several cycles belong to the same patient. In the patientwise dataset all cycles belonging to the same patient are assigned to one ID. TIME [day] either specifies the measurement times of WBC counts or the starting times of the Ara-C administrations. DV [G/L] is the dependent variable, in our case individual WBC count measurements. The column CMT specifies the compartment in which a dosing or observation event occurs. AMT [mg], RATE and DUR [day] define the Ara-C schedules. As every administration is an infusion lasting 3 hours, the entries of DUR and RATE are 0.125 [day] and -2. AMT defines the BSA-adjusted amount [mg] of Ara-C. 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. For more information about NONMEM datasets we refer the interested readers to [51].
We used all 42 CCs to personalise our mathematical models M1–M12 performing point estimations (individual approach). The point estimates were used to analyse the different modelling assumptions. Additionally we personalised the most relevant models M3, M10 (with I1), the model from Henrich et al. [22], the model from Mangas-Sanjuan et al. [23] and the model from Stiehl et al. [8] applying nonlinear mixed-effects modelling (population approach). The population approach is used on the one hand to qualitatively confirm our proposed model variation based on the set of population parameters. On the other hand we wanted to compare our set of population parameter values with recently published models and give a reason why the Friberg model serves as our basic model and not recently published models which are similar to the Friberg model, but with several extensions. A comment on the use of initial value approach I1 instead of I3 for the population approach is given in the S1 Appendix. Once the model parameters have particular values, the model is called personalised model (PM).
Prediction & cross-validation
The PMs were then used to predict (simulate) and cross-validate WBC counts for the last CC of 14 patients for whom at least two consecutive CCs are available. Additionally, we calculated predicted trec values from our 42 PMs applying D123 and D135 schedules and compared the descriptive statistics with published average trec values from a subset of data (367 CCs of 208 AML patients, no G-CSF support) of the AMLSG 07-04 trial in which the schedules D123 and D135 after 7+3 regimen were analysed [46]. The published AMLSG 07-04 [46] trial does not provide WBC counts to obtain new PMs, therefore we used the median of observed trec values for D123 and D135 Ara-C schedules. In the interest of a fair comparison (i.e., to avoid comparison with the value 0) we excluded five (d135:1 and D135:4) out of 42 PMs for which at least one out of the 42 predictions (M1–M12 with either D123 or D135) resulted in no WBC counts below the threshold value. This could occur as we personalised the models for a specific treatment plan, e.g. D135. Afterwards we applied a different treatment plan to the PMs, i.e D123, which may have resulted in a reduced cytotoxic effect. Not each out of the 42 predictions resulted in a nadir value below 1 G/L. Further, we predicted trec values for two Ara-C schedules in which a constant administration of Ara-C throughout days 1-5, with either 100 mg/m2 per day or 400 mg/m2 per day was given. These schedules, together with D135, have been clinically analysed for 1088 AML patients (median 52, 568 male) by Mayer et al. [3], and the superiority of D135 with respect to disease-free survival rates and remaining in continuous complete remission after four years has been shown but no trec values were reported. Finally, we analysed the effect of the inter-individual PK variability on the trec values derived by the models M3 and M10 (with I1). We applied schedules D123 and D135 with fixed population parameter values for B, slope, ktr, and γ and performed 500 simulations each with randomly chosen values from constructed inter-individual variability (IIV) for the PK parameters clearance CL and central volume VC.
All experiments were performed to analyse the 12 proposed models with respect to WBC count and trec predictability.
Schedule timing
After verifying the predictability performance of the PMs, we performed a simulation study in which we demonstrated a further possible application of the PMs in planning the start of consecutive CCs. We analysed the impact of the treatment timing on the individual nadir values. For each of the 14 patients, for whom at least two consecutive CCs were available, the nadir of the last CC was compared to 20 simulated nadirs. These nadirs resulted from simulations using the patient’s PMs (second row of Table 3) in which the timing of the last CC was varied daily with the maximal starting variation of 10 days earlier or later.
Results
Accuracy of PMs with fixed Ara-C schedule
Table 3 shows statistics about the accuracies of the PMs describing the clinical data, for a pure estimation (using all available WBC counts to personalise the model) and for a cross-validation (using all but the last CC for personalisation).
The accuracies depend strongly on the number of compartments and initial condition strategy (M1–M5), but do not differ much with respect to modelling assumptions of possible effects of Ara-C considered in M6–M12. These values were even better when the standard schedule D135 was applied in the estimated and predicted cycles. Regarding the root mean squared errors for M1-M5 the results implied that one transition compartment and initial approach I3 were the best choice for the structural model and hence served as a starting point to analyse different pharmacodynamic effects of Ara-C. As mentioned in the previous section, the number of transition compartments determines the MMT of the differentiating progenitor cells. Comparing the MMTs resulting from the population approach for M1-M3 we achieved a slight decrease from 154 h to 144 h to 128 h by using the corrected formula MMT = n/ktr[17] instead of the original formula MMT = (n + 1)/ktr [12]. During the administration of cytostatic drugs it is known that the cells are encouraged to rapidly differentiate such that a MMT of 128 h is reasonable. Furthermore, the MMT value from one transition compartment is closest to a previously published corrected MMT value of 106.4 [17, 30]. Studies with healthy volunteers reported MMTs of 153.6 h and 165.6 h [17]. But these values are difficult to compare as chemotherapy can speed up proliferation and differentiation. With this knowledge and the accuracy values of Table 3 we decided to fix the number of transition compartments to one. The original MMT formula from [12] would have resulted in 180 h, 193 h and 256 h. By using M5 as the reference model and analysing different hypotheses of Ara-C’s PD effect in M6-M12, all models could describe the clinical data equally well. Goodness-of-fit plots in Fig 3a and 3b and S1 Fig visually support the good match between model predictions and measured WBC counts (respectively observed vs. calculated trec values) around the nadir and a wider spread of large WBC counts. To analyse the reliability of the PMs to predict the WBC dynamics in subsequent CCs, Fig 3c and 3d indicate the involved model uncertainty from parameter uncertainty by means of Monte Carlo simulations. The model uncertainty was derived from 1000 randomly chosen parameter sets sampled from the variance-covariance matrix resulting from the individual parameter estimation problem (10) (for more information see S1 Appendix). The information from one CC and no available prior knowledge leads to a high uncertainty. The uncertainty reduces when more WBC counts are present, and the prediction accuracy for consecutive CCs and myelosuppression increases. Examining the accuracy of the PMs for each patient separately, the WBC counts around the nadir are explained well by all models for fixed Ara-C schedules (either D135 or D123), as shown in Fig 4a–4d for two exemplary patients and in S2 and S3 Figs for the other 12 patients with at least two consecutive CCs.
(a) Goodness-of-fit plot for M10. Shown are measured versus calculated white blood cell (WBC) counts. Models were personalised using complete data sets of one to three cycles from 23 patients. The measured counts around the nadir coincide well (RMSE = 0.740) with the calculated WBC counts. (b) As (a), but cross-validated: WBC counts from the last cycle of patients were not used for personalisation, but compared to predictions (RMSE = 0.927). The plot shows cross-validated WBC counts from the last cycle in red, others in blue. The plots are prototypical for M1–M12. (c) PMs based on M10 and either personalisation with WBC counts from one or from all three cycles. 1000 Monte Carlo simulations after personalisation with WBC counts from one cycle were used to indicate the propagated probability density function. (d) As (c), but using WBC counts from the first two cycles for personalisation. More measurements lead to higher prediction accuracy. The uncertainty tube tightens and the predicted trajectory gets closer to the solution that used all available WBC counts.
Patient with three D135 cycles (left) and patient PD123 with two D123 cycles (right), as indicated on the x-axis. The PMs exemplify reproducability (first row), predictability (second row) and simulation of a different schedule in prediction than estimation (third row). (a) Reproducability: all 12 PMs based on M1–M12 are able to explain the measured WBC counts. (b) As in (a), all PMs explain the measured WBC counts well, particularly around the nadirs. (c) Cross-validated prediction: all PMs explain the WBC counts well, also in the predicted third cycle. (d) As in (c), here with a slightly too slow predicted recovery time in the second cycle for all models. (e) Varied Ara-C schedule: prediction of D123 in the third cycle for a PM based on two D135 cycles shows faster WBC recovery for M9, M10, and M12. (f) Prediction of D135 in the second cycle for a PM based on one D123 cycle shows slower WBC recovery times for M9, M10, and M12.
Regarding the estimated parameter values, we only determine a slight change of the estimated fixed-effects parameter values for B, ktr and slope, the inter-individual variability for all four parameters and the residual error between models M3 and M10 whereas the estimated fixed-effects parameter value for slope significantly decreases when a second PD term is introduced (see).
Accuracy of PMs with altered Ara-C schedule
As we were interested to differentiate between the distinctive model hypothesis, we applied different chemotherapy schedules to the PMs and analysed their dynamical behaviour, especially the WBC recovery. Fig 4e and 4f show two cases where D135 was used for personalisation and D123 for prediction (and vice versa). Here, M9, M10, and M12 have a faster (slower) haematological recovery for D123 (D135). All three models assume that the proliferation speed γ depends on the Ara-C concentration. This modelling assumption is visualised in a different way for M5, M10 and M12 in S4, S5 and S6 Figs.
It is shown that the negative proliferation rate F of M10 and M12 compared to M5 has an altered dynamical behaviour during chemotherapy due to the increased γ value from the PD effect achieving a faster WBC recovery for D123 schedules. For this accelerated feedback relationship between WBCs and G-CSF biological interpretations are given in the next section.
In conclusion, the comparison of WBC recovery times between D123 and D135 treatments is a suitable criterium for model discrimination.
The next study was performed to compare the calculated trec values from M1-M12 with clinically collected values to figure out which of the models coincided with clinical findings. We used 444 PMs (using M1–M12 and clinical data from 37 cycles with schedules D135, d135, D123 and D12 from section Clinical Data & Personalisation) to predict the outcome of D135 and D123 schedules. The median values of the predicted trec were compared to the values from a subset of data (108 with D135 and 259 with D123 schedules) from the AMLSG 07-04 trial [46]. M9, M10, and M12 resulted in roughly 4 days faster trec for D123 compared to D135, similar to the clinical result from the literature and in contrast to the 1 day difference of M5 (compare Table 4). The individual results have been qualitatively confirmed by the predicted trec values from the population approach (see). The models from Henrich et al. and Mangas-Sanjuan et al. were not further considered, as both models simplified to the Friberg model after parameter estimation. For the model from Henrich et al., the estimated population parameter value ftr was 0.96, supporting the visual assessment that the patients’ nadirs are not decreasing during the CCs. The estimated parameter values kcycle = 0.0009 and Fprol = 0.941 of the model from Mangas-Sanjuan et al. yielded a non-existing stem cell cycle. A possible reason for the non-identifiability of the parameters might be the limited schedule variation. The authors state that a vast variation of schedules has to be available for parameter identification [23]. The model from Stiehl et al. provided the highest model accuracy with respect to the final objective function value, but exhibited disagreeing WBC recovery times and large relative standard errors (see S5 Table). Therefore we did not consider this model in our further studies. The simulation study analysing the effect of the PK variability on the resulting recovery times of schedules D123 and D135 for models M3 and M10 (with I1) revealed that model M10 was more sensitive to different high-dose Ara-C treatment schedules compared to model M3 despite the high inter-individual PK variability. This was verified in S8 Fig presenting boxplots of 500 simulated trec values for both models and schedules with constructed IIV on the PK.
Schedule timing
Similar to previous simulation studies dealing with varying and shortening cycle duration and finding the optimal number and timing of G-CSF administrations to reduce myelosuppression [52–54], we analysed the impact of different treatment starts of the last CCs with respect to obtained nadir values. A comparison to the clinically observed nadir values indicated a large potential for clinical improvement, i.e., a higher nadir value due to a different treatment timing (see Fig 5a). Fig 5b exemplarily shows the WBC dynamics for different treatment timings. Earlier (later) starts resulted in sequentially higher (lower) nadir values.
(a) Simulation study in which 20 simulated nadirs were compared with the true nadir of the last CC for the 14 patients who have more than one CC. The simulated nadirs were computed by using the patient’s PM (second row of Table 3) and varying the start of the last CC daily with the maximal starting variation of 10 days earlier or later. (b) Exemplary variation of the CC start for one patient. An earlier (later) start results in a larger (lower) nadir.
Discussion
High-density WBC counts from 23 AML patients were collected and used to personalise 12 mathematical models and analyse their prediction accuracy with respect to different modelling hypotheses and treatment schedules. The high prediction accuracies of the PMs, especially around the nadir, confirm previous claims [4, 24] that the general approach of in-silico studies can be used for clinical decision support. As clinical decision support we understand tools which help physicians to monitor and predict WBC dynamics and the duration and grade of myelosuppression. In combination with clinical expertise on the impact of schedules on relapse probabilities and their small scope determining the start of the next cycle due to subjective experience and the patients fitness, this might have an important clinical impact via altered treatment schedules which might eventually result in decreased depth and duration of myelosuppression.
Current drawbacks are the high model uncertainty, if insufficient information is available. This makes precise and reliable predictions difficult (compare Fig 3c and 3d). Furthermore, the lack of leukemic cell dynamics and the validation of the PMs based on one specific chemotherapy schedule might lead to not appropriate models concerning an optimisation of Ara-C dosage.
Comparing the estimated parameter values with published values, the estimated baseline value for the WBCs is within the normal human WBC range of 4.5−10 G/L but reduced by 2 G/L compared to published baseline values for the Friberg model being in the range of 7 to 7.8 G/L [12]. The mean maturation time of 128 h for M3 using the corrected formula [17] is reasonable and fits into the range of previously published values [55]. The estimated γ values are roughly two to three times higher compared to published values for the reason that we only use one transition compartment. The γ values for M2, containing three transition compartments, are in the same range then published values. The decrease of the slope parameter value from one to two PD effects occurs as the effect of Ara-C is distributed on two different sites of action. During the parameter estimation we observed for some CCs correlations (> 0.9) between γ and ktr and between γ and slope. But these correlations had no influence on the parameter identifiability. Further, it was shown that under certain assumptions, which we fulfilled (kprol = ktr), the Friberg model is structurally globally identifiable [55]. In future studies we propose to use global design measures from [56] to provide treatment schedules reducing global parameter sensitivity and undesired parameter correlation.
We showed that an analysis based on a fixed chemotherapy schedule cannot discriminate between different modelling hypotheses. The agglomerative nature of the mathematical models leads to a choice of model parameters that is not only personalised to the patient, but also to the applied schedule. Therefore, we used different schedules for personalisation and prediction to overcome this problem and to allow discrimination of the models. This approach allowed us to distinguish between the modelling hypotheses implemented in models M5–M12 and enabled us to find the suitable model assumption considered in M9, M10, and M12. In our opinion this procedure should be routinely applied, preferably using high density WBC counts for different schedules in the same patients. As an alternative to such a tedious clinical study we suggest to use average trec values as a discrimination criterion for competing models.
Comparing the trec values for D123 and D135 treatments from the PMs with our clinical data and the AMLSG 07-04 trial in Table 4 and S2 Table implies, that model M10 (based on Ockam’s razor in comparison to M9 and M12) is the best candidate among M1–M12 for future work on the simulation and optimisation of intermediate to high-dose Ara-C treatment schedules.
The 1 day shift in trec values between our clinical data and the AMLSG 07-04 trial can be explained by the age difference between patients in our clinical data (median 62 and 57 years for D123 and D135, respectively) and the subdata of the AMLSG 07-04 trial (median of all patients in the trial 49 years) and a related statistical analysis: Jaramillo et al. [46] found in a multivariable analysis a significantly longer WBC recovery for older patients (hazard ratio of a 10-year age difference, 0.89; P = 0.001) [46] and a significantly shorter WBC recovery for patients receiving D123 compared to the reference group D135 (hazard ratio, 1.94; P < 0.0001) [46] which coincides well with our findings.
Regarding the PK model, no published compartment models for high-dose Ara-C are available. Comparing our model with published low-dose Ara-C models, we showed that the published models do not reach the measured maximum Ara-C concentrations from high-dose schedules (see Fig 2) so that we rely on our derived model. As we logarithmised the collected Ara-C concentrations, lower values became more important during parameter estimation such that our fitted PK model slightly underpredicts the highest Ara-C concentrations (see Fig 2). Nevertheless, our model achieves higher values compared to the models from [20] and [47] providing a more reasonable PK behaviour of Ara-C. We estimated a smaller central volume leading to a reduced clearance activity derived from an almost equivalently estimated elimination rate constant value (Table 2). The distribution rate constants differ by a factor of 2 to 2.5 and the peripheral volume by a factor of almost 10. Future PK studies for high-dose Ara-C can be used for model verification or updating our model parameters. In a simulation study we analysed the influence of constructed PK variability on the WBC recovery time for models M5 and M10. We showed that model M10 is more sensitive to varied PK dynamics and reflects clinical findings more accurately, i.e. that the standard and dense treatment plans result in significantly different WBC recovery times. M5 was not able to match the clinical results. A critical part of the study is the constructed IIV. In Krogh et al. [20] IIV was analysed for low-dose Ara-C schedules. We used the published values as exemplary IIV values within our simulation study. Obviously, the results should be treated with care as IIV is related to the underlying study, treatment, model and population and thus cannot be applied to other studies in general. However, the IIV impacted only a small part of our sensitivity analysis and we do not expect qualitative changes for updated IIV values.
Summarising, we extended the gold-standard model for myelosuppression [12] to the most important component in consolidation therapy [2, 25], Ara-C, and showed that one modelling assumption was important for a faster WBC count recovery for D123 schedules. In models M9, M10, and M12 we assumed that the Ara-C concentration has a direct impact on the proliferation speed. As stated above, such a modelling assumption has an agglomerative nature and the underlying physiological processes are still unknown.
We speculate that the increased number of cell deaths following chemotherapy might play a role. Cell deaths induce phagocytosis and macrophage activation, which in turn might increase G-CSF secretion and hence proliferation speed [57] resulting in a positive circulus vitiosus. This might also be one explanation for the observation that G-CSF levels in neutropenia are increased after chemotherapy compared to primary diseases [58]. Independent from the underlying phyiological process, the dense treatment schedule D123 profits more than the standard plan D135 from the induced dynamics with respect to WBC recovery times as discussed in S4 and S5 Figs.
Future G-CSF concentration measurements for AML patients during consolidation cycles of D123 and D135 treatments and a comparison of our extended models with Quartino’s [30] integrated G-CSF-myelosuppression model or more sophisticated models from quantitative systems pharmacology [54] may shed light on these speculations.
Supporting information
S1 Table. Terminology and potentially confusing synonyms.
Expressions that are also used in the manuscript are in italic.
https://doi.org/10.1371/journal.pone.0204540.s005
(PDF)
S2 Table. Comparison of model predictions for low-dose treatment schedules.
As in Table 4, predicted nadir values for different treatment schedules are shown, based on underlying mathematical models M1–M12. Shown are the values of median, standard deviation, minimum and maximum (in brackets) for two low-dose schedules. Both assume a continuous infusion throughout days 1 to 5, with either 100 mg/m2 or 400 mg/m2 Ara-C per day. No clinical observations are available to compare these predictions, but they give additional insight on the possibility to discriminate models M1–M12 and a general trend showing that for 100 mg/m2 per day despite of M7-M9 almost all nadir values are above 1 G/L. The nadir values for the low-dose infusion with 400 mg/m2 Ara-C per day are in the same range compared to the results of the high-dose schedules (Two further personalised cycles were excluded because for some models no recovery after chemotherapy was observed). The simulated nadirs above 1 G/L for the low-dose schedule (100 mg/m2) reflect the lower toxic effects represented by required hospitalisation due to fever and neutropenia and platelet transfusions compared to the low-dose (400 mg/m2) and high-dose schedules explored in [3]. As M7-M9 are not able to reflect the lower toxic effects through higher nadir values, the simulation study serves as an indicator that the secondary effect of Ara-C may not be an Ara-C induced reduction of the transition rate.
https://doi.org/10.1371/journal.pone.0204540.s006
(PDF)
S3 Table. Model constants, patient-specific constants, and units of model parameters.
The values were used to obtain personalised mathematical models. The constants were determined from published data [39] and applied to all patients. To shorten notation we also used . The patient-specific infusion times and dosages that define a treatment schedule were modified for simulation and optimisation of different schedules. The range shows minimum and maximum values of all considered data in the clinical study.
https://doi.org/10.1371/journal.pone.0204540.s007
(PDF)
S4 Table. Objectives (final objective function values from FOCEi method (OBJ), population predicted 
and 
values), parameter and coefficient of variation (CV) estimates with relative standard errors (RSE) from nonlinear mixed-effects modelling of models M3 and M10 with initial condition approach I1 (with I1).
https://doi.org/10.1371/journal.pone.0204540.s008
(PDF)
S5 Table. Objectives (final objective function values from FOCEi method (OBJ), population predicted 
and 
values), parameter and coefficient of variation (CV) estimates with relative standard errors (RSE) from nonlinear mixed-effects modelling of model from Stiehl et al. [8].
https://doi.org/10.1371/journal.pone.0204540.s009
(PDF)
S1 Fig. Goodness-of-fit plot for all but three (because of WBC counts greater 1) measured and calculated trec values for M10 after model personalisation for each consolidation cycle.
The 39 measured trec values are slightly higher due to the coarser measurement grid.
https://doi.org/10.1371/journal.pone.0204540.s010
(PDF)
S2 Fig. Cross-validation of predicted white blood cell (WBC) counts from personalised models (PMs) M1-M12 and measured WBC counts for six patients treated with D135.
All but one cycle were used for personalisation and the last cycle for cross-validation. For patients (a)-(e) the PMs can predict the WBC count decrease after Ara-C administration in the last cycle where models M10 and M12 have a slower WBC recovery than M5. For patient (f) the WBC recovery from the PMs starts to early compared to the measured WBC counts.
https://doi.org/10.1371/journal.pone.0204540.s011
(PDF)
S3 Fig. Cross-validation of predicted white blood cell (WBC) counts from personalised models (PMs) M1-M12 and measured WBC counts for five patients (a)-(e) treated with d135 and one patient (f) treated with D12.
The PMs provide good predictions for patient (a) and (f) but show mismatches in recovery times and nadir values for patients (b)-(e).
https://doi.org/10.1371/journal.pone.0204540.s012
(PDF)
S4 Fig. Comparing personalised mathematical models (PMs) M5 and M10 for D123 and D135 schedules (exemplary patient I).
https://doi.org/10.1371/journal.pone.0204540.s013
(PDF)
S5 Fig. Comparing personalised mathematical models (PMs) M5 and M10 for D123 and D135 schedules (exemplary patient II).
https://doi.org/10.1371/journal.pone.0204540.s014
(PDF)
S6 Fig. Comparing personalised mathematical models (PMs) M10 and M12 for D123 and D135 schedules (exemplary patient I).
https://doi.org/10.1371/journal.pone.0204540.s015
(PDF)
S7 Fig. Visual predictive checks (VPCs), derived by 1000 simulations, for leukocytes [G/L] versus time [days] starting with the first measurement before dosing for model M3 (a) and M10 (with I1) (b).
Blue circles are the measured WBC counts of 23 AML patients described in section Clinical Data & Personalisation. One measurement was taken at timepoint 88.98 [days] with the value 7.18 [G/L] which is not shown in the VPCs. Red lines show the median (solid) and 5th and 95th percentiles (dashed) of measurements. The shaded areas represent the 95% confidence intervals around the 5th (blue), 50th (red) and 95th (blue) simulated percentiles of the model predictions. Regarding the VPCs, model M3 and M10 have an almost equivalent prediction accuracy. The 50% percentiles of measurements and model predictions perfectly overlap, thus supporting our individually based results from Table 3. The same applies to the start of the 5% and 95% percentiles until the nadir. After the nadir the 5% and 95% percentiles of the model predictions recover slightly faster/slower compared to the measurements. At day 30 the percentiles of measurements and model predictions coincide again.
https://doi.org/10.1371/journal.pone.0204540.s016
(PDF)
S8 Fig. Simulation study analysing the sensitivity of model M5 and M10 on inter-individual PK variability (IIV) when schedules D123 and D135 are applied.
(a) As Fig 2, but with 500 simulations of our fitted two-compartment PK model with IIV on the clearance and the central volume. (b) Recovery times (trec) from 500 simulations each of models M3 and M10 (with I1) applying schedules D123 and D135 with inter-individual variability given as coefficient of variation (CV) on PK parameters clearance (45%) and central volume (70%). Red lines within the boxes are the medians, the upper and lower box limits are the first (Q1) and third quartiles (Q3) of the data. The lower whiskers will extend to the first trec values greater than the first quartiles minus the 1.5-times the interquartile ranges (IQR) (Q1 − 1.5 * IQR). Equivalently, the upper whiskers will extend to last trec values less than Q3 + 1.5 * IQR. Beyond the whiskers, data are considered as outliers and are plotted as individual points (+). The simulation study revealed that model M10 was more sensitive to different high-dose Ara-C treatment schedules compared to model M3 despite the high inter-individual PK variability.
https://doi.org/10.1371/journal.pone.0204540.s017
(PDF)
Acknowledgments
The authors wish to thank Franziska Kluwe, Graduate Research Training Program PharMetrX and Department of Clinical Pharmacy and Biochemistry, Institute of Pharmacy, Freie Universitaet Berlin, for her support in the pharmacokinetics and nonlinear mixed-effects modelling analysis.
References
- 1. Malka R, Wolach B, Gavrieli R, Shochat E, Rom-Kedar V. Evidence for bistable bacteria-neutrophil interaction and its clinical implications. J Clin Invest. 2012;122(8):3002–3011. pmid:22820292
- 2. Döhner H, Estey E, Grimwade D, Amadori S, Appelbaum FR, Buchner T, et al. Diagnosis and management of AML in adults: 2017 ELN recommendations from an international expert panel. Blood. 2017;129(4):424–447. pmid:27895058
- 3. Mayer RJ, Davis RB, Schiffer CA, Berg DT, Powell BL, Schulman P, et al. Intensive postremission chemotherapy in adults with acute myeloid leukemia. Cancer and Leukemia Group B. N Engl J Med. 1994;331(14):896–903. pmid:8078551
- 4. Wallin JE, Friberg LE, Karlsson MO. Model-based neutrophil-guided dose adaptation in chemotherapy: evaluation of predicted outcome with different types and amounts of information. Basic Clin Pharmacol Toxicol. 2010;106(3):234–242. pmid:20050841
- 5. Jost F, Sager S, Le TTT. A feedback optimal control algorithm with optimal measurement time points. Processes. 2017;5(10):1–19.
- 6. Manesso E, Teles J, Bryder D, Peterson C. Dynamical modelling of haematopoiesis: an integrated view over the system in homeostasis and under perturbation. J R Soc Interface. 2013;10(80):20120817. pmid:23256190
- 7. Marciniak-Czochra A, Stiehl T, Ho AD, Jäger W, Wagner W. Modeling of asymmetric cell division in hematopoietic stem cells—regulation of self-renewal is essential for efficient repopulation. Stem Cells Dev. 2009;18(3):377–386. pmid:18752377
- 8. Stiehl T, Ho AD, Marciniak-Czochra A. Mathematical modeling of the impact of cytokine response of acute myeloid leukemia cells on patient prognosis. Sci Rep. 2018;8(1):2809. pmid:29434256
- 9. Østby I, Rusten LS, Kvalheim G, Grøttum P. A mathematical model for reconstitution of granulopoiesis after high dose chemotherapy with autologous stem cell transplantation. J Math Biol. 2003;47(2):101–136. pmid:12883857
- 10. Engel C, Scholz M, Loeffler M. A computational model of human granulopoiesis to simulate the hematotoxic effects of multicycle polychemotherapy. Blood. 2004;104(8):2323–2331. pmid:15226178
- 11. Craig M, Humphries AR, Mackey MC. A mathematical model of granulopoiesis incorporating the negative feedback dynamics and kinetics of G-CSF/neutrophil binding and internalization. Bull Math Biol. 2016;78(12):2304–2357. pmid:27324993
- 12. Friberg LE, Henningsson A, Maas H, Nguyen L, Karlsson MO. Model of chemotherapy-induced myelosuppression with parameter consistency across drugs. J Clin Oncol. 2002;20:4713–4721. pmid:12488418
- 13. Clairambault J. Modelling physiological and pharmacological control on cell proliferation to optimise cancer treatments. Math Model Nat Phenom. 2009;4(3):12–67.
- 14. Wang W, Stiehl T, Raffel S, Hoang VT, Hoffmann I, Poisa-Beiro L, et al. Reduced hematopoietic stem cell frequency predicts outcome in acute myeloid leukemia. Haematologica. 2017;102(9):1567–1577. pmid:28550184
- 15. Craig M. Towards quantitative systems pharmacology models of chemotherapy-induced neutropenia. CPT Pharmacometrics Syst Pharmacol. 2017. pmid:28418603
- 16. Fornari C, O’Connor LO, Yates JW, Cheung SA, Jodrell DI, Mettetal JT, et al. Understanding hematological toxicities using mathematical modeling. Clin Pharmacol Ther. 2018. pmid:29604045
- 17. De Souza DC, Craig M, Cassidy T, Li J, Nekka F, Bélair J, et al. Transit and lifespan in neutrophil production: implications for drug intervention. J Pharmacokinet Pharmacodyn. 2018;45(1):59–77.
- 18. Tetschke M, Lilienthal P, Pottgiesser T, Fischer T, Schalk E, Sager S. Mathematical modeling of RBC count dynamics after blood loss. Processes. 2018;6(9):157–185.
- 19. Quartino A, Friberg L, Karlsson M. A simultaneous analysis of the time-course of leukocytes and neutrophils following docetaxel administration using a semi-mechanistic myelosuppression model. Invest New Drugs. 2012;30:833–845. pmid:21153753
- 20. Krogh-Madsen M, Bender B, Jensen MK, Nielsen OJ, Friberg LE, Honoré PH. Population pharmacokinetics of cytarabine, etoposide, and daunorubicin in the treatment for acute myeloid leukemia. Cancer chemother and Pharmacol. 2012;69(5):1155–1163.
- 21. Pefani E, Panoskaltsis N, Mantalaris A, Georgiadis MC, Pistikopoulos EN. Chemotherapy drug scheduling for the induction treatment of patients with acute myeloid leukemia. IEEE Trans Biomed Eng. 2014;61(7):2049–2056. pmid:24686224
- 22. Henrich A, Joerger M, Kraff S, Jaehde U, Huisinga W, Kloft C, et al. Semimechanistic bone marrow exhaustion pharmacokinetic/pharmacodynamic model for chemotherapy-induced cumulative neutropenia. J Pharmacol Exp Ther. 2017;362(2):347–358. pmid:28600397
- 23. Mangas-Sanjuan V, Buil-Bruna N, Garrido MJ, Soto E, Trocóniz IF. Semimechanistic cell-cycle type–based pharmacokinetic/pharmacodynamic model of chemotherapy-induced neutropenic effects of diflomotecan under different dosing schedules. J Pharmacol Exp Ther. 2015;354(1):55–64. pmid:25948593
- 24. Netterberg I, Nielsen EI, Friberg LE, Karlsson MO. Model-based prediction of myelosuppression and recovery based on frequent neutrophil monitoring. Cancer Chemother Pharmacol. 2017; p. 1–11.
- 25. Löwenberg B, Pabst T, Vellenga E, van Putten W, Schouten HC, Graux C, et al. Cytarabine dose for acute myeloid leukemia. N Engl J Med. 2011;364(11):1027–36. pmid:21410371
- 26.
Perry MC, Doll DC, Freter CE. Chemotherapy source book. Wolters Kluwer Health/Lippincott Williams & Wilkins, Philadelphia; 2012.
- 27. Plunkett W, Liliemark JO, Adams TM, Nowak B, Estey E, Kantarjian H, et al. Saturation of 1-β-D-arabinofuranosylcytosine 5-triphosphate accumulation in leukemia cells during high-dose 1-β-D-arabinofuranosylcytosine therapy. Cancer Res. 1987;47(11):3005–3011. pmid:3471322
- 28. Li Z, Guo JR, Chen QQ, Wang CY, Zhang WJ, Yao MC, et al. Exploring the antitumor mechanism of high-dose cytarabine through the metabolic perturbations of ribonucleotide and deoxyribonucleotide in human promyelocytic leukemia HL-60 cells. Molecules. 2017;22(3):499.
- 29. Herold N, Rudd SG, Ljungblad L, Sanjiv K, Myrberg IH, Paulin CB, et al. Targeting SAMHD1 with the Vpx protein to improve cytarabine therapy for hematological malignancies. J Nat Med. 2017;23(2):256–263.
- 30. Quartino A, Karlsson M, Lindman H, Friberg L. Characterization of endogenous G-CSF and the inverse correlation to chemotherapy-induced neutropenia in patients with breast cancer using population modeling. Pharm Res. 2014;31(12):3390–3403. pmid:24919931
- 31. Kloft C, Wallin J, Henningsson A, Chatelut E, Karlsson MO. Population pharmacokinetic-pharmacodynamic model for neutropenia with patient subgroup identification: comparison across anticancer drugs. Clin Cancer Res. 2006;12(18):5481–5490. pmid:17000683
- 32.
Quartino A. Pharmacometric models for improved prediction of myelosuppression and treatment response in oncology [PhD thesis]. Acta Universitatis Upsaliensis; 2011.
- 33. Pefani E, Panoskaltsis N, Mantalaris A, Georgiadis MC, Pistikopoulos EN. Design of optimal patient-specific chemotherapy protocols for the treatment of acute myeloid leukemia (AML). Comput Chem Eng. 2013;57:187–195.
- 34.
Nock V. Pharmacometric modelling of processes in the haematopoietic system and blood [PhD thesis]. Freie Universität Berlin; 2013.
- 35.
Rinke K, Jost F, Findeisen R, Fischer T, Bartsch R, Schalk E, et al. Parameter estimation for leukocyte dynamics after chemotherapy. In: Proceedings of the Foundations of Systems Biology in Engineering (FOSBE) Conference. vol. 49. Magdeburg, Germany. Elsevier; 2016. p. 44–49.
- 36. Savic RM, Jonker DM, Kerbusch T, Karlsson MO. Implementation of a transit compartment model for describing drug absorption in pharmacokinetic studies. J Pharmacokinet Pharmacodyn. 2007;34(5):711–726. pmid:17653836
- 37. Brown G, Hughes PJ, Michell RH, Rolink AG, Ceredig R. The sequential determination model of hematopoiesis. Trends Immunol. 2007;28(10):442–448. pmid:17825625
- 38. Batinić D, Marusić M, Pavletić Z, Bogdanić V, Uzarević B, Nemet D, et al. Relationship between differing volumes of bone marrow aspirates and their cellular composition. Bone Marrow Transplant. 1990;6(2):103–107. pmid:2207448
- 39. Kern W, Schleyer E, Unterhalt M, Wörmann B, Büchner T, Hiddemann W. High antileukemic activity of sequential high dose cytosine arabinoside and mitoxantrone in patients with refractory acute leukemias. Cancer. 1997;79(1):59–68. pmid:8988727
- 40. Felmlee M, Morris M, Mager D. Mechanism-based pharmacodynamic modeling. Methods Mol Biol. 2012;929:583–600. pmid:23007443
- 41. Minami H, Sasaki Y, Saijo N, Ohtsu T, Fujii H, Igarashi T, et al. Indirect-response model for the time course of leukopenia with anticancer drugs. Clin Pharmacol Ther. 1998;64(5):511–521. pmid:9834043
- 42. Derendorf H, Meibohm B. Modeling of pharmacokinetic/pharmacodynamic (PK/PD) relationships: concepts and perspectives. Pharm Res. 1999;16(2):176–185. pmid:10100300
- 43. Hing J, Perez-Ruixo J, Stuyckens K, Soto-Matos A, Lopez-Lazaro L, Zannikos P. Mechanism-based pharmacokinetic/pharmacodynamic meta-analysis of trabectedin (ET-743, Yondelis) induced neutropenia. Clin Pharmacol Ther. 2008;83(1):130–143. pmid:17597713
- 44.
Hoffbrand AV, Higgs DR, Keeling DM, Mehta AB. Postgraduate haematology. John Wiley & Sons; 2016.
- 45. Takatani H, Soda H, Fukuda M, Watanabe M, Kinoshita A, Nakamura T, et al. Levels of recombinant human granulocyte colony-stimulating factor in serum are inversely correlated with circulating neutrophil counts. Antimicrob Agents Chemother. 1996;40(4):988–991. pmid:8849265
- 46. Jaramillo S, Benner A, Krauter J, Martin H, Kindler T, Bentz M, et al. Condensed versus standard schedule of high-dose cytarabine consolidation therapy with pegfilgrastim growth factor support in acute myeloid leukemia. Blood Cancer J. 2017;7(5):e564. pmid:28548643
- 47. Solans BP, Fleury A, Freiwald M, Fritsch H, Haug K, Trocóniz IF. Population pharmacokinetics of volasertib administered in patients with acute myeloid leukaemia as a single agent or in combination with cytarabine. Clin Pharmacokinet. 2018;57(3):379–392.
- 48.
Andersson J. A General-Purpose Software Framework for Dynamic Optimization [PhD thesis]. Arenberg Doctoral School, KU Leuven; 2013.
- 49. Jacqmin P, Snoeck E, Van Schaick E, Gieschke R, Pillai P, Steimer JL, et al. Modelling response time profiles in the absence of drug concentrations: definition and performance evaluation of the K–PD model. J Pharmacokinet Pharmacodyn. 2007;34(1):57–85. pmid:17051439
- 50.
Bock HG, Körkel S, Schlöder JP. In: Bock HG, Carraro T, Jäger W, Körkel S, Rannacher R, Schlöder JP, editors. Parameter estimation and optimum experimental design for differential equation models. Heidelberg: Springer; 2013. p. 1–30.
- 51.
Owen JS and Fiedler-Kelly J. Introduction to population pharmacokinetic/pharmacodynamic analysis with nonlinear mixed effects models. John Wiley & Sons; 2014.
- 52. Scholz M, Engel C, Loeffler M. Modelling human granulopoiesis under poly-chemotherapy with G-CSF support. J Math Biol. 2005;50(4):397–439. pmid:15614553
- 53. Vainas O, Ariad S, Amir O, Mermershtain W, Vainstein V, Kleiman M, et al. Personalising docetaxel and G-CSF schedules in cancer patients by a clinically validated computational model. Br J Cancer. 2012;107(5):814. pmid:22814580
- 54. Craig M, Humphries AR, Nekka F, Bélair J, Li J, Mackey MC. Neutrophil dynamics during concurrent chemotherapy and G-CSF administration: mathematical modelling guides dose optimisation to minimise neutropenia. J Theor Biol. 2015;385:77–89. pmid:26343861
- 55. Evans ND, Cheung SYA, Yates JWT. Structural identifiability for mathematical pharmacology: models of myelosuppression. J Pharmacokinet Pharmacodyn. 2018;45(1):79–90. pmid:29396780
- 56. Schenkendorf R, Xie X, Rehbein M, Scholl S, Krewer U. The impact of global sensitivities and design measures in model-based optimal experimental design. Processes. 2018;6(4):27.
- 57. Mendoza J, Caceres J, Santiago E, Mora L, Sanchez L, Corona T, et al. Evidence that G-CSF is a fibroblast growth factor that induces granulocytes to increase phagocytosis and to present a mature morphology, and that macrophages secrete 45-kd molecules with these activities as well as with G-CSF-like activity. Exp Hematol. 1990;18(8):903–910. pmid:1696904
- 58. Kavgaci H, Ozdemir F, Aydin F, Yavuz A, Yavuz M. Endogenous granulocyte colony-stimulating factor (G-CSF) levels in chemotherapy-induced neutropenia and in neutropenia related with primary diseases. J Exp Clin Canc Res. 2002;21(4):475–479.