Predictive nonlinear modeling of malignant myelopoiesis and tyrosine kinase inhibitor therapy

The primary challenge in developing mathematical models of normal and CML hematopoiesis is sorting through the combinatorial explosion of models that occurs when cell–cell signaling interactions are taken into account. Consider the model hematopoietic system shown in Figure 1A, which accounts for hematopoietic stem (HSC; S), multipotent progenitor (MPP; P), and two types of postmitotic, terminally differentiated cells—myeloid (TD_m) and lymphoid (TD_l). The HSC self-renew with fraction (e.g., probability) p₀ or differentiate with fraction 1-p₀. That is, the fraction of HSC that remain as HSC after division is p₀. The MPPs self-renew with fraction p₁ and differentiate into either lymphoid or myeloid cells with fractions q₁ and 1-p₁-q₁, respectively. The HSC and MPPs divide with rates η₁ and η₂ and the myeloid and lymphoid cells die at rates d_m and d_l, respectively. The ODEs that govern the dynamics of the cells are given in ‘Methods.’ Assuming that there is either positive or negative regulation of the self-renewal and differentiation probabilities and division rates of any cell type from any other cell type results in 59,049 models, counting each combination of regulated signaling as a separate model.

Figure 1

Download asset Open asset

To select the most physiologically accurate models, we first filtered the models using an automated approach (DSA) developed by Savageau and co-workers (Savageau et al., 2009; Fasani and Savageau, 2010; Lomnitz and Savageau, 2016) that enables models to be distinguished based on their range of qualitatively distinct behaviors, without relying on knowledge of specific values of the parameters. This method relies on identifying boundaries in parameter space that separate qualitative behaviors of a particular model, which is much more efficient than searching for model behaviors directly. The boundaries can be approximated from a sequence of inequalities that identify regions where one term on the right-hand side of each ODE (e.g., the rate of change) dominates all others in the sources (positive terms) and another dominates the sinks (negative terms). This is known as a dominant subsystem (S-system) of the model. The number of S-systems in each model depends on the number of combinations of positive and negative terms in the rates of change. If the equilibria of the S-systems, which are determined analytically, are not self-consistent (e.g., consistent with the assumed dominance of terms reflected in the inequalities) or the equilibria are not stable, then the S-system is rejected. If all the S-systems of a particular model are rejected, then that model is removed from further consideration. Models with at least one self-consistent and stable S-system are viable candidates for further analysis. DSA can be easily automated to make the analysis of very large numbers of equations feasible. Details are provided in ‘Methods’ and ‘Appendix 1’ (Section 1). The result of this procedure is the elimination of all but the four model classes shown in Figure 1B, which require negative regulation of the stem cell self-renewal fraction but differ by where this regulation arises. The models within the classes share at least one S-system and have common qualitative behaviors. The differences between models in a class lie in whether or not there is positive, negative, or no regulation on the rest of the parameters from any of the cell types. This reduces the number of possible models to 26,244.

Previous work has implicated several feedback mechanisms active in both normal and malignant hematopoiesis. Interleukin-6 (IL-6) is produced by differentiated myeloid cells and acts to bias MPPs toward a myeloid fate (Reynaud et al., 2011; Welner et al., 2015). Such negative feedback circuits, known as fate control, have been shown to provide an effective strategy for robust control of cell proliferation and reduction of oscillations in branched lineages (Buzi et al., 2015). The chemokine CCL3 (also known as macrophage inhibitory protein α [MIP-1α]), produced in BM by basophilic myeloid progenitors (Baba et al., 2016), acts to inhibit the proliferation and self-renewal of normal HSC (Broxmeyer et al., 1989; Staversky et al., 2018), but CML HSC are relatively resistant to its action (Eaves et al., 1993, Baba et al., 2013). In hypothesizing these regulatory networks, we arrived at a single model class as shown in Figure 1C. In this class, there are 729 model candidates, which differ only in how the HSC and MPP cell division rates and the MPP self-renewal fraction are regulated. These above results suggest that IL-6 is a candidate feedback factor expressed in the myeloid compartment (TD_m) with the ability to negatively regulate the fraction q₁ of MPPs that differentiate into lymphoid cells. CCL3 is a candidate factor mediating negative feedback from the MPP population onto HSC self-renewal. To further constrain the remaining models, we performed cell biological experiments in mice to glean information about cell–cell interactions by separately perturbing the stem cell and myeloid cell compartments.

As described in ‘Methods,’ healthy C57BL6/J (B6) mice were treated with low-dose (50 cGy) ionizing radiation, previously shown to be selectively toxic to HSC in the BM (Stewart et al., 1998). The BM stem/progenitor compartment was analyzed by flow cytometry in untreated mice, and on days 1, 3, and 7 post-irradiation, using the gating strategy in Figure 2A. These time points and the number of mice analyzed at each time point were informed by a Bayesian hierarchical framework for optimal experimental design of mathematical models of hematopoiesis (Lomeli et al., 2021). In particular, the Bayesian framework suggests combining early time points (soon after radiation was applied) with late time points because the early time points provide more information about division rates, while the late time points provide more information about the feedback parameters. One day after treatment, we observed an acute approximately twofold decrease in the relative size of the HSC compartment in the irradiated mice (Figure 2B), accompanied by approximately threefold increase in proliferation rates for both HSC and MPPs (Figure 2C). There was no significant change in MPP population, however, and the system returned to equilibrium by day 7. These results suggest that the HSC population exerts negative feedback on their own division rate (η₁) and inhibits the division of MPPs through a negative feedforward loop on η₂.

Figure 2

Download asset Open asset

B6 mice were treated with the anti-granulocyte antibody RB68C5 (50 μg), and their BM was analyzed 1 d after treatment (see ‘Methods’). This treatment resulted in an ~20% decrease in mature BM myeloid cells, as measured by CD11b expression (Figure 2D and E), and was accompanied by a concomitant increase in the size of the phenotypic MPP compartment (Figure 2E) and a decrease in the HSC compartment (Figure 2E). These results suggest that there is a negative feedback loop from the myeloid cells onto the MPP self-renewal fraction p₁.

Taking all these results into consideration, we arrive at the feedback-feedforward model shown in Figure 2F. The negative feedback loops shown in blue correspond to those suggested by previous experimental data, while the negative feedback and feedforward loops in red are suggested by the cell depletion experiments presented here. See ‘Methods’ for a detailed description of the corresponding ODEs. Although these validation data were derived from mice, we hypothesize that similar cell–cell signaling occurs in humans.

To determine biologically relevant parameters for the feedback-feedforward model in Figure 2F, a grid-search algorithm was employed. The full ODE model is given in ‘Methods’ and Appendix 1 (Section 2). The 12 model parameters (proliferation and death rates, self-renewal and branching fractions, feedback/feedforward gains) were sampled using a random uniform distribution for each parameter. See ‘Methods,’ Appendix 1 (Section 3), and Appendix 1—tables 2 and 3 for details and a full parameter list. Once parameter values were chosen, the model was simulated for long times. If a parameter set resulted in steady state values consistent with the range of values previously reported for a dynamic human hematopoiesis model (Manesso et al., 2013), that parameter set was accepted. Out of ~10⁶ possible parameter combinations, a total of 1493 parameter sets were accepted (Appendix 1—figure 4). We further restricted the candidate parameter sets by considering only those with sufficiently large feedforward gains on the MPP division rate (γ₅ > 0.01) in order to focus on the novel feedforward dynamics. This reduced the number of eligible parameter sets to 563, and their distributions are shown in Appendix 1—figure 5. Each of these parameter sets can be thought of as representing the ‘normal’ condition of a virtual patient by having different individual parameters, for example, due to genetic, epigenetic or environment factors, that nevertheless result in a ‘normal’ homeostatic hematopoietic system. The different parameter sets thus model a range of variability across individual CML patients. The values of the parameters used are given in Appendix 1, Section 3.

DSA can be used to determine qualitative model behaviors and how sensitive the model is to perturbations of key parameters. Here, we focused on the feedback gains γ₁ and γ₃ on the HSC and MPP self-renewal probabilities, respectively (see Appendix 1, Section 1.3 for details, and for sensitivity analyses for other parameters, see Figure 3—figure supplement 1). As indicated in Figure 3A, DSA identifies four regions (design space) in the γ₁ and γ₃ plane which the dynamics are governed by different S-systems. Using a parameter set in each design space region (indicated by white dots) as a base value, we performed a parameter sweep in which we vary γ₁ and γ₃ in a range within 0.9–1.1 times the magnitude of their original values. In Figure 3B, the evolution of each of the cell populations is shown, starting from an initial condition in which there are only a small number of HSC. The different graphs correspond to the parameter sets (Appendix 1—table 4) in the four regions of the design space although the dynamics are shown for the full ODE solutions. The solid curves denote results from the original (white dot) parameter set, and the shading denotes the range of behaviors when the parameters are varied. The black and blue curves correspond to the HSC and MPPs, respectively, while the dark-green and light-green correspond to the terminally differentiated myeloid and lymphoid cells. While the system tends to equilibrium for all parameter combinations, the approach to equilibrium is different. The dynamics in regions i and ii are monotonic while those in regions iii and iv are not (e.g., the equilibria in regions i and ii are stable nodes, while those in regions iii and iv are stable spirals). Further, the larger the γ₁, the faster the approach toward equilibrium. The cell numbers and proportions in each design space region are different as well. In regions i and ii, the HSC dominate while in regions iii and iv the differentiated myeloid cells dominate the population. Further, the number of cells in regions i and iii is larger than those in regions ii and iv. The equilibrium cell populations in region iii correspond more closely to the physiological populations identified by Manesso et al., 2013. Figure 3—figure supplement 2 depicts the effective parameters in region iii as it develops toward the physiological steady state.

Figure 3 with 2 supplements see all

Download asset Open asset

We next investigated the sensitivity of the model to perturbations about the equilibrium cell population. In Figure 3C, we present the results obtained by reducing the number of terminally differentiated myeloid cells from their equilibrium value by 10% (dot-dashed), 50% (dashed), and 90% (solid) and with parameters from design space region iii (Appendix 1—table 4). By initially depleting the myeloid cells, which is similar to the experiment in Figure 2D and E, the hematopoietic system is shifted away from its steady state. While the presence of the negative feedback loops introduces small magnitude oscillations of the HSC, MPPs, and lymphoid cells, the myeloid cell dynamics are monotonic and the system robustly returns to its steady state over times that are consistent with those established in previous experiments (Reynaud et al., 2011) for similar perturbation studies.

Following previous modeling studies, we modeled CML by introducing a parallel lineage of mutant leukemic cells (denoted by the superscript L) but with the behavior of that lineage coupled at many points to the behavior of non-mutant cells, and vice versa. In particular, the model for normal and CML cells shares the same lineage structure and feedback architecture with both normal and mutant cell types providing a source of regulating factors, and although all the leukemic parameters (Appendix 1—table 3) could be different from their normal counterparts (Appendix 1—table 2), we begin by assuming the only difference between the two lineages is a decrease in the feedback strength for leukemic HSC (HSC^L; S^L), as indicated by p₀^L in the schematic in Figure 4A. This makes the leukemic cells less responsive to negative feedback and enables leukemic cells to gain a competitive advantage for growth. One candidate mediator of this negative feedback is CCL3, previously shown to inhibit self-renewal and division of normal HSC but HSC^L are less sensitive to its inhibitory regulation (Eaves et al., 1993, Dürig et al., 1999, Baba et al., 2016; Staversky et al., 2018). An example of CML hematopoiesis is shown in Figure 4B, where it is seen that, after the introduction of a few HSC^L at equilibrium of the normal hematopoietic system, the CML cells (dashed curves) repopulate the BM at the expense of normal cells (solid curves). Because of negative feedback, the system will eventually reach a new equilibrium consisting solely of leukemic cells. See Appendix 1—table 5 for the leukemic parameter values, and Figure 4—figure supplements 1–5 for parameter sensitivity studies of systems containing both normal and CML cells.

Figure 4 with 6 supplements see all

Download asset Open asset

We then perturbed each of the leukemic parameters within 10% of the values in Appendix 1—table 5 and found that only the leukemic stem cell self-renewal parameters—the maximal HSC^L self-renewal fraction $p_{0,max}^L$ and the feedback gain ${\gamma }_1^L$ on the HSC^L self-renewal fraction—have the potential to significantly influence the results. The results are insensitive to changes in the other leukemic cell parameters (see Appendix 1, Section 9, Figure 6—figure supplements 2–4). These results are characteristic of even larger changes in the base parameters.

In this and subsequent parameter investigations, we constrained $p_{0,max}^L$ to be less than or equal to $p_{0,max}$ , the maximal self-renewal fraction of the normal HSC, motivated by the paucity of evidence that $p_{0,max}^L$ is larger than $p_{0,max}$ , coupled with experimental data suggesting that $p_{0,max}^L$ is less than or equal to $p_{0,max}$ . For example, CML long-term culture initiating cells (LTC-IC; thought to be phenotypically similar to stem cells) decrease significantly in in vitro cultures while the number of normal LTC-IC is unchanged, consistent with a relative decrease in self-renewal probability of the CML cells (Udomsakdi et al., 1992). In vivo, HSC self-renewal can be assessed directly through transplantation studies. In this regard, CML HSC engraft immunodeficient mice variably and inefficiently compared to normal human HSC (Wang et al., 1998) while HSC from BCR-ABL1 transgenic mice exhibit an engraftment defect upon secondary transplantation into syngeneic recipients (Schemionek et al., 2010). Both results are suggestive of a relative decrease in self-renewal capacity of BCR-ABL1+ stem cells.

Next, we performed a sweep through leukemic stem cell self-renewal parameters $p_{0,max}^L$ and ${\gamma }_1^L$ for each of the eligible parameter sets for normal hematopoiesis (see below). We found that for the terminally differentiated cell proportion to be at least 50% leukemic (darker regions), there are biological constraints upon the combination of $p_{0,max}^L$ and ${\gamma }_1^L$ (Figure 4C). As the heat map shows, in order for CML to dominate hematopoiesis (e.g., terminally differentiated cell proportion >50% leukemic), the CML stem cells should have $p_{0,max}^L$ sufficiently close to $p_{0,max}$ and ${\gamma }_1^L$ should be sufficiently small. As the ratio $p_{0,max}^L$ / $p_{0,max}$ decreases from 1, the system requires smaller feedback gains ${\gamma }_1^L$ to compensate and allow for CML to develop. Further, there are threshold values of the parameters required for CML hematopoiesis to prevail. Namely, the system is dominated by normal cells (CML cells do not ‘take over’) when $p_{0,max}^L$ / $p_{0,max}$ is sufficiently large or when ${\gamma }_1^L/{\gamma }_1$ is sufficiently small.

To further examine these biological constraints, we calculated characteristic effective self-renewal fractions for normal and leukemic stem cells, defined as ${\stackrel{-}{p}}_0^L=p_{0,max}^L/(1+{\gamma }_1^L\stackrel{-}{N})$ and ${\stackrel{-}{p}}_0=p_{0,max}/(1+{\gamma }_1\stackrel{-}{N})$, where $\stackrel{-}{N}={10}^5$ , a characteristic value for the size of the MPP population based on MPP steady state values (Manesso et al., 2013). The relative fitness of the CML cells defined by the ratio of characteristic values of the HSC^L and HSC self-renewal fractions: ${\stackrel{-}{p}}_0^L$ / ${\stackrel{-}{p}}_0$ . Here, all eligible parameter sets representing the states of the normal system are considered and the leukemic parameters $p_{0,max}^L$ / $p_{0,max}$ and ${\gamma }_1^L/{\gamma }_1$ are varied from 0.6 to 1.0 and 0.1–0.6, respectively. In Figure 4D, we examined the relative fitness of leukemic cells through the distribution of the ratio of characteristic values colored by leukemic cells outcompeting normal (orange) and normal cells maintaining majority (blue). As expected, the larger the relative fitness, the more likely that CML will take over the system and dominate hematopoiesis. For further analysis of the leukemic parameter combinations for CML hematopoiesis and under treatment, see Figure 4—figure supplement 6, Figure 6—figure supplements 2–4, Figure 7—figure supplement 2, Figure 8—figure supplement 2, Appendix 1, Sections 9–11, and Appendix 1—figures 11–17 for details.

To test whether our mathematical model recapitulates known features of CML biology, we simulated a published transplant experiment in a transgenic mouse model of CML that recapitulates the main features of human CML (Reynaud et al., 2011). In this experiment, either HSC^L or leukemic MPPs (MPP^L) were implanted into sublethally irradiated mice (Figure 5A). Transplantation of HSC^L enables engraftment and myeloid cell production that leads to CML. On the other hand, transplanting MPP^Ls does not allow for long-term engraftment but results in a larger fraction of donor-derived lymphoid cells after 35 d (Figure 5B). This study presented evidence that IL-6 produced by differentiated myeloid cells reprograms these MPP^L progenitors toward a myeloid fate (Reynaud et al., 2011). As described in ‘Methods’ and Appendix 1, Section 4, we modeled this experiment by reducing the number of cells in equilibrium to mimic the effects of sublethal radiation. We explored a range of possible reductions of HSC^L, MPP^L, and differentiated myeloid and lymphoid cells and tracked the outcomes when 4000 HSC^L or MPP^L were introduced after the decrements from equilibrium. We then discarded those parameter sets that did not yield results consistent with a simple majority of myeloid cells for HSC^L transplant and a simple majority of lymphoid cells for MPP^L transplant (Reynaud et al., 2011). In particular, 85 parameter sets were discarded, leaving a total of 478 parameter sets remaining. Characteristic results are presented in Figure 5C and D when the reductions for HSC^L, MPP^L, and terminally differentiated cells were 55, 35, and 10%, respectively, from their equilibrium values. See Figure 5—figure supplement 1 for results using other decrements, the removed parameter set criteria (Figure 5—figure supplement 2), and Figure 5—figure supplement 3 for the final parameter distributions. When HSC^L are transplanted (solid curves), the donor-derived MPP^L (Figure 5C, left) rapidly increased as did the terminally differentiated myeloid and lymphoid cells (Figure 5C, right). Consistent with the experiments, there is a larger fraction of donor-derived myeloid cells than lymphoid cells after 30 d (Figure 5D). In contrast, when MPP^L are introduced (dashed curves), their population decreases (Figure 5C, right) because the MPP^L do not stably engraft. Concomitantly, there is burst of donor-derived myeloid and lymphoid cells at early times (Figure 5C, right) as the transplanted MPP^L differentiate.

Figure 5 with 3 supplements see all

Download asset Open asset

The early time dynamics of the myeloid and lymphoid cells depend on the specific values of the MPP^L self-renewal (p₁) and fate control (q₁) fractions, whose values in turn depend on the number of myeloid cells through negative feedback regulation. In particular, if ${\displaystyle 1-p_1>2q_1}$, then more myeloid than lymphoid cells will be produced at early times, as in Figure 5C (right), whereas more lymphoid cells will be produced if ${\displaystyle 1-p_1<2q_1}$ . In both cases, because the MPP^Ls do not stably engraft and instead differentiate into lymphoid and myeloid cells, we observe that there is a larger fraction of donor-derived lymphoid cells after 30 d (Figure 5D), consistent with the experiments. This occurs because there is a decreasing flux of differentiating cells since there is no stable engraftment and the lymphoid cells are longer-lived (smaller death rate) than the shorter-lived myeloid cells, which have a larger death rate.

We next explored the effects of TKI therapy on CML in the model. While the overall size of the phenotypic HSC compartment is not increased in CML patients (Jamieson et al., 2004), the proportion of HSC^L in the BM can vary widely across newly diagnosed CML patients from a few percent to nearly 100% (Petzer et al., 1996; Diaz-Blanco et al., 2007; Abe et al., 2008; Thielen et al., 2016). We therefore investigated how the HSC^L load in the BM affects therapy outcomes. We used one eligible parameter set (see Appendix 1—table 4), out of all 478 parameter sets all of which are capable of characterizing the normal state of our simplified model of the hematopoietic system and one choice of leukemic parameters (see Appendix 1—table 4) in which the only difference between normal and leukemic cells is that the HSC^L are one half as sensitive to negative feedback regulation compared to the normal HSC (${\displaystyle \frac{{\gamma }_1^L}{{\gamma }_1}=0.5}$). The initial condition was obtained by simulating the development of CML, analogous to that shown in Figure 4B, prior to initiating therapy. TKI treatment was initiated at three different times to achieve varying leukemic stem cell load (6, 18, and 36 mo) and was simulated by introducing a death rate of HSC^L and MPP^L proportional to their proliferation rates, with the HSC^L proliferation rate lower than that of normal HSC (Jørgensen et al., 2006). While some studies have shown that primitive CML stem/progenitor cells are relatively resistant to killing by TKIs in vitro (Graham et al., 2002; Corbin et al., 2011), clinical studies suggest that long-term TKI therapy can decrement the CML stem cell compartment, at least in some patients (Etienne et al., 2017; Chen et al., 2019), consistent with mathematical modeling of patient BCR-ABL1 transcript data (Tang et al., 2011). This supports the concept that TKIs possess a measure of leukemia stem cell killing ability, and we therefore included this effect in our model. The TKI treatment parameters were the same across the three cases. See ‘Methods’ and Appendix 1 for details and Appendix 1—tables 4 and 5 for parameter values. Thus, these cases can be thought of as representing the response of one virtual patient to TKI therapy implemented at different times after disease initiation.

At an early treatment time with lower (<90%) initial HSC^L fractions (HSC^L, Figure 6A), the numbers of MPP^L (blue-dashed), leukemic terminally differentiated lymphoid (light-green-dashed), and myeloid (dark-green-dashed) decrease rapidly at the early stages of treatment and are accompanied by a rapid increase in HSC^L due to the loss of negative feedback from the MPP^L. This loss of negative feedback from the MPP^L also results in a rapid increase in the number of normal HSC (black solid curves) that subsequently drives an increase in the normal MPPs (blue solid). The increased number of HSC and HSC^L decreases their division rates due to the autocrine negative regulatory loop as well as the division rates of the MPPs and MPP^L through feedforward negative regulation. This decreases the flux into the terminally differentiated cell compartments (both normal and leukemic), thereby decreasing their numbers at early times. At later times, both the HSC^L and MPP^L gradually decrease in response to TKI-induced cell death, which drives an accompanying decrease in the leukemic differentiated cells. A small, transient increase in MPP^L is observed before the gradual decline. This is driven primarily by the increase in flux into the MPP^L compartment by differentiating HSC^L, although there is also a small contribution from the feedforward regulation of the MPP^L division rate, which lowers the effectiveness of TKI therapy on the MPP^L. Both of these effects are reduced as the HSC^L numbers are decreased by TKI therapy. This in turn increases the effectiveness of TKI therapy in killing leukemic cells at later times and enables the normal cells (solid curves) to rebound toward their pre-leukemic equilibrium values.

Figure 6 with 4 supplements see all

Download asset Open asset

At intermediate treatment time with larger (90–99%) HSC^L fractions (Figure 6B), the responses of the leukemic and normal cells to TKI treatment at early times are qualitatively similar to those observed in the previous case although the effects are more pronounced. The increase in HSC^L is much larger than the previous case because there are fewer normal cells to compete with in the BM. This significantly decreases the HSC/HSC^L and MPP/MPP^L division rates through the negative feedback/feedforward regulation and correspondingly the rates of TKI-induced cell death. Accordingly, at later times the MPP^L population rebounds, driven by the flux of differentiating HSC^L, and eventually the system reaches a state in which both normal and leukemic cells coexist. The stem cell compartment is dominated by HSC^L which are largely quiescent, while the multipotent progenitor and terminally differentiated cell compartments have a higher fraction of normal cells. This is consistent with experimental results from mouse models (Reynaud et al., 2011) and our own unpublished data. In this scenario, BCR-ABL1 transcript levels in the peripheral blood are ~1–9%, but the patient would not respond further to TKI treatment and hence would not reach MR3; this has been observed clinically including one of the patients in our study (see below). The small flux of differentiating normal and leukemic stem and progenitor cells, combined with the negative feedback loops on the self-renewal and branching fractions, supports nearly steady populations of differentiated cells.

When given years to develop and a late time to treatment, the HSC^L fraction is nearly 100% (Figure 6C), and there are so few normal stem cells that the leukemic cells easily maintain nearly 100% of each cellular compartment even in the presence of TKI therapy. Aside from a short-lived, transient decrease in MPP^L (and differentiated leukemic cell) numbers, the leukemic cells are largely unresponsive to TKI therapy because the feedback/feedforward negative regulation of stem and progenitor cell division rates makes these rates so low that the TKIs are largely ineffective in killing the leukemic stem and progenitor cells. As in the previous case, the negative feedback regulation and the small fluxes of differentiating leukemic stem and progenitor cells enables the system to approach a steady state containing only leukemic cells.

In Figure 6D, we plot the simulated BCR-ABL1 transcript levels over time for the three scenarios. As described in ‘Methods,’ the transcripts are modeled using a relative ratio of leukemic and normal terminally differentiated myeloid and lymphoid cell numbers. The solid blue curve corresponds to CML using the treatment time from Figure 6A, which responds to TKI therapy. Just as in the clinical data (symbols), the response to TKI therapy produces a biphasic exponential decrease in BCR-ABL1 transcripts, which decreases below 10^–1, representing a so-called major molecular response (MMR or MR3), which represents a major goal of therapy in CML as the risk of relapse and leukemia-related death is virtually nonexistent once this milestone is achieved (Hochhaus et al., 2017). Consistent with previous interpretations, the rapid initial decrease in BCR-ABL1 transcripts is due to TKI-induced cell death of MPP^L and the increase in normal HSC and MPPs, which induce corresponding changes in the myeloid and lymphoid cells (Figure 6A). The long-term, slower depletion of leukemic cells and the stable normal cell populations result in the second phase of the biphasic response. The simulated results compare well with clinical data from the DAISISON study of imatinib versus dasatinib in patients with newly diagnosed CML (Cortes et al., 2016) where the data corresponds to mean BCR-ABL1 transcripts, with standard deviations, adapted from Glauche et al., 2018 for patients who received imatinib (blue) or dasatinib (red).

The two other curves in Figure 6D correspond to the treatment times from Figure 6B (solid orange) and C (dashed orange). In these cases, the BCR-ABL1 transcripts do not decrease below the MR3 threshold, indicating that neither of these virtual patients responds adequately to TKI therapy. There is a partial response in patient from Figure 6B as the transcripts initially decrease due to TKI-mediated death of MPP^L, but this effect soon saturates because the leukemic stem cells are able to drive the regrowth and persistence of leukemic progenitor and differentiated cells. For the virtual patient with parameters from Figure 6C, there is essentially no change in the BCR-ABL1 transcripts when therapy is applied. These behaviors are consistent with those observed in CML patients with primary resistance to TKI therapy (Zhang et al., 2009; Yeung et al., 2012; Pietarinen et al., 2016).

The fundamental difference between these three virtual patients is the number of leukemic stem cells at the start of therapy, which occurs because treatment is initiated at different times following the development of CML (early—6 mo after CML initiation ~93% initial HSC^L fraction, intermediate—18 mo after CML initiation ~99% initial HSC^L fraction, late—36 mo after CML initiation ~99.99% initial HSC^L fraction). Our results suggest that the higher the HSC^L fraction at the start of therapy, the less effective the therapy. This follows from the feedback/feedforward regulation where increased numbers of HSC and HSC^L decrease their own proliferation rates as well as those of the MPPs and MPP^L (see Figure 6—figure supplement 1 and Appendix 1—figure 6 for further explorations of feedback/feedforward regulation of parameters). This reduces the effectiveness of TKI therapy as evidence suggests TKIs preferentially target dividing leukemic cells (Graham et al., 2002; Corbin et al., 2011) and suggests a mechanism why some patients are destined to do poorly with TKI therapy.

To test this hypothesis, we created BM chimeric mice containing both normal and leukemic (BCR-ABL1⁺) HSC by transplantation of BM from conditional BCR-ABL1 transgenic mice (Koschmieder et al., 2005) into unirradiated congenic recipient mice. Following stable engraftment, BCR-ABL1 expression is induced in transgenic HSCs by withdrawal of doxycycline (see ‘Methods’). These chimeras represent a novel and physiologically accurate in vivo model of early CML development that reflects interactions between normal and CML cells in a BM microenvironment unperturbed by radiation (Rodriguez et al., 2022). Two cohorts of chimeric mice bearing either a high HSC^L burden (94 ± 1.5% of the HSC population) or an intermediate HSC^L burden (58 ± 12%) were treated with dasatinib (25 mg/kg daily by oral gavage). Both cohorts showed a hematological response to TKI therapy, with decreased peripheral blood leukocyte counts (Figure 6E) and a decreased percentage of circulating granulocytes (Figure 6F). By contrast and consistent with the predictions of the quantitative model, while mice bearing smaller populations of HSC^L showed a decrease in the percentage of circulating BCR-ABL1⁺ granulocytes in response to TKI therapy, mice with the highest HSC^L burden showed virtually no decrease in circulating leukemic cells (Figure 6H). Because the level of circulating granulocytes reflects the proportion of BM HSC (Wright et al., 2001; data not shown), these results demonstrate that TKI therapy is unable to decrement the HSC^L compartment in mice with predominantly BCR-ABL1⁺ HSC at the start of treatment.

While analyses of clinical data also show that patients with lower leukemic stem cell burden are more likely to respond to TKI treatment (Thielen et al., 2016), some patients with a high percentages of leukemic stem cells at the start of treatment are nonetheless still capable of responding to TKI therapy (e.g., see Figure 3 in Thielen et al., 2016). This suggests that leukemic stem cell burden alone does not predict the molecular response to TKIs. To investigate this further, we tested the response to TKI treatment for each of our 478 parameter sets. In Figure 7A, we present the results using only one choice of leukemic parameters (see Appendix 1—table 5). Other choices of leukemic parameters give similar results (see Appendix 1—figure 13). The model outcomes bear a striking resemblance to the clinical data of Thielen et al., 2016. The leukemic stem cell fraction does influence TKI response, but treatment outcomes are seen to vary among virtual patients within the same initial leukemic stem cell load. We then asked what characteristics (e.g., parameter sets) distinguish whether a virtual patient achieves a MR3 response within 50 mo. We also varied the HSC^L parameters, taking into account several studies suggesting that CML stem cells are at least 5–10 times less sensitive to CCL3-mediated inhibition of self-renewal (Eaves et al., 1993, Chasty et al., 1995; Wark et al., 1998; Dürig et al., 1999); for example, $\frac{{\gamma }_1^L}{{\gamma }_1}$ should be less than 0.2. Further, since 10–15% of patients do not respond to TKI treatment even in the absence of BCR-ABL1 mutations (Hanfstein et al., 2012; Marin et al., 2012), we estimate $p_{0,max}^L$ / $p_{0,max}\approx 0.8$ from Appendix 1—figure 11 to roughly match this proportion of nonresponding patients. Taken together, this suggests that the effective leukemic stem cell fitness would be ${\displaystyle \frac{{\overline{p_0}}^L}{p_{0,max}}\approx 0.7}$. We thus varied the HSC^L parameters accordingly. In Figure 7B, we plot the results for ${\displaystyle \frac{{\overline{p_0}}^L}{p_{0,max}}\ge 0.7}$ as a bivariate histogram for ${\displaystyle \frac{{\overline{p_0}}^L}{p_{0,max}}}$ and $p_{0,max}$ with proportion of response (blue) and nonresponse (orange) for every parameter set. Surprisingly, we found that although the fitness ${\displaystyle \frac{{\overline{p_0}}^L}{p_{0,max}}}$ impacts response, the parameter that clearly distinguished responders from nonresponders was the maximal self-renewal fraction p_0,max for normal stem cells, shown in Figure 7B (marginal y-axis). See Figure 7—figure supplement 1 for the distributions as a function of the other parameters using the single leukemic parameter set from Figure 7A (see Appendix 1—table 5), and Appendix 1—figures 15–16 in Appendix 1, Section 10 for different bivariate distributions corresponding to different choices of the minimum fitness ${\displaystyle \frac{{\overline{p_0}}^L}{p_{0,max}}}$ . In particular, larger values of p_0,max and γ₁ are correlated with a decreased response to TKI therapy after leukemia develops. Although these parameters are associated with normal HSC, the self-renewal fraction p₀^L of HSC^L and feedback strength ${\gamma }_1^L$ depends on these parameters since we assumed the fitness of the CML stem cells ${\displaystyle \frac{{\overline{p_0}}^L}{p_{0,max}}}$ is larger than a minimum threshold. Therefore, increasing the self-renewal fraction of the normal stem cells has the effect of also increasing the fitness of the CML stem cells.

Figure 7 with 2 supplements see all

Download asset Open asset

To understand further the differences between response and nonresponse to TKI therapy, we took the parameter set from Figure 6A as a representative patient for response and selected an arbitrary nonresponsive parameter set (Appendix 1—table 6) to be a representative patient for nonresponse. In Figure 7C, we show that the effective p₀^L (the fraction of HSC^L self-renewal after feedback) for nonresponders (orange) is larger after TKI therapy is applied than for responders (blue). In particular, as TKI treatment kills the leukemic progenitors, this increases the effective self-renewal fraction for both normal and leukemic stem cells because of the release of negative feedback. When the maximum self-renewal p_0,max is larger, the leukemic stem cells experience an acute increase in self-renewal, resulting in their dominance over normal stem cells that then leads to a decreased response to TKIs.

Clinical data provide support for this mechanism of resistance. Patients with clonal hematopoiesis, in which there is a dominant clone driving hematopoiesis, exhibit predominantly normal hematopoiesis but frequently have mutations in the genes, such as TET2, DNMT3A, and ASXL1, that are known to increase stem cell self-renewal (Steensma, 2018). Clinical data shows that CML patients whose blood cells have mutations in TET2 and ASXL1, some of which may exist prior to development of CML, frequently exhibit a poor response to TKI therapy (Kim et al., 2017; Marum et al., 2017). Taken together, these data suggest that patients with higher stem cell self-renewal fare worse when their CML is treated using TKIs than patients with lower stem cell self-renewal.

Several measures of the response of CML patients to TKI therapy have been developed, based on BCR-ABL1 transcript levels in peripheral blood. Here we test a new, model-driven criterion for predicting patient response and compare the results with several other criteria currently used in the clinic. A major focus has been on the predictive value of the decline in transcripts over the first 3 mo of treatment, principally the so-called ‘early molecular response’ or EMR (defined as BCR-ABL1 transcripts <10% at 3 mo and <1% at 6 mo), where patients with >10% transcripts had significantly lower probability of achieving cytogenetic remission and decreased overall survival (Hanfstein et al., 2012; Marin et al., 2012). Subsequently, there was an effort to improve the predictive power by focusing on the velocity of reduction in transcripts (Branford et al., 2014; Hanfstein et al., 2014; Pennisi et al., 2019). Because the best predictor of patient response to TKIs, the self-renewal fraction of normal stem cells, is very difficult to measure clinically, we searched for an alternative criterion that could accurately predict patient response and could still be measured using the data collected in standard practice. Therefore, we focused on alternative time frames and calculation methods for assessing BCR-ABL1 transcript levels (Figure 7D). It is important to be able to predict the long-term TKI response early after starting treatment in order to enable changes in therapy. However, since both responders (blue) and nonresponders (orange) may show significant decreases in the transcript levels in the first months of treatment, it was difficult to distinguish between the two at relatively early time points. By contrast, responses in the 3–6-month time frame make it easier to identify the different behaviors of responders and nonresponders (Figure 7D).

By calculating the relative changes of the transcript levels from 3 to 6 mo, we developed a prognostic formula: $PF(3,6)=\left(BCRABL1\left(6\right)-BCRABL1\left(3\right)\right)/BCRABL1\left(6\right)$. We found that optimizing for sensitivity (TPR, the true positive rate) and specificity (1-FPR, with FPR being the false positive rate) resulted in a prognostic threshold of $PF\approx -3.2,$ with sensitivity of ${\displaystyle \approx 0.91}$ and specificity of ${\displaystyle \approx 0.91}$ (Figure 7E, orange curves) compared to the optimal sensitivity and specificity of the velocity-based prognostic (${\displaystyle \approx 0.89}$ and ${\displaystyle \approx 0.88,}$ respectively) and ${\displaystyle \approx 0.77}$ and ${\displaystyle \approx 0.98}$ for EMR 1% with our parameter sets. This demonstrates that this prognostic tool had higher sensitivity and specificity than previously developed predictive criteria in separating responders (${\displaystyle <-3.2}$) from nonresponders (${\displaystyle >-3.2}$), where response is defined as achievement of MR3 within a clinically relevant timeframe of 18 mo. We also tested the various prognostics at the 0–3-month interval as is the current clinical practice, but that resulted in lower predictive power (Figure 7E, blue curves). These results highlight the importance of including the 3–6-month TKI response in predicting the long-term outcome of treatment, instead of considering only the first 3 mo. See Appendix 1—figures 7 and 17 in Sections 8 and 12, respectively, for further discussion, comparison of additional prognostic criteria, and the effect of leukemic parameters.

We then applied our prognostic criterion to anonymized CML patient data (see ‘Methods’) to determine clinical significance and utility. The prognostic tests shown in Figure 7E were calculated for both the first 3 mo and the subsequent 3–6-month period after the start of therapy, for the patients who were treated with the same TKI and dosage for the full 6 mo. All the prognostic tests achieved a more accurate prediction of patient outcome using the 3–6-month data compared to the same test applied to the first 3 mo (Figure 7F). To expand clinical utility, the prognostics were calculated for cases where TKI therapy was changed (due either to toxicity or an inadequate response) but then maintained for a subsequent 6-month period, which were added to the data from Figure 7F. The aggregated data (Figure 7G) reaffirms the improved accuracy in prediction using the 3–6-month transcript data compared to that from 0 to 3 mo. Over the first 3 mo, all the prognostic criteria performed similarly. Although the number of patients was small, the results suggest that our prognostic criteria may perform better than the EMR and velocity-based prognostics that are in current clinical use. For comparisons between the prognostic criteria and time frames with patient data, see Appendix 1, Section 7, Appendix 1—figures 8–10.

Our model suggests that combination therapy to modulate the stem cell self-renewal rate, in addition to directly targeting the leukemic HSC and MPPs with TKI therapy, might counteract TKI treatment resistance mediated by high stem cell self-renewal. Such pro-differentiation therapy could be accomplished through either direct stimulation of differentiation or through suppression of self-renewal. In our modeling experiments, we explored the impact of this approach through the suppression of self-renewal (see ‘Methods’ for details). To begin the exploration of the combined TKI-differentiation therapy, we performed this combination therapy on each of our 478 parameter sets, which represents a population of CML patients with person-to-person variability. We then recorded which parameter sets achieved MR3 within 50 mo for each strength of the differentiation therapy (Δ), where Δ is a dimensionless constant greater than 0 that constantly suppresses stem cell self-renewal (both normal and CML) in the setting of combination therapy (see ‘Methods’ and Appendix 1, Equations 40 and 41). Using these data, Figure 8A depicts the proportion of parameter sets achieving MR3 given a strength of differentiation therapy of Δ. As differentiation therapy strength increases from zero, the proportion of parameter sets that achieve MR3 increases before leveling off between Δ = 0.2–0.3, with maximum efficacy occurring at a strength of differentiation Δ of about 0.24. The efficacy of combination therapy then begins to decline rapidly, and with too great a strength of differentiation treatment, the combination therapy becomes inferior to TKI therapy alone.

Figure 8 with 2 supplements see all

Download asset Open asset

To investigate how combination therapy effectively targets resistance, and the mechanism of the decreased efficacy of combination therapy in achieving MR3 when Δ is large, we returned to examining parameter distributions. Figure 8B depicts the same distribution of p_0,max as in Figure 7B, but overlaid with a second histogram (hatched regions) to denote the effect of the differentiation therapy at the point of maximum efficacy (see Figure 8—figure supplement 1 for all the parameter distributions). The two types of hatching reveal important factors that determine under which conditions combination therapy improves or impairs response. The orange hatching represents transition from response to nonresponse by combination therapy; this occurs in individuals with the lowest p_0,max. In these cases where stem cell self-renewal is already close to the ideal effective self-renewal fraction of 0.5, differentiation therapy pushes too many normal cells into differentiation, causing the normal cell populations to deplete themselves and decreasing the efficacy as Δ increases beyond 0.24. In contrast, the blue hatching shows the desired scenario of nonresponding individuals with high p_0,max becoming responders and achieving MR3 within 50 mo due to the combination therapy.

To understand further the mechanisms underlying the efficacy of combination therapy, we explored treatment dynamics (changes in BCR-ABL1 transcript levels) and the rates of change of both normal and leukemic stem cell populations for the nonresponsive individual from Figure 7D. By applying two different strengths of differentiation therapy (Δ = 0.24 and 0.5) in combination with TKI therapy, for this individual both strengths are able to achieve MR3 at ~18 mo (Figure 8C) in contrast to TKI monotherapy, which resulted in a failure to reach MR3 (Figure 7B, orange). Figure 8D and E show how the rates of change in the size of the normal and leukemic stem cell compartments vary with respect to time for the three different Δ values. For TKI therapy alone (Δ = 0), rates of growth of both the normal and leukemic stem cell populations show an increase as a result of the loss of negative feedback due to TKI-induced killing of MPP^L, but the leukemic stem cells experience a much greater numerical increase and outcompete normal cells, resulting in a system that exhibits resistance to the TKI therapy. Under conditions of maximum efficacy (Δ = 0.24), the normal stem cell population rate of change still increases rapidly but to a maximum level below that for TKI monotherapy before decreasing more rapidly to zero as the system re-equilibrates. In contrast, under the same conditions the rate of change of the leukemic stem cell population is greatly reduced and becomes negative after normal stem cells begin to outcompete the leukemic stem cells. Under conditions of stronger differentiation therapy (Δ = 0.5), although the accumulation rate of the normal stem cells is substantially reduced, the growth rate of the HSC^L is immediately negative. This enables the normal HSC to easily outcompete the leukemic cells and restore the system to the normal state. Effectively, for large values of HSC self-renewal and corresponding feedback gains, the differentiation promoter acts to bring the self-renewal fraction of the normal HSC closer to that of the HSC^L, which then enables the TKI therapy to disadvantage the leukemic cells, and allow for repopulation and dominance by normal cells. For the effect of combination therapy on different combinations of leukemic parameters, see Figure 8—figure supplement 2 and Appendix 1, Section 11 for further analysis.

The classical depiction of hematopoiesis is a hierarchy of cell types starting with the hematopoietic stem cell at the top, followed by progenitors and ultimately ending with mature cells located in the peripheral blood. Therefore, we model hematopoiesis using a lineage ODE model (Roeder et al., 2006; Komarova and Wodarz, 2007; Horn et al., 2008; Foo et al., 2009; Lander et al., 2009; Marciniak-Czochra et al., 2009; Manesso et al., 2013; Buzi et al., 2015; Hähnel et al., 2020; Pedersen et al., 2021) to describe cellular growth dynamics. The modeling allows us to follow the similar hierarchical structure by creating an order of differentiation. Our branched lineage model of hematopoiesis is simplified and only models HSCs, progenitor cells (MPPs), and two types of terminally differentiated cells (myeloid and lymphoid cells). The model can easily include additional cell types, such as committed progenitor cell types, which will provide an additional level of detail. The model consists of two dividing cell types consisting of S (HSC) and P (MPP) cells with a division rate associated to the cells (η₁ and η_2, respectively). The S cells have the ability to self-renew with fractions (p₀) or differentiate (1 p₀). The P cells have the ability to self-renew with fraction (p₁) or differentiate into either TD_l (lymphoid) or TD_m (myeloid) cells (q₁ or 1-p₁-q_1, respectively). Both S and P cells do not die within the model. The terminal cells form the majority of the hematopoietic system and consist of TD_l and TD_m cells. TD_m and TD_l cells are postmitotic and die at rates d_m and d_l, respectively. The following equations (Equations 1–4) describe the dynamics of the system:

Further expanded forms of the equations are shown in Appendix 1, Equations 26–29, with the addition of feedback regulation for each of the parameters.

We use an automated method developed by Savageau and collaborators (Savageau et al., 2009; Fasani and Savageau, 2010; Lomnitz and Savageau, 2013; Lomnitz and Savageau, 2016) that separates models by distinct qualitative behaviors at steady state. The strategy is to deconstruct the model of interest at steady state to focus on cases where one production term and one loss term dominate, which gives a dominant subsystem (S-System). This implies that particular inequalities hold in order to ensure the production and loss terms chosen are larger than the others. The inequalities are evaluated at the S-system’s steady state to assess self-consistency. If the inequalities are satisfied, the system is self-consistent and the regions where equality holds form boundaries that pertain to a particular qualitative behavior associated with the system. The interior region (where strict inequality holds) is termed a domain in design space. If all the S-systems associated with a model do not have any equilibria that are self-consistent or equilibria that are stable, then the model is rejected. The benefits of this method are that it does not require prior knowledge about parameter values, and it can enumerate the different types of qualitative dynamics a certain system may have. By eliminating subsets of parameters for which the equilibrium is unstable, this approach will automatically select models that are robust to parameter variation due to stability. When we applied this method to the ODE system in Figure 2F (Appendix 1, Equations 26–29), we found that only the four model classes shown in Figure 1B were accepted. See Appendix 1, Section 1, for further details.

To approximate biologically relevant parameters for the model a grid-search algorithm was employed. Parameters were sampled using a random uniform distribution for each parameter (see Appendix 1, Section 3.1). Once parameter values were chosen, the model was simulated to steady state. If a parameter set resulted in steady-state values consistent with the order of magnitude in Manesso et al., 2013, the parameter set was accepted, otherwise it was rejected. Specifically, these inequalities had to be satisfied 10⁴ < HSC < MPPs with MPPs fixed at 10⁵ and MPPs < TD_l < T_Dm. For 10⁶ iterations, a sample of 1493 parameter sets were accepted. The distribution for these parameter sets is shown in Appendix 1—figure 4. To further explore the effect of the feedforward interaction, these parameter sets were reduced to the 563 sets with γ₅ > 0.01. The distribution for these sets is shown in Appendix 1—figure 5. The parameter sets used in Figures 3—7 are provided in Appendix 1—table 4.

To model CML development in the presence of normal hematopoietic cells, we introduce a new leukemic cell type for each compartment. Each compartment is then composed of both normal and leukemic subcompartments, which exhibit feedback together as a single compartment. We assume the only difference between the two cell lineages is the feedback strength for leukemic HSC self-renewal. This small difference gives the leukemic lineage a competitive advantage for growth, consistent with the ability for leukemic HSC to initiate CML (Reynaud et al., 2011; Holyoake and Vetrie, 2017) and the differential response of the normal and leukemic cells to CCL3 (Eaves et al., 1993, Baba et al., 2013), which negatively regulates stem cell self-renewal. The full equations used in the model are shown in Appendix 1, Equations 30–37.

The model was tested by simulating the transplant experiments (Figure 4C and D) of Reynaud et al., 2011 where HSC^L or MPP^L were implanted into sublethally irradiated mice and terminal cell counts were measured after 35 d. We used two parallel lineages of leukemic cells with identical parameters to mirror the two leukemic cell populations of the experiment. To mimic the effects of sublethal radiation, we reduced the cell populations from their equilibrium values by variable amounts. The HSC^L depletions varied between 50 and 70% and the MPP^L depletions varied between 30 and 50% while both terminal cells were depleted by 10%. After depletion, an additional 4000 cells of either stem or progenitor types were transplanted in accordance with the experiment. We then discarded the 85 parameter sets that were not consistent with the results from Reynaud et al., 2011, leaving 478 eligible parameter sets. The results shown in Figure 4 used depletions of 55% for HSC^L, 35% for MPP^Ls. See Appendix 1, Section 4 and Figure 5—figure supplements 1–3 for results using other decrements, the discarded parameter sets, and the final parameter distributions.

To account for the treatment by TKIs, additional proliferation-dependent death terms are added to the equations for leukemic stem cells and leukemic progenitor cells shown in Appendix 1, Equations 38 and 39 (parameter values are given in Appendix 1—figure 5). These represent the ability of TKIs to induce cell death in the leukemic cells. Both cell types have unique death rates, to reflect TKIs having different efficacy in killing stem cells and progenitors. The death rates were selected using a single parameter set to ensure a reasonable biphasic curve for BCR-ABL1 transcript levels compared to patient transcript levels from Glauche et al., 2018. The same death rates were then used across every parameter set to ensure consistency. In addition to these changes upon initiation of TKI therapy, the leukemic stem cell division rate is reduced. This reduction models the ability of TKIs to drive leukemic stem cells to quiescence (Jørgensen et al., 2006).

To approximate the BCR-ABL1 transcript levels, we used a method based upon (Michor et al., 2005). We use the cell counts of both normal and leukemic terminal cells for both myeloid and lymphoid lineages. The terminal cells are used as in our model they are the closest to peripheral blood in which transcript levels are measured clinically. This results in the following measure for BCR-ABL1 transcript levels: $\frac{BCR-ABL1}{BCR}=\frac{T{D}_L^L+T{D}_M^L}{T{D}_L^L+T{D}_M^L+2(T{D}_L+T{D}_M)}$ .

Combination therapy consists of simultaneously employing TKI therapy, described in ‘Methods’ and Appendix 1, Sections 2–3, and the addition of a new differentiation therapy. To model differentiation therapy, we altered the form of p₀ by including a new constant repressive force that affects both normal and leukemic self-renewal, resulting in $p_{0,new}=\frac{p_{0,max}}{1+{\gamma }_{1}P+\mathrm{\Delta }}$ and $p_{0,new}^L=\frac{p_{0,max}^L}{1+{\gamma }_1^{L}P+\mathrm{\Delta }}$, where $P=x_{P^L}+x_P$ and $\mathrm{\Delta }$ is the differentiation therapy strength. We then performed combination therapy using our existing parameters, swept through differentiation therapy strengths and recorded which parameter sets achieved MR3 within 50 mo. We then determined that a differentiation therapy strength of D = 0.24 resulted in the highest proportion of parameter sets that achieved MR3 response. The full equations for combined TKI and differentiation therapy are shown in Appendix 1, Equations 40 and 41.

Data from newly diagnosed CML patients (n = 21) treated with TKI therapy at UCI Health were obtained under an honest broker mechanism from the electronic health record under Exemption 4 for human subjects research.

C57BL/6J female mice (Jackson Laboratories), 6–12 wk of age were used for irradiation and myeloid depletion experiments. Conditional BCR-ABL1 double transgenic mice (Koschmieder et al., 2005) were obtained from Dr. Emmanuel Passegue (Columbia University). All protocols in mice were approved by the Institutional Animal Use and Care Committee of University of California, Irvine.

To achieve selective depletion of HSCs, a 50 cGy dose of irradiation from X-ray source (Precision X-rad 320) was applied. Control mice did not receive irradiation. The distribution of time points at which observations were made (days 1, 3, and 7 post-irradiation), and the number of mouse replicates to use at each time point (between 2 and 7, totaling 13 mice), were informed by our Bayesian hierarchical framework for optimal experimental design (Lomeli et al., 2021).

RB6-8C5, an anti-Gr1 antibody (catalog # BE0075, BioXCell) or isotype control (catalog # BE0090, BioXCell) was injected intravenously, 50 µg per mouse, and mice sacrificed 24 hr later.

In irradiation experiments, mice were pulsed with BrdU by IP injection of 200 µl of 10 mg/ml BrdU in DPBS. BrdU flow kit (552598) from BD Biosciences was used for detection of BrdU labeling in hematopoietic cells by flow cytometry.

BM cells from femur and tibia of control and dosed mice were isolated by flushing bones. Following lysis of red blood cells (RBC lysis buffer, eBiosciences), leukocytes were stained with CD34 antibody for 1 hr and subsequently incubated with a cocktail of biotinylated antibodies directed against lineage markers (CD3, Gr-1, B220, Ter119) and stem/progenitor markers (c-Kit, Sca-1, CD48) for 30 min. Streptavidin (SA)-conjugated fluorochrome was utilized to detect biotinylated antibodies. Following fixation, permeabilization, and DNase digestion, anti-BrdU antibody was used to assess BrdU incorporation. Events were acquired on FACS Arial II and analyzed with Flowjo v.10 software.

Monoclonal antibodies for flow cytometry were biotinylated mouse lineage panel (559971, BD Biosciences), PE-CF594 Streptavidin (562318, BD Biosciences), anti-mouse CD48 (561242, BD Biosciences), anti-mouse CD34 eFluor450 (48-0341-82, eBiosciences), anti-mouse Sca-1-PE (108108, BioLegend), anti-mouse c-Kit-APC (17-1171-82, eBiosciences), and FITC BrdU flow kit (559619, BD Biosciences).

The full details of the CML mouse model will be published elsewhere (Jena et al., in preparation). Briefly, BM cells from conditional BCR-ABL1 double transgenic mice (CD45.2⁺) (Koschmieder, Gottgens et al. 2005) (40 million cells) were transplanted intravenously into unirradiated C57BL/6J recipients (CD45.1⁺CD45.2⁺) maintained on doxycycline to suppress BCR-ABL1 expression. Two months post-transplant, doxycycline was removed to allow induction of CML-like leukemia. Chimerism was assessed by percentage of CD45.1^– CD45.2⁺ granulocytes in peripheral blood. To generate chimeric mice with high (>90%) leukemic stem cell burden, the donor and recipient pair was reversed, with double transgenic mice transplanted with normal B6 BM. In mice with established CML-like leukemia (peripheral blood leukocytes > 20,000/μl and >40% circulating granulocytes), TKI treatment was initiated with dasatinib (25 mg/kg daily by oral gavage).

Code used to generate the figures, determine the S-Systems and files containing parameter sets are found at the following GitHub: https://github.com/jonatdr/CML_Treatment (copy archived at Jonatdr, 2023).

The complete ODE model of normal hematopoiesis is composed of the following equations:

The ODE model for CML hematopoiesis tracks the dynamics of both the normal and CML cells (superscript $L$), and assumes that both cell types provide and respond to feedback signaling, although the CML stem cells are slightly less responsive to negative feedback regulation, which gives them a fitness advantage. The complete system is given by

When TKI therapy is applied, Equations 34 and 35 are replaced with

where $TK{I}_{HSC}$ and $TK{I}_{MPP}$ denote TKI-induced death rates of the CML stem and MPP cells. With the introduction of differentiation therapy, Equations 30 and 38 instead become

The parameters for the model of normal hematopoiesis, and their descriptions, are listed in Appendix 1—table 2. The additional parameters needed to model the CML cell population, and the application of TKI therapy are given in Appendix 1—table 3.

3.1 Parameter distributions

A grid-search algorithm was used to find parameter values that demonstrate steady-state cell counts consistent with Manesso et al., 2013 and time of recovery to steady state from cellular perturbations consistent with Reynaud et al., 2011. The resulting distributions of the 1493 parameters are displayed along the diagonal of Appendix 1—figure 4. In the distributions, the feedback gain of the feedforward regulation is skewed toward smaller values across all parameter sets. To investigate the observed feedforward loop, we selected only the parameter sets with sufficiently large feedforward gain ${\displaystyle {\gamma }_5>0.01}$. This resulted in a reduced distribution of 563 parameter sets and is shown in Appendix 1—figure 5. The parameter ranges for the gridsearch were found using the method shown in the following pseudo-Python.

for i in range(1e6):

param_gridsearch_dist = dict(p0max = np.random.uniform(0.5, 1.0),
p1max = np.random.uniform(0.0,0.5), q1max = np.abs(np.random.uniform(0.0,0.49)),
eta1max = np.random.uniform(0,0.5), eta2max = 10**np.random.uniform(–2,1.5),gam2=10**np.random.uniform(–6,0), gam3=10**np.random.uniform(–6,0),
gam4=10**np.random.uniform(–6,0), gam5=10**np.random.uniform(–6,0),dL = 10**np.random.uniform(–4,1), dM = 10**np.random.uniform(–4,1))
y=hematopoiesis(param_gridsearch_dist)

if y[0,–1]>.01 and y[0,–1]<y[1,-1]<y[2,-1]<y[3,-1]:master_params.append(param_gridsearch_dist)

The specific parameters used to model the normal hematopoietic system in Figures 3—7 in the main text are given in Appendix 1—table 4. When CML cells are introduced, the additional parameters of a representative responder associated with the leukemic cell ODEs are given in Appendix 1—table 5. A representative parameter set for a nonresponder is shown in Appendix 1—table 6.

Appendix 1—figure 4

Download asset Open asset

Appendix 1—figure 5

Download asset Open asset

Appendix 1—figure 6

Download asset Open asset

Appendix 1—table 4

Parameter	Value
$p_{0,max}$	0.756641
$p_{1,max}$	0.357913
$q_{1,max}$	0.032241
${\eta }_{1,max}$	0.197639
${\eta }_{2,max}$	0.47088
${\gamma }_1$	0.513281
${\gamma }_2$	0.197639
${\gamma }_3$	0.47088
${\gamma }_4$	0.513281
${\gamma }_5$	0.543987
$d_L$	0.000165
$d_m$	0.000428

Appendix 1—table 5

Parameter	Value
$p_{0,max}^L$	0.756641
$p_{1,max}^L$	0.357913
$q_{1,max}^L$	0.032241
${\eta }_{1,max}^L$	0.197639
${\eta }_{2,max}^L$	0.47088
${\gamma }_1^L$	0.2566405
${\gamma }_2^L$	0.197639
${\gamma }_3^L$	0.47088
${\gamma }_4^L$	0.2566405
${\gamma }_5^L$	0.543987
$d_L^L$	0.000165
$d_m^L$	0.000428
$TK{I}_{HSC}$	0.201311
$TK{I}_{MPP}$	0.024757

Appendix 1—table 6

Parameter	Value
$p_{0,max}$	0.838481
$p_{1,max}$	0.009776
$q_{1,max}$	0.418627
${\eta }_{1,max}$	0.226025
${\eta }_{2,max}$	0.247591
${\gamma }_1$	0.676962
${\gamma }_2$	0.000219
${\gamma }_3$	0.000326
${\gamma }_4$	0.367564
${\gamma }_5$	0.168132
$d_L$	0.071273
$d_m$	0.155649
$p_{0,max}^L$	0.838481
$p_{1,max}^L$	0.009776
$q_{1,max}^L$	0.418627
${\eta }_{1,max}^L$	0.226025
${\eta }_{2,max}^L$	0.247591
${\gamma }_1^L$	0.338481
${\gamma }_2^L$	0.000219
${\gamma }_3^L$	0.000326
${\gamma }_4^L$	0.367564
${\gamma }_5^L$	0.168132
$d_L^L$	0.071273
$d_m^L$	0.155649
$TK{I}_{HSC}$	0.25
$TK{I}_{MPP}$	30