Jonathan Rodriguez,Abdon Iniguez,Nilamani Jena,Prasanthi Tata,Zhong-Ying Liu,Arthur D Lander,John Lowengrub,Richard A Van Etten
Predictive nonlinear modeling of malignant myelopoiesis and tyrosine kinase inhibitor therapy
cancer,nonlinear dynamics,hematopoiesis,treatment dynamics,mathematical model,chronic myeloid leukemia,Mouse
Chronic myeloid leukemia (CML) is a blood cancer characterized by dysregulated production of maturing myeloid cells driven by the product of the Philadelphia chromosome, the BCR-ABL1 tyrosine kinase. Tyrosine kinase inhibitors (TKIs) have proved effective in treating CML, but there is still a cohort of patients who do not respond to TKI therapy even in the absence of mutations in the BCR-ABL1 kinase domain that mediate drug resistance. To discover novel strategies to improve TKI therapy in CML, we developed a nonlinear mathematical model of CML hematopoiesis that incorporates feedback control and lineage branching. Cell–cell interactions were constrained using an automated model selection method together with previous observations and new in vivo data from a chimeric BCR-ABL1 transgenic mouse model of CML. The resulting quantitative model captures the dynamics of normal and CML cells at various stages of the disease and exhibits variable responses to TKI treatment, consistent with those of CML patients. The model predicts that an increase in the proportion of CML stem cells in the bone marrow would decrease the tendency of the disease to respond to TKI therapy, in concordance with clinical data and confirmed experimentally in mice. The model further suggests that, under our assumed similarities between normal and leukemic cells, a key predictor of refractory response to TKI treatment is an increased maximum probability of self-renewal of normal hematopoietic stem cells. We use these insights to develop a clinical prognostic criterion to predict the efficacy of TKI treatment and design strategies to improve treatment response. The model predicts that stimulating the differentiation of leukemic stem cells while applying TKI therapy can significantly improve treatment outcomes.
Predictive nonlinear modeling of malignant myelopoiesis and tyrosine kinase inhibitor therapy

Chronic myeloid leukemia (CML) is a blood cancer characterized by dysregulated production of maturing myeloid cells driven by the product of the Philadelphia chromosome, the BCR-ABL1 tyrosine kinase. Tyrosine kinase inhibitors (TKIs) have proved effective in treating CML, but there is still a cohort of patients who do not respond to TKI therapy even in the absence of mutations in the BCR-ABL1 kinase domain that mediate drug resistance. To discover novel strategies to improve TKI therapy in CML, we developed a nonlinear mathematical model of CML hematopoiesis that incorporates feedback control and lineage branching. Cell–cell interactions were constrained using an automated model selection method together with previous observations and new in vivo data from a chimeric BCR-ABL1 transgenic mouse model of CML. The resulting quantitative model captures the dynamics of normal and CML cells at various stages of the disease and exhibits variable responses to TKI treatment, consistent with those of CML patients. The model predicts that an increase in the proportion of CML stem cells in the bone marrow would decrease the tendency of the disease to respond to TKI therapy, in concordance with clinical data and confirmed experimentally in mice. The model further suggests that, under our assumed similarities between normal and leukemic cells, a key predictor of refractory response to TKI treatment is an increased maximum probability of self-renewal of normal hematopoietic stem cells. We use these insights to develop a clinical prognostic criterion to predict the efficacy of TKI treatment and design strategies to improve treatment response. The model predicts that stimulating the differentiation of leukemic stem cells while applying TKI therapy can significantly improve treatment outcomes.

Chronic myeloid leukemia (CML) is a myeloproliferative neoplasm of the hematopoietic system, which normally produces billions of mature myeloid and erythroid cells on a daily basis, is tightly regulated, and accommodates massive increases in the production of individual cell types in response to physiological and pathological stresses. The hematopoietic system is organized hierarchically as a collection of progressively more differentiated cells starting from a hematopoietic stem cell (HSC) located in the bone marrow (BM) and ending with postmitotic terminally differentiated myeloid and lymphoid cells (Rieger and Schroeder, 2012; Liggett and Sankaran, 2020). CML is characterized by an overproduction of myeloid cells including mature granulocytes (neutrophils, basophils, and eosinophils) and their immediate precursors (metamyelocytes, myelocytes, and promyelocytes), and of myeloid progenitors (Jamieson et al., 2004) including multipotential progenitors (MPPs) and committed progenitors (common myeloid progenitors [CMP], granulocyte-macrophage progenitors [GMPs], and megakaryocyte-erythroid progenitors [MEPs]). Untreated, the disease has three distinct phases (Chereda and Melo, 2015). In the initial ‘chronic’ phase, the differentiation of myeloid progenitors is essentially normal, resulting in excessive levels of mature postmitotic neutrophils and their immediate precursors. In later stages of the disease (accelerated phase and blast crisis), differentiation is reduced and expansion of immature progenitors is observed. Additional clonal karyotypic abnormalities are typically only observed during the accelerated and blast crisis phases (Hehlmann et al., 2020). CML has one of the simplest cancer genomes. It is driven by a single genetic abnormality arising somatically in an HSC, the Philadelphia (Ph) chromosome, the result of a balanced translocation between chromosomes 9 and 22 that creates a fusion of the genes for BCR and ABL1. The product of the BCR-ABL1 fusion gene is a dysregulated cytoplasmic protein-tyrosine kinase, BCR-ABL1. CML thus represents a natural model of dysregulated granulocytopoiesis (Quintás-Cardama and Cortes, 2009). Cell biological studies have shown that Ph+ cells expressing markers of normal HSC are capable of engrafting immunodeficient mice (Sirard et al., 1996; Lewis et al., 1998), implying that these cells are leukemia-initiating or leukemic ‘stem’ cells (LSCs). More mature committed progenitors in CML, like normal progenitors, lack sustained self-renewal capacity and cannot stably engraft immunodeficient mice nor generate hematopoietic colonies in vitro upon serial replating (Huntly et al., 2004). The proportion of LSCs in the BM is highly variable across CML patients at diagnosis and can range from a few percent to nearly 100% (Petzer et al., 1996; Diaz-Blanco et al., 2007; Abe et al., 2008; Thielen et al., 2016), perhaps reflecting different periods of time patients spend in chronic phase before they are diagnosed, different rates of disease progression, or both. There is persuasive experimental evidence of significant feedback regulation of different cell compartments in the dynamics of myeloid cell production in both normal and CML hematopoiesis, including signaling between the normal and CML cells (Jiang et al., 1999; Devireddy et al., 2005; Vicente-Dueñas et al., 2009; Naka et al., 2010; Reynaud et al., 2011; Zhang et al., 2012; Krause et al., 2013; Walenda et al., 2014; Welner et al., 2015). For instance, experiments in a mouse model of CML provided evidence that IL-6, produced by leukemic neutrophils, blocked MPP differentiation toward a lymphoid fate, implying feedback from the myeloid lineage onto MPPs (Reynaud et al., 2011). Surprisingly, our knowledge of the details of feedback regulation in hematopoiesis is still incomplete, especially for granulopoiesis, where even late-stage feedback interactions are poorly understood. For example, two cytokines, granulocyte colony-stimulating factor (G-CSF) and granulocyte-macrophage colony-stimulating factor (GM-CSF), can pharmacologically increase neutrophil production, but mice lacking both cytokines maintain baseline neutrophil levels and can still increase neutrophil production in response to infection (Basu et al., 2000). In many cases, it is not known which cell types are providing and receiving the feedback, what signals are used, and what aspects of proliferative cell behavior they influence (i.e., proliferation rates, renewal probability, or progeny fate choice). In spite of this knowledge deficit, CML can be treated quite effectively using selective small-molecule tyrosine kinase inhibitors (TKIs) of the BCR-ABL1 kinase. TKIs such as imatinib, dasatinib, and nilotinib, which inhibit proliferation and increase apoptosis of Ph+ cells, have dramatically lowered CML death rates (Gambacorti-Passerini et al., 2011). The response to TKI therapy in CML is monitored primarily by determining the level of BCR-ABL1 mRNA transcripts in peripheral blood, normalized to a control RNA and expressed as a percentage on an International Scale (Arora and Press, 2017). BCR-ABL1 transcript levels, an approximation of the proportion of circulating malignant cells at any given time, generally decrease exponentially in patients responding to TKI therapy resulting in at least two distinct slopes when plotted semi-logarithmically—an initial rapid decline attributed to TKI-induced killing of more mature myeloid cells, and a subsequent slower decline postulated to represent lower death rates of more primitive leukemic stem/progenitor cells (Michor et al., 2005). Clinical resistance to TKI therapy in CML is a significant problem and is classified as acquired resistance (increasing BCR-ABL1 transcript levels following a substantial decrease) or primary resistance (lack of an adequate initial response). Many patients with acquired resistance have developed mutations in the BCR-ABL1 kinase domain that mediate pharmacological resistance to the TKI (Ernst and Hochhaus, 2012). By contrast, 10–15% of newly diagnosed CML patients fail to achieve an ‘early molecular response,’ defined as the level of BCR-ABL1 transcripts being less than 10% at 3 mo (Hanfstein et al., 2012; Marin et al., 2012). Clinical data indicate that switching TKIs may not benefit these patients (Yeung et al., 2012; Yeung et al., 2015), suggesting that this group is destined to do poorly regardless of the specific inhibitor used. BCR-ABL1 mutations are generally not present in this group of patients (Zhang et al., 2009; Pietarinen et al., 2016), and thus the mechanism(s) underlying this primary resistance is unclear. We hypothesized that these variable patient responses to TKI therapy arise from nonlinearity introduced by non-cell-autonomous interactions between normal and CML cells. To test this hypothesis, we developed a novel mathematical model of CML hematopoiesis and TKI treatment that incorporates lineage branching and interactions between normal and CML cells through feedback and feedforward regulation. Mathematical modeling of leukemia has a long history aimed at understanding disease progression and improving treatment response using single and combination targeted therapies and immunotherapy (Whichard et al., 2010; Pujo-Menjouet, 2015; Brunetti et al., 2021; Kuznetsov et al., 2021; Roeder and Glauche, 2021). Further, recent efforts have been made to integrate mathematical modeling in clinical decision-making to design personalized therapies (Hoffmann et al., 2020; Engelhardt and Michor, 2021). Many models of leukemia have utilized simplified lineage architectures and minimal feedback (Roeder et al., 2006; Komarova and Wodarz, 2007; Horn et al., 2008; Foo et al., 2009; Hähnel et al., 2020; Pedersen et al., 2021). While these models can be made to fit the multiphasic disease response data of CML to TKI treatment, the simplicity of the models can make these fitted parameters of limited clinical value. More physiologically accurate, nonlinear models that account for cell–cell signaling and lineage branching are expected to improve clinical relevance. Mathematical models that incorporate feedback signaling have been developed in normal (Engel et al., 2004; Marciniak-Czochra et al., 2009; Mahadik et al., 2019, Mon Père et al., 2021) and diseased (Wodarz, 2008; Sachs et al., 2011; Krinner et al., 2013; Stiehl et al., 2014; Stiehl et al., 2015; Crowell et al., 2016; Woywod et al., 2017; Jiao et al., 2018; Stiehl et al., 2018; Zenati et al., 2018; Park et al., 2019; Sharp et al., 2020) hematopoiesis. Because of the vast number of possible ways in which feedback models of normal hematopoiesis and leukemia can be configured, mathematical models tend to greatly simplify the lineage architectures and the feedback interactions among the cell types. For example, Manesso et al., 2013 developed a hierarchical ordinary differential equation (ODE) model of normal hematopoiesis containing multiple cell types and branch points (16 cell types and 4 branch points) in the lineage tree. Limiting the feedback loops to involve only local, negative regulation (e.g., regulation by self and immediate progenitor/progeny in the lineage tree) results in about 106 models, which enabled the use of a stochastic optimization algorithm to obtain parameters consistent with homeostasis and a requirement for a rapid return to equilibrium following system perturbations. In the context of leukemia, the model architectures are typically much simpler. Generally, models of leukemia introduce a parallel mutant lineage with the same structure as that used to model the normal hematopoietic cells but with different parameters. For example, Wodarz developed an unbranched lineage ODE model of normal and leukemia stem and differentiated cells in which feedback from the differentiated cells controlled whether the stem cells divided symmetrically or asymmetrically, and demonstrated this provides a mechanism for blast crisis in CML to occur without additional mutations (Wodarz, 2008). Krinner et al. incorporated positive and negative feedback regulation of differentiation and proliferation in an unbranched lineage model that combined a discrete agent-based model for the stem cell compartment with an ODE system for the progenitor and differentiated cells to provide a detailed view of the stem cell dynamics and to test the effect of therapies (Krinner et al., 2013). Stiehl et al., 2015 developed an unbranched lineage ODE model of normal and leukemic cells in which only negative feedback regulation of stem and progenitor cell self-renewal fractions was considered, and this was further limited to arise only from factors produced by the postmitotic, mature normal and leukemic cells. Later work extended this approach to investigate clonal selection and therapy resistance (Stiehl et al., 2014), the role of cytokines on leukemia progression (Stiehl et al., 2018), combination treatment strategies (Banck and Görlich, 2019), and niche competition (Stiehl et al., 2020). Clonal competition was also considered in an ODE feedback model of CML (Woywod et al., 2017) and a stochastic model with feedback (Dinh et al., 2021). Simpler unbranched lineage models of normal and leukemic cells in which only the normal cells respond to feedback but normal and leukemic cells compete for space in the BM have been used to investigate regimes of coexistence of normal and leukemic cells (Crowell et al., 2016; Jiao et al., 2018) and design combination therapies using optimal control algorithms (Sharp et al., 2020). Here, we develop a nonlinear ODE model of normal and CML hematopoiesis using a general approach that integrates an automated method, design space analysis (DSA; Fasani and Savageau, 2010), with data gleaned from previously published experiments, and from two new in vivo experiments presented here that separately decrement the number of stem cells and terminally differentiated myeloid cells in the BM of mice. This approach enables us to systematically select among plausible model architectures and signaling interactions without a priori knowledge of which cells are providing and receiving signaling stimuli. We start with a model for normal hematopoiesis that accounts for stem, multipotent progenitor cells, and two types of terminally differentiated cells representing the myeloid and lymphoid lineage branches. This approach allows us to reduce the potential model space from about 60,000 models to a single model class and reveals the existence of feedforward and feedback mechanisms. We then extend the model to incorporate CML hematopoiesis by introducing a parallel lineage of CML cells with the same model architecture but with different parameters. The model captures the dynamics of CML at various stages of the disease and exhibits variable response to TKI treatment consistent with that observed in clinical data. The model suggests biomarkers of primary resistance, identifies the underlying mechanisms governing the response to TKI therapy, and suggests new treatment strategies.

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 (TDm) and lymphoid (TDl). The HSC self-renew with fraction (e.g., probability) p0 or differentiate with fraction 1-p0. That is, the fraction of HSC that remain as HSC after division is p0. The MPPs self-renew with fraction p1 and differentiate into either lymphoid or myeloid cells with fractions q1 and 1-p1-q1, respectively. The HSC and MPPs divide with rates η1 and η2 and the myeloid and lymphoid cells die at rates dm and dl, 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. 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 (TDm) with the ability to negatively regulate the fraction q1 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 (η1) and inhibits the division of MPPs through a negative feedforward loop on η2.

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 p1. 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 ~106 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 (γ5 > 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 γ1 and γ3 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 γ1 and γ3 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 γ1 and γ3 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 γ1, 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. 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 (HSCL; SL), as indicated by p0L 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 HSCL 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 HSCL 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. 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 HSCL self-renewal fraction and the feedback gain on the HSCL 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 to be less than or equal to , the maximal self-renewal fraction of the normal HSC, motivated by the paucity of evidence that is larger than , coupled with experimental data suggesting that is less than or equal to . 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 and 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 and (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 sufficiently close to and should be sufficiently small. As the ratio / decreases from 1, the system requires smaller feedback gains 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 / is sufficiently large or when is sufficiently small. To further examine these biological constraints, we calculated characteristic effective self-renewal fractions for normal and leukemic stem cells, defined as and , where , 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 HSCL and HSC self-renewal fractions: / . Here, all eligible parameter sets representing the states of the normal system are considered and the leukemic parameters / and 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 HSCL or leukemic MPPs (MPPL) were implanted into sublethally irradiated mice (Figure 5A). Transplantation of HSCL enables engraftment and myeloid cell production that leads to CML. On the other hand, transplanting MPPLs 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 MPPL 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 HSCL, MPPL, and differentiated myeloid and lymphoid cells and tracked the outcomes when 4000 HSCL or MPPL 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 HSCL transplant and a simple majority of lymphoid cells for MPPL 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 HSCL, MPPL, 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 HSCL are transplanted (solid curves), the donor-derived MPPL (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 MPPL are introduced (dashed curves), their population decreases (Figure 5C, right) because the MPPL do not stably engraft. Concomitantly, there is burst of donor-derived myeloid and lymphoid cells at early times (Figure 5C, right) as the transplanted MPPL differentiate. The early time dynamics of the myeloid and lymphoid cells depend on the specific values of the MPPL self-renewal (p1) and fate control (q1) fractions, whose values in turn depend on the number of myeloid cells through negative feedback regulation. In particular, if , 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 . In both cases, because the MPPLs 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 HSCL 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 HSCL 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 HSCL are one half as sensitive to negative feedback regulation compared to the normal HSC (). 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 HSCL and MPPL proportional to their proliferation rates, with the HSCL 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 HSCL fractions (HSCL, Figure 6A), the numbers of MPPL (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 HSCL due to the loss of negative feedback from the MPPL. This loss of negative feedback from the MPPL 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 HSCL decreases their division rates due to the autocrine negative regulatory loop as well as the division rates of the MPPs and MPPL 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 HSCL and MPPL gradually decrease in response to TKI-induced cell death, which drives an accompanying decrease in the leukemic differentiated cells. A small, transient increase in MPPL is observed before the gradual decline. This is driven primarily by the increase in flux into the MPPL compartment by differentiating HSCL, although there is also a small contribution from the feedforward regulation of the MPPL division rate, which lowers the effectiveness of TKI therapy on the MPPL. Both of these effects are reduced as the HSCL 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. At intermediate treatment time with larger (90–99%) HSCL 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 HSCL is much larger than the previous case because there are fewer normal cells to compete with in the BM. This significantly decreases the HSC/HSCL and MPP/MPPL division rates through the negative feedback/feedforward regulation and correspondingly the rates of TKI-induced cell death. Accordingly, at later times the MPPL population rebounds, driven by the flux of differentiating HSCL, and eventually the system reaches a state in which both normal and leukemic cells coexist. The stem cell compartment is dominated by HSCL 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 HSCL 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 MPPL (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 MPPL 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 MPPL, 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 HSCL fraction, intermediate—18 mo after CML initiation ~99% initial HSCL fraction, late—36 mo after CML initiation ~99.99% initial HSCL fraction). Our results suggest that the higher the HSCL fraction at the start of therapy, the less effective the therapy. This follows from the feedback/feedforward regulation where increased numbers of HSC and HSCL decrease their own proliferation rates as well as those of the MPPs and MPPL (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 HSCL burden (94 ± 1.5% of the HSC population) or an intermediate HSCL 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 HSCL showed a decrease in the percentage of circulating BCR-ABL1+ granulocytes in response to TKI therapy, mice with the highest HSCL 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 HSCL 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 HSCL 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, 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 / 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 . We thus varied the HSCL parameters accordingly. In Figure 7B, we plot the results for as a bivariate histogram for and with proportion of response (blue) and nonresponse (orange) for every parameter set. Surprisingly, we found that although the fitness impacts response, the parameter that clearly distinguished responders from nonresponders was the maximal self-renewal fraction p0,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 . In particular, larger values of p0,max and γ1 are correlated with a decreased response to TKI therapy after leukemia develops. Although these parameters are associated with normal HSC, the self-renewal fraction p0L of HSCL and feedback strength depends on these parameters since we assumed the fitness of the CML stem cells 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. 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 p0L (the fraction of HSCL 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 p0,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: . 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 with sensitivity of and specificity of (Figure 7E, orange curves) compared to the optimal sensitivity and specificity of the velocity-based prognostic ( and respectively) and and 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 () from nonresponders (), 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. 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 p0,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 p0,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 p0,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 MPPL, 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 HSCL 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 HSCL, 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.

In this work, we developed a nonlinear mathematical model of normal and CML hematopoiesis that incorporated feedback control, lineage branching, and signaling between normal and CML cells. Using ODEs, we modeled the dynamics of the stem, multipotent progenitor, and terminally differentiated cell populations. To filter through the combinatorial explosion of models that occurs when cell–cell signaling interactions are taken into account, we focused first on normal hematopoiesis. We used DSA (Savageau et al., 2009; Fasani and Savageau, 2010; Lomnitz and Savageau, 2013; Lomnitz and Savageau, 2016), an approach that enables models to be distinguished based on their range of qualitatively distinct behaviors without relying on knowledge of specific values of the parameters, to perform an automated search for regions of stability in thousands of proposed models and efficiently eliminate unphysiological, unstable models. When combined with previous observations and new in vivo data to further constrain cell–cell interactions, we arrived at a new feedback-feedforward model (Figure 2F). Using cell perturbation experiments in mice, we validated several features of the model, including feedback from differentiated myeloid cells on MPP self-renewal, and feedforward regulation by stem cells on proliferation of stem and MPP compartments. We postulate that these regulatory loops may also regulate human blood cell production. While there are some known differences between mouse and human hematopoiesis (Parekh and Crooks, 2013), many signaling pathways are conserved between species. For example, the role of IL-6 in regulating lymphoid differentiation (e.g., γ4 in Figure 2F) has been validated in mice (Reynaud et al., 2011) and human samples (Welner et al., 2015), while CCL3 mediates negative feedback from progenitors onto stem cell self-renewal (e.g., γ1 in Figure 2F) in both mice (Staversky et al., 2018) and humans (Broxmeyer et al., 1989). TGF-β, produced by HSC, differentiated myeloid cells, and BM stroma, is a candidate factor regulating negative feedback of HSC onto their own division rate and that of the MPPs (Zhao et al., 2014b; Naka and Hirao, 2017), while IL-6 may inhibit MPP self-renewal and increase myeloid differentiation (Zhao et al., 2014a), at least under stress conditions. The role of these candidate hematopoietic regulators could be tested directly in our mouse model via a genetic approach. Moving forward, it will be important to validate results from mouse model systems in human studies whenever possible. We used a grid-search algorithm to determine a set of approximately 500 biologically relevant parameter sets for our new model. While we could have used other approaches such as the Latin hypercube algorithm to sample our multidimensional parameter space (Read et al., 2018), we chose to perform a gridsearch because of the ease of implementation and the fact that our goal was not to exhaustively search the full parameter space, but rather to obtain a set of biologically relevant parameter values consistent with normal hematopoietic homeostasis. In particular, using each parameter set in the model yields steady states that are consistent with normal ranges of hematopoietic cells. These parameter sets model a population of individuals with normal cell counts but person-to-person variability of parameters due, for example, to genetic, epigenetic and/or environmental differences. We then extended the model to incorporate CML hematopoiesis by introducing a mutant lineage with the same structure as the normal system. We incorporated one of the central features of CML pathophysiology, that the leukemic stem cell clone, hypothesized to arise from a single HSC that acquires a Ph chromosome, has a competitive advantage over normal HSC and with time comes to dominate the stem cell compartment (Dingli et al., 2010; Thielen et al., 2016; Holyoake and Vetrie, 2017; Majeti et al., 2022). This competitive advantage could be a consequence of positive feedback (autocrine or paracrine) on the HSCL population or negative feedback with different strengths for normal and leukemic stem cells. Candidate mediators of such positive and negative feedback include interleukin-3 (Jiang et al., 1999) and CCL3 (Baba et al., 2016), respectively. Our current model incorporated differential negative feedback of MPPs on HSC (Figure 4A) with the HSCL being less sensitive to the negative feedback than are the normal HSC, which is consistent with CCL3 (Eaves et al., 1993, Baba et al., 2013). This one difference provided leukemic cells with a competitive advantage for growth, and in the absence of treatment, the leukemic cells will take over the BM at the expense of normal cells (Figure 4B). Upon exploration of the leukemic parameter space, we found that only the leukemic cell parameters for the leukemic stem cells (HSCL)—the maximal HSCL self-renewal fraction , the feedback gain on the HSCL self-renewal fraction, and the TKI-induced HSCL death rate —have the potential to significantly influence the results. The results are insensitive to changes in the other leukemic cell parameters (see Figure 4—figure supplement 6, Figure 6—figure supplements 2–4, and Appendix 1—figure 11). When combined with TKI therapy, the feedback/feedforward model exhibited variable responses to TKI treatment, consistent with those observed in CML patients. That is, although our 500 parameter sets were consistent with normal hematopoietic cell counts, the responses to TKI treatment were highly variable, with some sets responding to treatment while others did not. The model predicted that a contributor to primary TKI resistance is the overall proportion of HSC that are leukemic, consistent with experimental data in mice (Figure 6G) as well as patient data (Thielen et al., 2016). However, leukemic stem cell burden alone does not predict the molecular response to TKIs, as observed both clinically (Thielen et al., 2016) and in our data (Figure 7A), since some patients with high HSCL fractions in their BM nonetheless still respond to TKIs. The model suggested that a key predictor of reduced response to TKI treatment is an increased tendency of normal hematopoietic stem cells to self-renew, which in turn influences self-renewal of the leukemic stem cells since they were estimated to be sufficiently fit with respect to the normal stem cells. This is also consistent with clinical data that suggest that CML patients whose normal and leukemic cells share mutations in genes such as TET2 and ASXL1, which are known to increase stem cell self-renewal (Steensma, 2018), tend to have inferior outcomes under TKI therapy (Kim et al., 2017; Marum et al., 2017). This is illustrated in Figure 8F (red panel), where the high initial HSCL population and the subsequent decline of progenitor cells reveals the effect that high stem cell self-renewal has on driving TKI resistance. In our model, the presence of a TET2 or ASXL1 mutation in both normal and leukemic stem cells that led to a proportional increase in self-renewal in both populations would tend to cause resistance to TKI therapy, provided that the HSCL are sufficiently fit in the presence of the mutations, which we would expect. The self-renewal-driven resistance we describe herein challenges the prevailing paradigm that TKI resistance is proliferation-driven and a consequence of HSCL quiescence (Graham et al., 2002; Corbin et al., 2011). Because stem cell self-renewal is hard to quantify experimentally, we developed a clinical prognostic criterion to predict TKI response based on the relative changes in the BCR-ABL1 transcripts over a 3-month period. Using the synthetic data from our normal and leukemic parameter sets, we found that using changes in transcripts from 3 to 6 mo was very effective in predicting the long-term outcome of treatment (e.g., reaching MR3 within 18 mo). In contrast, using transcript data from 0 to 3 mo resulted in less accurate predictions. This observation also holds for prognostic criteria based on EMR and transcript halving time, which are currently used in the clinic. We then tested the prognostic criteria on data obtained from small number of anonymized CML patients and found the same conclusions hold. Our results suggest that the relative change prognostic criterion more accurately predicts patient response than EMR and the halving time, although more data are needed to confirm this. Our cohort of patients was small due to the variable nature of patient treatment and inconsistent data collection, for example, patients were frequently switched from one TKI to another or one dosage to another (sometimes multiple times), and the patients’ BCR-ABL1 transcript levels were not always consistently recorded. However, we believe that this pilot study demonstrates the feasibility of our approach. Moving forward, we aim to apply our approach to larger datasets and hope to convince others to do the same. Two strategies can be postulated to overcome TKI resistance. One approach could be to decrease stem cell self-renewal either by inhibiting self-renewal directly (e.g., by augmenting TET2 function using ascorbate; Agathocleous et al., 2017; Cimmino et al., 2017) or by promoting differentiation (e.g., using retinoids; Drumea et al., 2008). By applying combined TKI and pro-differentiation therapy, the self-renewal fractions of the normal and leukemic stem cells can be decreased and brought closer together, which ultimately disadvantages the leukemic cells because of TKI-induced cell death (Figure 8F, blue panel). An alternative or complementary approach would be to increase stem cell proliferation via pro-proliferative stimuli such as IFN-alpha (Essers et al., 2009) to increase efficacy of TKIs in killing HSCL. It is apparent that the feedback/feedforward interactions incorporated in our model, which are necessarily somewhat restricted, may be further constrained by spatial characteristics of the BM microenvironment. Nonetheless, our model still displays consistent and biologically relevant behaviors, and although further refinement of the model behaviors is possible, based upon our findings the key behaviors (feedback mechanisms and importance of stem cell self-renewal) would be expected to remain much the same. To explore experimentally observed phenomena not captured by our current model such as treatment-free remission, where a low level of HSCL persists in the absence of TKI pressure without myeloid cell expansion, improvement of the model is necessary. For example, it may be necessary to incorporate features of the BM microenvironment such as stem cell–niche interactions (MacLean et al., 2014; Lai et al., 2022) and interactions with immune cells (Hähnel et al., 2020). The inclusion of a quiescent stem cell state and additional cellular compartments (such as committed progenitors) coupled with appropriately constrained cell–cell signaling would also make the model more physiologically accurate. In summary, the feedback/feedforward model we have presented here, while a simplified version of normal and CML hematopoiesis, makes novel and testable predictions regarding the origins of non-genetic primary resistance, which patients will respond to TKI treatment and suggests a combination therapy that can overcome primary resistance. Although preliminary evidence was presented to support model predictions, future work should focus on designing targeted experiments and collecting patient outcomes to generate data to more thoroughly test the model.

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 (η1 and η2, respectively). The S cells have the ability to self-renew with fractions (p0) or differentiate (1 p0). The P cells have the ability to self-renew with fraction (p1) or differentiate into either TDl (lymphoid) or TDm (myeloid) cells (q1 or 1-p1-q1, 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 TDl and TDm cells. TDm and TDl cells are postmitotic and die at rates dm and dl, 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 104 < HSC < MPPs with MPPs fixed at 105 and MPPs < TDl < TDm. For 106 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 γ5 > 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 HSCL or MPPL 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 HSCL depletions varied between 50 and 70% and the MPPL 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 HSCL, 35% for MPPLs. 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: .

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 p0 by including a new constant repressive force that affects both normal and leukemic self-renewal, resulting in and , where and 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).