In-silico comparison of two induction regimens (7 + 3 vs 7 + 3 plus additional bone marrow evaluation) in acute myeloid leukemia treatment

Acute myeloid leukemia (AML) is a rare malignant disease of the blood cell formation and is the most common acute leukemia among adults leading to most death events caused by leukemias [1]. In particular, AML is constituted of genetically different hematopoietic neoplasms that collectively stem from various multistep mutations affecting the myeloid cell line resulting in accumulating neoplastic precursor cells [2]. The intrinsic origin of AML is a small subset of leukemic stem cells (LSC) leading to a self-perpetuating proliferation of clonal progenitor cells also referred to as blasts [3]. Fast growing number of inoperative and undifferentiated blasts induce a disorder of normal hematopoiesis located in the bone marrow with further systemic implications in blood and other tissues [1]. As part of a clonal evolution in one patient, genetically different AML clones exist, develop and are specifically responsible for diagnosis or potential relapse because of a presumed selection by chemotherapy [4].

This complex pathogenesis and additional resistance mechanisms lead to various treatment strategies and respective various patient outcomes [5, 6]. Independently from new more tailored treatment approaches (e.g. CAR-T cells), an established but relatively unspecific combination chemotherapy of cytarabine and anthracycline sets still the standard aiming for a first clinical remission during induction therapy [2, 5].

A widely used therapeutic approach is the 7 + 3 regimen (starting with seven days cytarabine and supplementary for the first three days anthracycline). In clinical practice variation of this regimen exists differing e.g. in dose and/or schedule [5, 7]. Different 7 + 3 regimens are preferred depending on respective region, e.g. an evaluation process with potential re-induction in the United States compared to a preferred double induction in Europe [8].

Our scientific objectives were to compare different intensities of two 7 + 3 chemotherapy regimens using a mathematical model, which characterizes the dynamics of AML using ordinary differential equations. We primarily intend to increase the efficiency of this known induction therapy by detecting more disease-specific treatment conditions. A single-induction 7 + 3 regimen was compared with a 7 + 3 regimen plus an additional bone marrow (BM) evaluation on day 14 and/or 21 after treatment start with a potential second induction cycle.

In total, we analyzed about ten thousand different intensity combinations, which is more than in tangible experiments (in vivo or in vitro) or clinical studies are feasible [9–11].

To evaluate each scenario, we computed the time from treatment start to complete remission (CR) and the following CR duration as two essential clinical parameters that enable a rational comparison [12].

For our analysis, we extended the published two-compartment AML model by Stiehl et al. [13]. The two compartments represent hematopoietic stem cells (HSC) within the bone marrow (first compartment) which can differentiate by cell division into non-proliferating (differentiated) cells (second compartment). Healthy and non-healthy (leukemic stem cells, LSC) cells are modelled separately and differ in their parameter values. The model is able to explain the dynamics of cell population abundance adequately [13–16]. Normal HSC and pathological LSC are represented by a set of two ordinary differential equations. Cellular abundance (in cells/kg body weight) at day t is denoted by c₁(t) for HSC, c₂(t) for healthy differentiated cells, l₁(t) for LSC and l₂(t) for non-proliferating leukemic cells, respectively. HSC, LSC and non-proliferating leukemic cells are considered to reside in the bone marrow, whereas healthy differentiated cells belong to the blood stream. This model assumption corresponds with prior results by Stiehl [13], who showed that a model extension by a bone marrow exit did not lead to qualitative changes of the selected cell properties via chemotherapy. The proliferation rates p^c (HSC) and p^l (LSC) and the self-renewal rates, (a^c for HSC, a^l for LSC), are used to characterize healthy hematopoiesis and leukemia growth kinetics, respectively. In the model the term of self-renewal represent this self-sustaining fraction as a proportion (0–1). We considered proliferation rates in the range 0–2 and self-renewal rates between 0 and 1.

Depletion (e.g. apoptosis or migration into, for our model, negligible states) of non-proliferating cells is modelled as constant death rate \( {d}_2^c \) and\( {d}_2^l \). Resident cells (here c₁, l₁, and l₂) will also deplete, if bone marrow cell number exceed a density threshold value, i.e. the physiological equilibrium value of bone marrow cell count. The function d(x(t)) is an additional death rate describing the fraction of bone marrow cells dying because of overcrowding. A feedback regulation denoted by s(t) is integrated to represent cellular communication. Affecting self-renewal, feedback regulation leads to the result that increasing number of differentiated healthy cells causes a reduced number of HSC and LSC (and vice versa).

Intensive induction chemotherapy of AML contains a combination of two or more chemotherapeutics applied in a specific therapy regimen [5, 7]. We extended the model to implement the 7 + 3 scheme based on 7 days cytarabine and 3 days anthracycline treatment. Cytarabine acts as an antimetabolic agent and attacks primarily on cells during their synthesis phase (S-phase) by inhibiting the DNA-polymerase α [17, 18]. The modelled chemotherapy mechanism (k_cyt · p^c · c₁(t)) can be considered as a cytarabine-like chemotherapy acting on proliferating cells. Anthracycline affects proliferating and non-proliferating cells via various mechanisms (e.g. inhibition of topoisomerase II or free radical generation) [19]. A second chemotherapy mechanism has been introduced acting on non-proliferating cells. We assume that the effect of anthracycline on mitotic cells is limited to the proliferation phase. This is justified by experimental observation that anthracycline’s toxicity on mitotic cells is mainly due to active proliferative state [20–22]. The extended Stiehl model [13] is defined by

The model represents an intermediate state between Model 1 and Model 2 from [13]. It is known that the bulk of leukemic cells express granulocyte colony-stimulating factor (G-CSF) receptors [23]. G-CSF is the main mediator for hematopoietic feedback regulation and can stimulate as well leukemic cells [24]. G-CSF-feedback regulation is mainly directed by, not fully understood, transcriptional signaling processes with STAT3/SOCS proteins [25–27]. A relevant amount of AML subtypes shows a significant dysregulation of STAT3/SOCS-related pathways [28–31]. Thus, we assume for a wider part of AML no negative G-CSF feedback regulation by leukemic cells. As an implementation the leukemic cell self-renewal depends on the feedback s(t) (eq. (6)), but the feedback s(t) does not depend on the leukemic cell count in our model.

We implemented the model in the statistical software R [32]. Numerical solutions for the ordinary differential equations were calculated using the R package ‘deSolve’ [33]. A more detailed description, including the parametrization of the model, is given in Additional file 1. The R syntax is given in Additional file 2.

To analyze the model, we performed the following steps

Simulation of two combination therapy regimens

Two chemotherapy regimens of AML have been modelled and analyzed (Fig. 1).

All simulations start with a small amount (1 cell / kg) in the L1 compartment and trajectories are observed. The initial state of the simulations is given in Table 1. When the blasts percentage hits the diagnostic threshold, therapy is triggered.

Compartment	Initial state
c1(t = 0)a	2 × 109 cells/kg
c2(t = 0)a	3.9 × 109 cells/kg
l1(t = 0)	1 cell/kg
l2(t = 0)	0 cells/kg

^aCompartment c₁ and c₂ start in the steady state of the healthy systems (without leukemic cells or therapy)

The “standard arm” resembles a single 7 + 3 induction course [38] with a cytarabine-like chemotherapy for 7 days and anthracycline-like chemotherapy for 3 days. No further therapy is applied and the course of the leukemia is observed until the end of the simulation (2000 days, 5000 days for the slow pace leukemia).

A variation of the 7 + 3 chemotherapy regimen that is based on the 2017 guidelines of the National Comprehensive Cancer Network (NCCN) [1] has also been implemented. Analogously to classical 7 + 3 the prescribed combination therapy will be applied, if blast fraction exceeds 20%. A first evaluation of therapy success by assessing blast fraction is made on day 14 after treatment start. A second induction therapy will be performed, if blast fraction ≥5.5%. However, second induction therapy will be intensity-reduced (5 days of cytarabine-like chemotherapy and 2 days of anthracycline-like chemotherapy (5 + 2)). If blast fraction is < 5.5%, a second evaluation will be made on day 21 after treatment start. In that case, a blast fraction ≥5% also leads to a second therapy course (5 + 2). For blast fractions < 5% no further treatment will be applied. Within one simulation run intensity values (k_cyt, k_anthra) will be preset and fixed for 7 + 3 and 5 + 2.

To identify a realistic range of chemotherapy intensities (k_cyt and k_anthra, in unit cell deaths per day) we sampled therapy combinations starting with no therapy (k_cyt= k_anthra=0) and increased intensities until we reached the area of overtreatment (i.e. no complete remission can be achieved and all compartments are completely depleted by intensive therapy). By intention, we included the mono-treatment scenarios. While cytarabine mono-treatment is common practice in the pre-phase, consolidation and treatment of the elderly [1] anthracycline mono-treatment is uncommon in clinical practice. Finally, we simulated all scenarios between intensity values between 0 and 10 (step size 0.1). We could show for the selected fast and intermediate leukemia that for k_cyt>8.8 (mono-therapy) in the standard arm and k_cyt> 5.8 in the evaluation regimen leads, regardless of k_anthra, leads to overtreatment (complete depletion of all compartments). Under-treatment exists for low k_cytand k_anthra rates that lead to no CR. As internal validation, we also measured absolute reduction of leukemic cells (l₁ and l₂) per chemotherapy combination on day 29 after treatment start (compared to leukemic cell abundance at diagnosis). For the effective regions (no under/over-treatment) we observed reductions between 10⁵and 10⁹ cells, resembling realistic values [39].

An exemplary simulation of a fast leukemia under a specific treatment dose is shown in Fig. 2, which depicts the cell number trajectories and blast percentages over time.

We modelled three leukemias by varying proliferation and self-renewal rates leading from initial mutation to three different times of diagnosis. Results from AML pathogenesis research show that growth properties (such as proliferation rate) affect different survival outcomes [41, 42]. This finding is also supported by mathematical modeling results [40]. Quantifying leukemias by growth kinetics (e.g. leukemic pace as time-to-diagnosis) is relevant. Various mutations were identified and specific cytogenetics are linked to different patient outcomes [43–45]. However, there is no information about the time to diagnosis for specific leukemia types and how chemotherapy influences dynamics of hematopoiesis over time. For example, bone marrow examinations are performed for diagnosis and 7 to 10 days after induction chemotherapy [5]. In between, no continuous data is collected to reduce patient burden. When modern cytometric techniques get more available this gap can be closed. These procedures, especially if based on bone marrow samples, cannot be used for continuous monitoring [46]. Mathematical models are able to bridge this diagnostic gap and can reveal suspected useful therapeutic implications by its dynamical approach.

Here, we presented results based on a homogenous approach (one type of leukemic cells) to characterize the dynamic behavior of certain AML subtypes. Nevertheless, it is known that AML is a multiclonal disease [4, 35]. We modelled a combination chemotherapy attacking this one leukemia and summarized the simulated therapy combinations. In reality, a mixture of leukemic cells with different properties is observed at diagnosis and a clonal evolution leading to relapse can be demonstrated [4, 47]. Despite of this variety, in most cases a dominant clone induces AML onset [48]. We only focused on induction therapy’s impact on this dominant clone. Post-remission therapy like consolidation chemotherapy is intentionally not considered and is subject of future research. As a consequence, our comparison of two induction regimens aims exclusively at an improvement of CR achievement without considering effects on potential relapse in the course of clonal evolution. Our presented results must be assessed in the specific model context and a direct comparison with usual clinical outcome or endpoint parameters, which are based on patient populations, cannot be made instantly. In our focused consideration, relapses by different clones are not possible, in order that quasi-infinite CRs are obtained by induction treatment by complete depletion of a single leukemic clone. In multiclonal models relapses are expected.

With respect to realistic therapy concepts several limitations occur. The leukemia’s sensibility to chemotherapy usually is influenced by drug resistance mechanisms [49, 50]. The therapy effectiveness (k_cyt, k_anthra) can be considered as combinations of therapy intensity (dose) and the leukemia’s resistance (only affected by proliferation rate and cell numbers) to the therapy, e.g. due to its specific genetic characteristics. At present, the model cannot simulate an AML-type specific resistance. Future model extensions will aim for data-derived proliferation and self-renewal parameters (representing specific genotypes) and respective resistance mechanisms. Besides, in clinical practice, chemotherapy intensity is applied in units of mg/m² adapted to body surface to take account of side effects [5]. Currently, we cannot compare model parameters to clinical therapy intensities directly. The model is at least able to review different doses qualitatively on a course-grained scale (high vs low). Another future prospect will be to link the model parameters in a pharmacodynamics model to therapy doses, to e.g. replay studies intensifying induction by dose increase [51–53]. With respect to the selected leukemia parameter combinations via leukemic pace our selection might bias our conclusions if leukemic proliferation and self-renewal parameters influences our selected outcome measures significantly. Additional simulations (data not shown) indicate that time to CR is not influenced significantly by leukemia parameters. For lower intensity therapy combinations CR duration is related to the leukemias proliferation rate, while self-renewal has no effect. There seems to exist a threshold proliferation rate. Below this threshold, only very short CRs can be observed. While we focused on therapy intensity in this publication, a more elaborate analysis of the interplay between leukemia characteristics and therapy outcomes will be necessary and should be investigated in future research.

To assess the modelled chemotherapy values compared to realistic used intensities, we used the established criterion of 3 log₁₀ cytoreduction, which is at the minimum required for a reduction of leukemic cells under 5% in bone marrow [39]. In addition, a reduction of transcription products of more than 3 log₁₀ levels is also used as a prognostic factor in minimal residual disease (MRD) monitoring after induction therapy [54–56]. Therefore, in relation to an adequate chemotherapy intensity a log₁₀ reduction ≥ 3 of leukemic cells can be considered as a predictor for treatment success. All chemotherapy intensities leading to CR feature a leukemic cell reduction > 3 log₁₀ levels. In fact, in the model reductions often exceed this criterion. Referring to minimal detection level of minimal residual disease (MRD) with sensitivities between 10^− 4 and 10^− 5 [46], our model provides a starting basis for further and novel MRD investigations by demonstrating cell trajectories (with precise blast percentages) over time (Fig. 2). Common medical diagnostics cannot enable a comparable continuous view.

In the next steps, patient data like healthy steady state stem and progenitor cell numbers must be integrated into the model. In future, exact cell number analysis of a patient will be challenging, especially the knowledge transfer from mice models to manageable in-vivo analysis [43]. At the same time, a determination of self-renewal must be derived from this patient individual data [13, 45]. Availability of personalized parameter values lead to a further specialized model would be able to translate full spectrum of genetic change into specific values of proliferation and self-renewal [57]. Subsequently, each individual AML as an own genetic entity and effects of therapy could be highly effectively modelled and assessed. Regarding this, the categorization of AML types by current classification systems (e.g. ELN, MRC, WHO, FAB) considering cytomorphological, genetical and immunological properties is complex because of known heterogeneity of AML [5]. These classifications especially regard rather static properties like mutations or the immunophenotype. Derived aggregations lead to risk groups that comprises similar patient outcomes, but they do not precisely describe how fast the AML proliferates, nor which resistance mechanisms exist and which consecutive dynamical impact on the hematopoietic system is generated. Hereby, the mathematical model provides a functional perspective, which allows a more individual analysis of the AML pathogenesis and therapy effects.

Varying therapy regimens are used worldwide, that mostly differ in time and duration of chemotherapy administering [8]. The presented treatment model is also suitable for a planned comparison of different double induction concepts like TAD-HAM vs. S-HAM [7]. In addition to our evaluation regimen investigation, we will analyze different evaluation time points to find out, when bone marrow assessment should be optimally done. Indeed, evaluation timing is still an open issue and our dynamical model could be helpful to provide additional valuable insights [58–60].

Risk-stratification concepts for AML therapy are currently gaining in importance and are particularly established in post remission therapy [61–63]. Hereby, treatment strategy is especially influenced by estimated outcome [64]. We could observe that personalized treatment intensities in induction chemotherapy produce relevant advantages (lower minimum necessary chemotherapy intensities). Therefore, prospective stratification concepts might also include inherent properties of AML.

Complementary to randomized controlled trials (RCT), modelling can be understood as a tool that can add a holistic point of view to classical reductionist medicine [65]. Moreover, clinical relevant modelling consists of a hypothesis-driven research, that connects in-silico experimental results with established experimental facts in a scientific interacting cycle [66, 67]. Concerning an effective personalized medicine in AML treatment, we are convinced that this interdisciplinary approach will be inevitable and offers great potential. At present, our model can derive clinical relevant conclusions despite the prescribed limitations, because our integral dynamical approach enables new insights in AML-hematopoiesis and optimal chemotherapy effect relating to a certain AML types.

Our results suggest, that the “7 + 3” regimen results in CR more often. In addition, more therapy combinations result in quasi-infinite CR. This holds for the fast and the intermediate paced leukemia (in the model the slow paced leukemia is not treatable to achieve CR). The results support the current scientific view that “7 + 3” regimen is a standard of care independent from existing diverse regimen variations that are applied in study groups all over the world [7]. Nevertheless, a more extensive evaluation and comparison with many more established therapy schemata is necessary.

We assume that genetic heterogeneity of every leukemic clone determines unique characteristics that require corresponding unique therapy concepts. This assumption is based on significantly different survival outcomes that are strictly dependent on specific genetic constitution [5]. Induction therapy is not routinely adapted to genetic disposition and patients are treated with standardized induction doses only adapted to body surface [5]. Concepts of higher doses or adding a third agent were implemented in several randomized studies, but the comparison proved difficult and dose increases were not precisely adapted to individual patient [7]. Regarding this, our model suggests that a whole spectrum of effective (that means CR as well as prolonged CR) chemotherapy intensities from relatively low to high exists and “7 + 3” regimen offers a larger effective spectrum that is respectively adjusted to characteristics of considered clone type. On the hypothesis that these insights hold true concerning real life, “7 + 3” regimen could be lead to higher cure probability of a standardized dose applied to heterogeneous diseases because of more existing effective dose combinations for each AML type. Nevertheless, a more substantial model based comparison of “7 + 3” and other therapy regimens is lacking.

In our simplifying model, the essential criterion for an optimal clinical outcome is ultimately duration of CR. As a result, we initially consider different intensity combinations as equivalent, as far as achieved CR duration is similar. However, side effects of chemotherapy are modelled via a cytotoxic effect on blood cells (indeed model is calibrated to number of neutrophil granulocytes as most frequent leukocytes [13]). Thus, lower and higher intensities leading to same CR duration only differ in absolute cytoreduction without affecting the defined outcome. Other significant side effects concerning for example the gastrointestinal tract or residual blood system (or anyway related infections) are not factored in. For that reason, we consider lower intensities leading to same CR duration as superior. There is growing evidence that a relevant part of progress in AML outcome is due to improved supportive therapy [68–70]. Finding lowest as possible and at the same time effective therapy intensities seems to be eminently important [71]. To that effect, we observe for both regimens that spectrum of effective therapy intensities leading to CR as well as persistent CR only differs in the region of lowest intensities. We can derive from this model results that lower and still effective induction therapy intensities may exist depending on different AML clone type. In the concrete case of our model, a fast leukemia can be treated efficiently with lower intensities than the intermediate leukemia. We record that the evaluation regimen enables lowest intensities leading to prolonged CR (i.e. most efficient treatment) for both selected parameter combinations.

It should be noted that evaluation regimen provides these lower effective intensities especially for the intermediate paced leukemia, because hereof largest reduction of essential intensity (in relative comparison to “7 + 3”) is obtained. We emphasize that in our model evaluation regimen does not offer more effective therapy combinations but the most efficient regarding optimal outcome and respective minimum intensity. Hereby, in our model the evaluation approach is particularly worthwhile to minimize therapy intensities and consecutive side effects with regard to medium-fast proliferating leukemia. The main potential of the evaluation approach is presumably present for AML of poorer risk categories. This insight complies inherently with current results that discusses necessity of bone marrow assessment and recommends a more individualized decision of evaluation [72].

An early response to the first induction cycle is a known prognostic factor, but its impact on evaluation process stays unclear [58, 73]. Our results show that therapy combinations enabling fastest CR (in our model quasi-instantly after chemotherapy) constitute the majority of achievable CR independent of administered regimen and leukemia pace. In consideration of intensity and optimal outcome, we detected therapy combinations leading to CR within 20 to 45 days. Concerning this efficient subset, in general, “7 + 3” enables faster CR and in addition, the intermedium paced leukemia takes longer time to CR. Persistent CR under minimal dose (the most efficient situation) is obtained by evaluation regimen regardless of whether other therapy combinations lead to even faster CR.

Published literature proposes that CR should be reached as fast as possible [74–76]. According to our modeling results, therapy combinations leading to CR within the first 10 days does not always enable a long-lasting CR. A therapy intensification can lead to longer times until onset of CR, but then lead also to a more stable CR. In summary, in our model fast as possible CR achievement is not inevitably optimal.

Enabling persistent CR with examined minimal dose, cytarabine-like monotherapy turns out to be optimal for both regimens. This model result conflicts with clinical reality, which ascribes relevant importance to a combination chemotherapy for decades [38]. Our model considered one type of leukemic cells per patient (see model limitations above). We know from recent studies, that in one patient several subtypes exist and combination chemotherapy leads to a selection process [13, 47, 57]. Some subtypes are more resistant to this chemotherapy for several reasons (e.g. a lower proliferation rate). In our homogenous model, we do not examine the prescribed selection process, because we focus on the treatment effect regarding one specific AML. The success of classical “7 + 3” cannot be considered as resounding because of existing numerous AML subgroups with poor survival outcomes [5]. Therefore, multi-layered approaches of targeted therapy (e.g. immunotherapy or pathway inhibition) are recently under investigation without major breakthrough until now [7].