Predictive nonlinear modeling of malignant myelopoiesis and tyrosine kinase inhibitor therapy
Abstract
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 BCRABL1 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 BCRABL1 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 BCRABL1 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 selfrenewal 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.
Editor's evaluation
This is an important study that investigates the impact of tyrosine kinase inhibitors (TKIs) in chronic myeloid leukemia. Through a combination of preclinical in vivo measurements, clinical data, and computational modeling, the authors present solid evidence regarding the heterogeneous effects of TKIs in patients and how the response to treatment may be improved. This study is of interest to those working in the fields of mathematical oncology and cancer biology.
https://doi.org/10.7554/eLife.84149.sa0Introduction
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], granulocytemacrophage progenitors [GMPs], and megakaryocyteerythroid 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 BCRABL1 fusion gene is a dysregulated cytoplasmic proteintyrosine kinase, BCRABL1. CML thus represents a natural model of dysregulated granulocytopoiesis (QuintásCardama 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 leukemiainitiating or leukemic ‘stem’ cells (LSCs). More mature committed progenitors in CML, like normal progenitors, lack sustained selfrenewal 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; DiazBlanco 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; VicenteDueñ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 IL6, 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 latestage feedback interactions are poorly understood. For example, two cytokines, granulocyte colonystimulating factor (GCSF) and granulocytemacrophage colonystimulating factor (GMCSF), 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 smallmolecule tyrosine kinase inhibitors (TKIs) of the BCRABL1 kinase. TKIs such as imatinib, dasatinib, and nilotinib, which inhibit proliferation and increase apoptosis of Ph^{+} cells, have dramatically lowered CML death rates (GambacortiPasserini et al., 2011). The response to TKI therapy in CML is monitored primarily by determining the level of BCRABL1 mRNA transcripts in peripheral blood, normalized to a control RNA and expressed as a percentage on an International Scale (Arora and Press, 2017). BCRABL1 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 semilogarithmically—an initial rapid decline attributed to TKIinduced 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 BCRABL1 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 BCRABL1 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 BCRABL1 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. BCRABL1 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 noncellautonomous 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; PujoMenjouet, 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 decisionmaking 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; MarciniakCzochra 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 10^{6} 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 agentbased 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 selfrenewal 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.
Results
Model of normal hematopoiesis
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 selfrenew with fraction (e.g., probability) p_{0} or differentiate with fraction 1p_{0}. That is, the fraction of HSC that remain as HSC after division is p_{0}. The MPPs selfrenew with fraction p_{1} and differentiate into either lymphoid or myeloid cells with fractions q_{1} and 1p_{1}q_{1}, respectively. The HSC and MPPs divide with rates η_{1} and η_{2} 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 selfrenewal 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 coworkers (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 righthand 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 (Ssystem) of the model. The number of Ssystems in each model depends on the number of combinations of positive and negative terms in the rates of change. If the equilibria of the Ssystems, which are determined analytically, are not selfconsistent (e.g., consistent with the assumed dominance of terms reflected in the inequalities) or the equilibria are not stable, then the Ssystem is rejected. If all the Ssystems of a particular model are rejected, then that model is removed from further consideration. Models with at least one selfconsistent and stable Ssystem 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 selfrenewal fraction but differ by where this regulation arises. The models within the classes share at least one Ssystem 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. Interleukin6 (IL6) 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 α [MIP1α]), produced in BM by basophilic myeloid progenitors (Baba et al., 2016), acts to inhibit the proliferation and selfrenewal 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 selfrenewal fraction are regulated. These above results suggest that IL6 is a candidate feedback factor expressed in the myeloid compartment (TD_{m}) with the ability to negatively regulate the fraction q_{1} of MPPs that differentiate into lymphoid cells. CCL3 is a candidate factor mediating negative feedback from the MPP population onto HSC selfrenewal. 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.
Depletion of HSC increases HSC and MPP proliferation
As described in ‘Methods,’ healthy C57BL6/J (B6) mice were treated with lowdose (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 postirradiation, 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}.
Depletion of mature myeloid cells increases the MPP population
B6 mice were treated with the antigranulocyte 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 selfrenewal fraction p_{1}.
Taking all these results into consideration, we arrive at the feedbackfeedforward 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.
Parameter estimation for feedbackfeedforward model of hematopoiesis
To determine biologically relevant parameters for the feedbackfeedforward model in Figure 2F, a gridsearch algorithm was employed. The full ODE model is given in ‘Methods’ and Appendix 1 (Section 2). The 12 model parameters (proliferation and death rates, selfrenewal 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^{6} 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.
Sensitivity analyses of hematopoiesis model
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 selfrenewal 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 Ssystems. 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 darkgreen and lightgreen 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% (dotdashed), 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.
Extension of the hematopoiesis model to CML
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 nonmutant 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_{0}^{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 selfrenewal 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.
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 selfrenewal parameters—the maximal HSC^{L} selfrenewal fraction ${p}_{0,max}^{L}$ and the feedback gain ${\gamma}_{1}^{L}$ on the HSC^{L} selfrenewal 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 selfrenewal 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 longterm culture initiating cells (LTCIC; thought to be phenotypically similar to stem cells) decrease significantly in in vitro cultures while the number of normal LTCIC is unchanged, consistent with a relative decrease in selfrenewal probability of the CML cells (Udomsakdi et al., 1992). In vivo, HSC selfrenewal 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 BCRABL1 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 selfrenewal capacity of BCRABL1+ stem cells.
Next, we performed a sweep through leukemic stem cell selfrenewal 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 selfrenewal 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 selfrenewal 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.
Validation of the CML model
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^{L}s does not allow for longterm engraftment but results in a larger fraction of donorderived lymphoid cells after 35 d (Figure 5B). This study presented evidence that IL6 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 donorderived 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 donorderived 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 donorderived myeloid and lymphoid cells at early times (Figure 5C, right) as the transplanted MPP^{L} differentiate.
The early time dynamics of the myeloid and lymphoid cells depend on the specific values of the MPP^{L} selfrenewal (p_{1}) and fate control (q_{1}) fractions, whose values in turn depend on the number of myeloid cells through negative feedback regulation. In particular, if $1{p}_{1}>2{q}_{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 $1{p}_{1}<2{q}_{1}$ . In both cases, because the MPP^{L}s do not stably engraft and instead differentiate into lymphoid and myeloid cells, we observe that there is a larger fraction of donorderived 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 longerlived (smaller death rate) than the shorterlived myeloid cells, which have a larger death rate.
Leukemic stem cell load influences TKI therapy outcomes
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; DiazBlanco 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 ($\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 longterm 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 BCRABL1 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} (bluedashed), leukemic terminally differentiated lymphoid (lightgreendashed), and myeloid (darkgreendashed) 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 TKIinduced 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 preleukemic equilibrium values.
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 TKIinduced 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, BCRABL1 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 selfrenewal 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 shortlived, 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 BCRABL1 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 BCRABL1 transcripts, which decreases below 10^{–1}, representing a socalled major molecular response (MMR or MR3), which represents a major goal of therapy in CML as the risk of relapse and leukemiarelated death is virtually nonexistent once this milestone is achieved (Hochhaus et al., 2017). Consistent with previous interpretations, the rapid initial decrease in BCRABL1 transcripts is due to TKIinduced 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 longterm, 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 BCRABL1 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 BCRABL1 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 TKImediated 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 BCRABL1 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).
HSC^{L} load influences the response to TKI therapy in a mouse CML model
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 (BCRABL1^{+}) HSC by transplantation of BM from conditional BCRABL1 transgenic mice (Koschmieder et al., 2005) into unirradiated congenic recipient mice. Following stable engraftment, BCRABL1 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 BCRABL1^{+} 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 BCRABL1^{+} HSC at the start of treatment.
HSC selfrenewal as an additional determinant of TKI response
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 CCL3mediated inhibition of selfrenewal (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 BCRABL1 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 $\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 $\frac{{\overline{{p}_{0}}}^{L}}{{p}_{0,max}}\ge 0.7$ as a bivariate histogram for $\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 $\frac{{\overline{{p}_{0}}}^{L}}{{p}_{0,max}}$ impacts response, the parameter that clearly distinguished responders from nonresponders was the maximal selfrenewal fraction p_{0,max} for normal stem cells, shown in Figure 7B (marginal yaxis). 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 $\frac{{\overline{{p}_{0}}}^{L}}{{p}_{0,max}}$ . In particular, larger values of p_{0,max} and γ_{1} are correlated with a decreased response to TKI therapy after leukemia develops. Although these parameters are associated with normal HSC, the selfrenewal fraction p_{0}^{L} of HSC^{L} and feedback strength ${\gamma}_{1}^{L}$ depends on these parameters since we assumed the fitness of the CML stem cells $\frac{{\overline{{p}_{0}}}^{L}}{{p}_{0,max}}$ is larger than a minimum threshold. Therefore, increasing the selfrenewal 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 p_{0}^{L} (the fraction of HSC^{L} selfrenewal 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 selfrenewal fraction for both normal and leukemic stem cells because of the release of negative feedback. When the maximum selfrenewal p_{0,max} is larger, the leukemic stem cells experience an acute increase in selfrenewal, 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 selfrenewal (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 selfrenewal fare worse when their CML is treated using TKIs than patients with lower stem cell selfrenewal.
Predicting longterm response to TKI treatment
Several measures of the response of CML patients to TKI therapy have been developed, based on BCRABL1 transcript levels in peripheral blood. Here we test a new, modeldriven 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 socalled ‘early molecular response’ or EMR (defined as BCRABL1 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 selfrenewal 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 BCRABL1 transcript levels (Figure 7D). It is important to be able to predict the longterm 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–6month 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 (1FPR, with FPR being the false positive rate) resulted in a prognostic threshold of $PF\approx 3.2,$ with sensitivity of $\approx 0.91$ and specificity of $\approx 0.91$ (Figure 7E, orange curves) compared to the optimal sensitivity and specificity of the velocitybased prognostic ($\approx 0.89$ and $\approx 0.88,$ respectively) and $\approx 0.77$ and $\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 ($<3.2$) from nonresponders ($>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–3month 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–6month TKI response in predicting the longterm 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–6month 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–6month 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 6month period, which were added to the data from Figure 7F. The aggregated data (Figure 7G) reaffirms the improved accuracy in prediction using the 3–6month 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 velocitybased 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.
Improving response to therapy: Combining TKIs with interventions that promote differentiation
Our model suggests that combination therapy to modulate the stem cell selfrenewal rate, in addition to directly targeting the leukemic HSC and MPPs with TKI therapy, might counteract TKI treatment resistance mediated by high stem cell selfrenewal. Such prodifferentiation therapy could be accomplished through either direct stimulation of differentiation or through suppression of selfrenewal. In our modeling experiments, we explored the impact of this approach through the suppression of selfrenewal (see ‘Methods’ for details). To begin the exploration of the combined TKIdifferentiation therapy, we performed this combination therapy on each of our 478 parameter sets, which represents a population of CML patients with persontoperson 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 selfrenewal (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 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 selfrenewal is already close to the ideal effective selfrenewal 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 BCRABL1 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 TKIinduced 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 reequilibrates. 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 selfrenewal and corresponding feedback gains, the differentiation promoter acts to bring the selfrenewal 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.
Discussion
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 feedbackfeedforward model (Figure 2F).
Using cell perturbation experiments in mice, we validated several features of the model, including feedback from differentiated myeloid cells on MPP selfrenewal, 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 IL6 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 selfrenewal (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 IL6 may inhibit MPP selfrenewal 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 gridsearch 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 persontoperson 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 HSC^{L} population or negative feedback with different strengths for normal and leukemic stem cells. Candidate mediators of such positive and negative feedback include interleukin3 (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 HSC^{L} 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 (HSC^{L})—the maximal HSC^{L} selfrenewal fraction ${p}_{0,max}^{L}$ , the feedback gain ${\gamma}_{1}^{L}$ on the HSC^{L} selfrenewal fraction, and the TKIinduced HSC^{L} death rate $TK{I}_{HSCL}$—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 HSC^{L} 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 selfrenew, which in turn influences selfrenewal 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 selfrenewal (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 HSC^{L} population and the subsequent decline of progenitor cells reveals the effect that high stem cell selfrenewal 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 selfrenewal in both populations would tend to cause resistance to TKI therapy, provided that the HSC^{L} are sufficiently fit in the presence of the mutations, which we would expect. The selfrenewaldriven resistance we describe herein challenges the prevailing paradigm that TKI resistance is proliferationdriven and a consequence of HSC^{L} quiescence (Graham et al., 2002; Corbin et al., 2011).
Because stem cell selfrenewal is hard to quantify experimentally, we developed a clinical prognostic criterion to predict TKI response based on the relative changes in the BCRABL1 transcripts over a 3month 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 longterm 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’ BCRABL1 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 selfrenewal either by inhibiting selfrenewal 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 prodifferentiation therapy, the selfrenewal fractions of the normal and leukemic stem cells can be decreased and brought closer together, which ultimately disadvantages the leukemic cells because of TKIinduced cell death (Figure 8F, blue panel). An alternative or complementary approach would be to increase stem cell proliferation via proproliferative stimuli such as IFNalpha (Essers et al., 2009) to increase efficacy of TKIs in killing HSC^{L}.
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 selfrenewal) would be expected to remain much the same. To explore experimentally observed phenomena not captured by our current model such as treatmentfree remission, where a low level of HSC^{L} 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 nongenetic 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.
Methods
Mathematical model of hematopoiesis
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; MarciniakCzochra 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 selfrenew with fractions (p_{0}) or differentiate (1 p_{0}). The P cells have the ability to selfrenew with fraction (p_{1}) or differentiate into either TD_{l} (lymphoid) or TD_{m} (myeloid) cells (q_{1} or 1p_{1}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.
Design space analysis
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 (SSystem). 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 Ssystem’s steady state to assess selfconsistency. If the inequalities are satisfied, the system is selfconsistent 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 Ssystems associated with a model do not have any equilibria that are selfconsistent 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.
Parameter estimation
To approximate biologically relevant parameters for the model a gridsearch 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 steadystate 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^{4} < HSC < MPPs with MPPs fixed at 10^{5} and MPPs < TD_{l} < T_{Dm}. For 10^{6} 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.
Modeling CML development
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 selfrenewal. 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 selfrenewal. The full equations used in the model are shown in Appendix 1, Equations 30–37.
Modeling transplant experiments
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^{L}s. 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.
Modeling TKI therapy
To account for the treatment by TKIs, additional proliferationdependent 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 BCRABL1 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 BCRABL1 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 BCRABL1 transcript levels: $\frac{BCRABL1}{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})}$ .
Modeling combined TKI and differentiation therapy
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_{0} by including a new constant repressive force that affects both normal and leukemic selfrenewal, 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.
CML patient data
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.
Mice
C57BL/6J female mice (Jackson Laboratories), 6–12 wk of age were used for irradiation and myeloid depletion experiments. Conditional BCRABL1 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.
Irradiation of mice
To achieve selective depletion of HSCs, a 50 cGy dose of irradiation from Xray source (Precision Xrad 320) was applied. Control mice did not receive irradiation. The distribution of time points at which observations were made (days 1, 3, and 7 postirradiation), 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).
Myeloid cell depletion
RB68C5, an antiGr1 antibody (catalog # BE0075, BioXCell) or isotype control (catalog # BE0090, BioXCell) was injected intravenously, 50 µg per mouse, and mice sacrificed 24 hr later.
BrdU injections
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.
Flow cytometry analysis of cell populations
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, Gr1, B220, Ter119) and stem/progenitor markers (cKit, Sca1, CD48) for 30 min. Streptavidin (SA)conjugated fluorochrome was utilized to detect biotinylated antibodies. Following fixation, permeabilization, and DNase digestion, antiBrdU antibody was used to assess BrdU incorporation. Events were acquired on FACS Arial II and analyzed with Flowjo v.10 software.
Antibodies
Monoclonal antibodies for flow cytometry were biotinylated mouse lineage panel (559971, BD Biosciences), PECF594 Streptavidin (562318, BD Biosciences), antimouse CD48 (561242, BD Biosciences), antimouse CD34 eFluor450 (48034182, eBiosciences), antimouse Sca1PE (108108, BioLegend), antimouse cKitAPC (17117182, eBiosciences), and FITC BrdU flow kit (559619, BD Biosciences).
Generation and TKI treatment of chimeric BCRABL1 mice
The full details of the CML mouse model will be published elsewhere (Jena et al., in preparation). Briefly, BM cells from conditional BCRABL1 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 BCRABL1 expression. Two months posttransplant, doxycycline was removed to allow induction of CMLlike 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 CMLlike leukemia (peripheral blood leukocytes > 20,000/μl and >40% circulating granulocytes), TKI treatment was initiated with dasatinib (25 mg/kg daily by oral gavage).
Data availabality
Code used to generate the figures, determine the SSystems and files containing parameter sets are found at the following GitHub: https://github.com/jonatdr/CML_Treatment (copy archived at Jonatdr, 2023).
Appendix 1
1 Design space analysis
1.1 Introduction
In this section, we provide details on the application of design space analysis (DSA), developed in Savageau et al., 2009; Fasani and Savageau, 2010; Lomnitz and Savageau, 2013; Lomnitz and Savageau, 2016 for chemical reaction networks, to a simplified version of the cell lineage model considered in the main text (see Figure 1). This enables us to analyze nonlinear dynamical systems near steady state to identify the regions in parameter space where common qualitative behaviors occur. Applying this analysis allows us to (1) ignore specific parameter values, (2) obtain analytical steady states, (3) reduce the search area in parameter space by searching boundaries separating regions of common behaviors, and (4) easily automate this process. In this section, we review how one can construct a design space given an ODE model. We then apply this analysis to an ODE model of cell lineages and show how this model can be used to define regions of stability in parameter space.
1.2 Boundaries of design space for general models
In order to apply DSA, the ODE must be a generalized mass action system, shown below in Equations 1 and 2,
for $i$ = 1, $n$. Here, $n$ corresponds to the number of dependent variables and $m$ corresponds to the number of independent variables. The ${\alpha}_{ik}$ and ${\beta}_{ik}$ parameters correspond to rate constants of the differential equation and $r$ corresponds to the number of associated rate constants.
DSA takes advantage of the above form by creating a system of deconstructed ODEs where one source term and one sink term in the differential equation dominate, known as a subsystem (Ssystem). The following is the generalized form of an Ssystem:
where $p$ and $q$ correspond to the number of positive and negative terms of the differential equation, respectively. We are interested in solving these solutions at steady state, and therefore solve the system by setting the time derivative to zero. We take advantage of the form shown in Equation 3 and take the log of the system:
thus, making this a linear solve in log space. In defining the Ssystems, we must make assumptions about the model and its parameters. To satisfy the Ssystems, we impose inequality constraints to satisfy the dominating source and sink terms of the Ssystems by the following:
We can then log transform these inequalities to obtain
These inequalities become our boundaries in parameter space in which qualitative behavior is shared once the ${X}_{j}$ ’s are evaluated at steady state, where the steady states are obtained from solving the loglinear Ssystem. Not all Ssystems will have a unique solution or satisfy the inequality constraints. The Ssystems that do not satisfy these constraints will be discarded. The analysis is summarized in Appendix 1—figure 1.
1.3 Analysis of a fourcell lineage model
We next apply this analysis to a lineage model with four cell types (see Appendix 1—figure 2). The model consists of two dividing cell types consisting of $S$ (HSC) and $P$ (MPP) cells with a division rate associated to the cells (${\eta}_{1}$ and ${\eta}_{2}$, respectively). The $S$ cells have the ability to selfrenew with fraction (${p}_{0}$) or differentiate ($1{p}_{0}$). The $P$ cells have the ability to selfrenew with fraction (${p}_{1}$) or differentiate into either $T{D}_{l}$ (lymphoid) or $T{D}_{m}$ (myeloid) cells (${q}_{1}$ or $1{p}_{1}{q}_{1}$, respectively). $T{D}_{m}$ and $T{D}_{l}$ cells only have the ability to die at rates ${d}_{m}$ and ${d}_{l}$ , respectively. We add negative feedback on the selfrenewal fraction of the stem cells from the differentiated cells. Appendix 1—figure 2 shows a schematic of the lineage with the parameters. The corresponding differential equations are
where ${p}_{0}={\overline{p}}_{0}/\left(1+\gamma {x}_{P}\right)$, and $\overline{{p}_{0}}$ is defined as the maximum stem cell selfrenewal fraction and $\gamma $ is the feedback strength.
We begin by rewriting the equations in the form shown in Equation 1.
Note that we introduce a new variable ${x}_{new}=1+\gamma {x}_{P}$ , which is necessary to achieve the form of Equation 1. We next find all combinations in which one source term and one sink term dominates the system of differential equations.
From Equations 14–18, we obtain 24 Ssystem combinations, a subset of which is shown in Appendix 1—table 1. We continue the analysis with Ssystem 2 from Appendix 1—table 1. Using Ssystem 2, we set the time derivatives to zero and rearrange the equations such that we obtain $A\overline{x}=\overline{b}$ :
such that $\overline{{x}_{S}}$ , $\overline{{x}_{P}}$ , ${\overline{x}}_{T{D}_{l}}$ , and ${\overline{x}}_{T{D}_{m}}$ are the steadystate solutions for the Ssystem and ${\overline{x}}_{new}$ is the solution of the newly defined variable at steady state. Solving the linear equation gives us
We can now construct the boundaries for this Ssystem. We substitute our steadystate solutions obtained above into the logged inequality constraints from Appendix 1—table 1. Thus, the inequalities in log space become
Out of the 24 possible Ssystems, only Ssystem 2 has a unique steady state that satisfies the constraints. We can plot the design space by varying $\overline{{p}_{0}}$ and $\gamma $ (see Appendix 1—figure 3). The design space shows one region where $\overline{{p}_{0}}>0.5$, a requirement for a positive steady state in the full system. The domain in parameter space corresponds to Ssystem 2 in Appendix 1—table 1.
It is possible to relate the Ssystem back to the true ODE system with classical techniques. For example, it was shown in Savageau et al., 2009; Fasani and Savageau, 2010; Lomnitz and Savageau, 2013; Lomnitz and Savageau, 2016 that the Ssystem and full ODE system have the same linear stability behavior in the parameter regime appropriate for the Ssystem. Thus, parameter sensitivity analyses of the Ssystem provide insight on the behavior of the full system.
Performing a linear stability analysis on the Ssystem that corresponds to the region in parameter space shown in Appendix 1—figure 3, we obtain the following eigenvalues:
which suggests a stable spiral or stable node depending on the value of d. We also perform the linear stability analysis on the full system and obtain the following eigenvalues:
and in the parameter space for the Ssystem, we obtain the same behavior. Using the full system, we can also plot the dynamics in the domain, which are shown as insets in Appendix 1—figure 3. When plotting the full system in Appendix 1—figure 3, we observe that the stem and progenitor populations oscillate before reaching steady state, characteristics of a stable spiral. The Ssystem and the true system eigenvalues are not identical, but using the same parameters in the appropriate Ssystem parameter space, the systems will have the same behavior. However, outside of the domain of validity of the Ssystem, we cannot conclude anything about the full system’s dynamics from the Ssystem. In fact, from Appendix 1—figure 3, we observe being outside of the domain of validity of Ssystem 2 (e.g., ${\overline{p}}_{0}<0.5$) yields to solutions that tend to the zero steady state.
We use DSA to select among all possible ODE models for normal hematopoiesis consistent with the lineage diagram shown in Figure 1A in the main text in which each parameter ${p}_{0}$ , ${p}_{1}$ , ${q}_{1}$ , ${\eta}_{1}$ , and ${\eta}_{2}$ is either unregulated (constant) or is subject to positive or negative regulation from most one cell type in the lineage. As described in the main text, there are 59,049 possible models, counting each combination as a model. We implement an automated implementation of DSA that enables an efficient exploration of this large space of models. Eliminating models with no valid Ssystems and those with unstable equilibria, we eliminate all but the four model classes shown in Figure 1B in the main text.
2 Mathematical model
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 TKIinduced death rates of the CML stem and MPP cells. With the introduction of differentiation therapy, Equations 30 and 38 instead become
3 Parameter estimation
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 gridsearch algorithm was used to find parameter values that demonstrate steadystate 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 ${\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 pseudoPython.
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.
4 Simulations of transplant experiment
As described in the main text, we simulated a transplant experiment in a transgenic mouse model of CML performed in Reynaud et al., 2011. In this experiment, either leukemic stem HSC ${\mathrm{}}^{\mathrm{L}}$ or leukemic MPP ${\mathrm{}}^{\mathrm{L}}$ cells were implanted into sublethally irradiated mice. Transplantation of HSC ${\mathrm{}}^{\mathrm{L}}$ enables engraftment and myeloid cell production that leads to CML. On the other hand, transplanting MPP ${\mathrm{}}^{\mathrm{L}}$ cells does not allow for longterm engraftment but results in a larger fraction of donorderived lymphoid cells after 35 d. We modeled this experiment by reducing the number of cells in equilibrium to mimic the effects of sublethal radiation (see ‘Methods’). Here, we present results of a range of possible reductions of HSC ${\mathrm{}}^{\mathrm{L}}$ and MPP ${\mathrm{}}^{\mathrm{L}}$ cells, and tracked the outcomes when 4000 HSC ${\mathrm{}}^{\mathrm{L}}$ or MPP ${\mathrm{}}^{\mathrm{L}}$ were introduced after the decrements from equilibrium. We then determine which of our parameter sets are consistent with the experimental outcomes found in Reynaud et al., 2011 using a simple majority of myeloid cells for HSC ${\mathrm{}}^{\mathrm{L}}$ transplant and a simple majority of lymphoid cells for MPP ${\mathrm{}}^{\mathrm{L}}$ transplant as consistency criteria. The results are summarized in Figure 5—figure supplements 1–2. The pairwise parameter distributions of the 478 remaining parameter sets are shown in Figure 5—figure supplement 3.
5 Effective parameters
Here, we present the effective proliferation rates and selfrenewal and branching factors—that is, the values of these parameters that takes the feedback regulation into account. That is, the effective stem cell proliferation rate ${\eta}_{1}={\eta}_{1,max}/\left(1+{\gamma}_{2}{x}_{S}\right)$ in normal hematopoiesis and ${\eta}_{1}={\eta}_{1,max}/\left(1+{\gamma}_{2}\left({x}_{S}+{x}_{{S}^{L}}\right)\right)$ when CML stem cells are present. The other effective parameters are defined analogously.
In Figure 3—figure supplement 2, the effective HSC and MPP proliferation rates (${\eta}_{1}$ , ${\eta}_{2}$), the effective HSC and MPP selfrenewal fractions (${p}_{0}$ , ${p}_{1}$) and branching fraction ${q}_{1}$ are shown. The corresponding effective parameters are shown in Figure 6—figure supplement 1 when CML cells are introduced when the normal hematopoietic model system is at steady state and in response to treatment by TKIs (starts at the time labeled $t=0$, as indicated by the vertical line) for the cases shown in Figure 6 in the main text. In Appendix 1—figure 6, we plot the effective proliferation, selfrenewal, and branching parameters for 50 parameter sets under treatment at early times.
6 Distributions of model parameters grouped by response to TKI treatment using synthetic data (478 parameter sets)
In Figure 7—figure supplement 1, we plot the distributions of all the unregulated proliferation, selfrenewal, and branching parameters, as well as the feedback gains, grouped by response to TKI therapy. In particular, the blue color indicates achievement of MR3 by 50 mo (termed as responders) while orange indicates that MR3 is not achieved within 50 mo (termed as nonresponders). The only parameters that clearly delineate the responders from nonresponders are ${p}_{0,max}$ and ${\gamma}_{1}$ , with responders occurring at the lower values and nonresponders at the higher values.
7 Comparisons of prognostic criteria for predicting response to TKI therapy
7.1 Performance of prognostic criteria using synthetic data (478 parameter sets)
Using data generated from our 478 parameter sets, we tested whether prognostic criteria could correctly identify patients who achieve MR3 within 18 mo after therapy starts. We tested the performance of our prognostic criterion (relative change of transcript levels) against two existing clinical prognostics: halving time (the time it takes for the BCRABL1 transcripts to reach onehalf of their pretreatment value) and early molecular response (EMR) in which the fraction of BCRABL1 transcripts are 10% or less after 3 mo of treatment. We also tested a prognostic criterion based on the ratio of transcript levels. These different prognostic criteria were calculated for both the first and second 3 mo after the start of therapy (e.g., 0–3 mo and 3–6 mo). We calculated the corresponding receiving operating characteristic (ROC) curves and found the optimal threshold by maximizing the difference between true and false positive rates. The results are presented in Appendix 1—figure 7 and reveal that generally the ratio and relative change prognostics offer similar performance, but that both demonstrate somewhat better performance compared to the traditional prognostics. In addition, all the prognostic criteria are more accurate when applied 3–6 mo after the start of therapy than when applied during the 0–3month period. We did not calculate the ROC curves for EMR but rather we only plotted the point that corresponds to 10% transcript levels at 3 mo (black circle) and 1% at 6 mo (open black circle). The EMR prognostic criterion has fewer false positives but also fewer true positives than the other prognostics.
7.2 Performance of prognostic criteria using patient data
The prognostic criteria were tested on anonymized patient data obtained from Dr. Van Etten’s clinical practice, again asking whether MR3 at 18 mo after treatment could be correctly predicted. We used data in which patients were kept on the same therapy for 6 mo either from the start of therapy or after a change of therapy. For patients that achieved MR3 within 3 mo, we did not include their data at 6 mo. In Appendix 1—figure 8, the first two columns correspond to results when the same therapy is applied for 6 mo after patient diagnoses. The last two columns correspond to patients who have had a change of therapy, but the new therapy is maintained for 6 mo. We do not use the EMR as a prognostic in the cases when therapy is changed. The figure demonstrates that the prognostics are more accurate for the 3–6month period, as predicted from the synthetic data. Although the numbers of patients are small, the relative ratio prognostic criterion is at least as accurate as, or more accurate than, the other criteria. In Appendix 1—figure 9, the predictions of the prognostic criteria are grouped by whether the patients achieve or do not achieve MR3 by 18 mo and by time period. In Appendix 1—figure 10, the prognostic criteria data are aggregated into 3month windows, where 0–3 mo contains both 0–3 mo after start of therapy and after a therapy change. The 3–6month data is similarly aggregated. It is clear that the predictions using the 3–6month data are more accurate than those using the 0–3month window.
8 Combination therapy parameter distributions by response
Finally, in Figure 8—figure supplement 1, we show how the distributions of model parameters from Figure 7—figure supplement 1 change when TKI therapy is combined with a differentiation promoter. The blue hatching indicates a nonresponder that becomes a responder while the yellow hatching indicates that a responder becomes a nonresponder. Here, response is defined as achieving MR3 in 50 mo. We observe that differentiation therapy is very effective in driving nonresponders at large ${p}_{0,max}$ and ${\gamma}_{1}$ to become responders, but also drives parameter sets that responded to TKI monotherapy at small values of ${p}_{0,max}$ and ${\gamma}_{1}$ to no longer achieve MR3 at 50 mo (nonresponders).
9 Impact of intrinsic differences between normal and leukemic cells
In Figure 4—figure supplements 2–4, we explored perturbations ranging from 90 to 110% of the original parameter values from Appendix 1—tables 4 and 5 to explore the effect of intrinsic differences between normal and leukemic cells. One exception is ${p}_{0,max}^{L}$ , which is limited to a perturbation range of 90–100% due to biological constraints as described in the main text. For the single parameter set for normal cells from Appendix 1—table 4, perturbations in most leukemic parameters yielded insignificant differences at the cellular and response dynamics levels. The three parameters that did see significant sensitivity to perturbation were the leukemic stem cellspecific parameters ${p}_{0,max}^{L}$ , $\frac{{\gamma}_{1}^{L}}{{\gamma}_{1}}$ , and $TK{I}_{HSC}$ . To determine whether this sensitivity applies to the entire population of parameter sets, we explored sweeps of the newly added leukemic parameters and their associated feedback gain shown in the heat maps of Appendix 1—figure 11 and Figure 4—figure supplement 6. For a parameter combination to be considered useful, there must be a region of the left plots that is of a lower value, such that in a majority of cases leukemic cells can dominate the system. Along the bottom right plot, the region should be neither fully light or fully dark to give regions where there are both responsive parameter sets and nonresponsive parameter sets. Through these parameter combination studies, we find that even on the broader parameter set population ${p}_{0,max}^{L}$ and $\frac{{\gamma}_{1}^{L}}{{\gamma}_{1}}$ are the only cases with significant differences in overall qualitative outcomes. Additionally, we find that the domains of ${p}_{0,max}^{L}$ and $\frac{{\gamma}_{1}^{L}}{{\gamma}_{1}}$ are relatively restricted with the only viable values of ${p}_{0,max}^{L}$ being roughly equivalent with ${p}_{0,max}$ with necessary decreases in $\frac{{\gamma}_{1}^{L}}{{\gamma}_{1}}$ for lower $\frac{{p}_{0,max}^{L}}{{p}_{0,max}}$ values to ensure parameter combination viability.
10 Effect of leukemic stem cell parameters on response to TKI therapy
By finding similar qualitative regions of the leukemic stem cell selfrenewal heat map in Appendix 1—figure 11 (lower right), we examined the distributions of response for a few combinations of parameters shown in Appendix 1—figure 12. We found that it is possible to achieve similar distributions with different combinations of parameters. When ${p}_{0,max}^{L}/{p}_{0,max}=0.8$ and ${\gamma}_{1}^{L}/{\gamma}_{1}=0.2$ , in 33% of the cases the treatment is unsuccessful in achieving MR3 in 50 mo. When ${p}_{0,max}^{L}/{p}_{0,max}=0.8$ and ${\gamma}_{1}^{L}/{\gamma}_{1}=0.1$ , in 20% of the cases the treatment is unsuccessful in achieving MR3 in 50 mo. When ${p}_{0,max}^{L}/{p}_{0,max}=0.9$ and ${\gamma}_{1}^{L}/{\gamma}_{1}=0.3$ , in 30% of the cases the treatment is unsuccessful in achieving MR3 in 50 mo. When ${p}_{0,max}^{L}/{p}_{0,max}=1.0$ and ${\gamma}_{1}^{L}/{\gamma}_{1}=0.5$, in 39% of the cases the treatment is unsuccessful in achieving MR3 in 50 mo.
When we analyzed the role of the pretreatment leukemic stem cell proportion on response to TKI therapy (Figure 7A), we find the results to agree qualitatively across the leukemic parameter combinations (Appendix 1—figure 13). From this, to determine whether the $p}_{0$ or $p}_{0}^{L$ is the true predictor of response to TKI therapy $p}_{0$, we first calculated a quantity we termed a characteristic effective selfrenewal fraction to attempt to group these combinations by similarity. The characteristic selfrenewal fractions for leukemic and normal stem cells are defined as ${\overline{p}}_{0}^{L}={p}_{0,max}^{L}/\left(1+{\gamma}_{1}^{L}\overline{N}\right)$ and ${\overline{p}}_{0}={p}_{0,max}/\left(1+{\gamma}_{1}\overline{N}\right)$ . We take $\overline{N}={10}^{5}$ to be a characteristic value of the size of the MPP population. We then analyzed the behavior as a function of the maximal HSC selfrenewal fraction $p}_{0,max$ and $\overline{p}}_{0,max}^{L}/{p}_{0,max$ . The results are shown in Appendix 1—figure 14 where all 478 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, using blue and yellow colors to denote responders and nonresponders. We observe that when the fitness of the $HS{C}^{L}$ (as measured by $\overline{p}}_{0,max}^{L}/{p}_{0,max$) is sufficiently low (e.g., ${\overline{p}}_{0,max}^{L}/{p}_{0,max}<0.65$), all the systems respond to TKI therapy. When the $HS{C}^{L}$ increase in fitness, the number of nonresponders increases but nonresponders are only observed when $p}_{0,max$ is above a critical threshold, which depends on $\gamma}_{1}^{L$ through $\overline{p}}_{0,max}^{L}/{p}_{0,max$ . Note that the fitness of the HSC^{L} can be increased either by increasing $p}_{0,max}^{L$ or $\gamma}_{1}^{L$ or both.
The overall fitness of the leukemic stem cells relative to that of the normal cells determines whether CML will develop and whether treatments will succeed or fail. This is shown in Figure 7—figure supplement 2. The relative fitness of the CML cells is measured by the ratios of characteristic values of the $HS{C}^{L}$ and $HSC$ selfrenewal fractions: $\overline{p}}_{0}^{L}/{\overline{p}}_{0$ . Here, all 478 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. The larger the relative fitness, the more likely that CML will develop and take
over the system (Figure 4D) and that the system will be refractory to TKI treatment (Figure 7—figure supplement 2).
As we described in the main text, we estimated the relative fitness of $HS{C}^{L}$ as ${\overline{p}}_{0}^{L}/{\overline{p}}_{0}\approx 0.70$ . In Appendix 1—figures 15–16 and Figure 7B, we plot the bivariate histogram for response to TKI treatment for combinations of leukemic stem cell parameters such that ${\overline{p}}_{0}^{L}/{\overline{p}}_{0}\ge 0.50$ (Appendix 1—figure 15), ${\overline{p}}_{0}^{L}/{\overline{p}}_{0}\ge 0.60$ (Appendix 1—figure 16), and ${\overline{p}}_{0}^{L}/{\overline{p}}_{0}\ge 0.7$ (Figure 7B). Through exploration of bivariate and marginal distributions, we find that $p}_{0,max$ is capable of separating response, while the fitness $\overline{p}}_{0}^{L}/{p}_{0,max$ does not have a clear delineation (Appendix 1—figures 15–16 and Figure 7B). Additionally, in our virtual patient population we consider variations of ${p}_{0,max}$ for individual biological variation with CML cells operating in a similar capacity across virtual patients to be more meaningful.
11 The effect of leukemic stem cell parameters on combination therapy and prognostic criterion
We checked the combinations of $\frac{{p}_{0,max}^{L}}{{p}_{0,max}}$ and $\frac{{\gamma}_{1}^{L}}{{\gamma}_{1}}$ to ensure that the effectiveness of combination therapy and the accuracy of our prognostic criterion are largely unchanged by varying the leukemic stem cell parameters. In Figure 8—figure supplement 2, we see combination therapy is still successful at improving the proportion of responders. We find that the optimal values from Figure 7E, where a single set of leukemic parameters was used, need to be modified when all the leukemic parameter combinations that have significant takeover proportions ($>60\mathrm{\%}$) are considered. Nevertheless, using synthetic data we find that the 3–6month time frame has a better predictivity of TKI response than does the 0–3month time frame (Appendix 1—figure 17).
Data availability
Modelling code and parameter set data are available in a Github repository. Patient data is unavailable publicly as it could be used to potentially identify the patients. Deidentified raw patient transcript data will be made available to qualified researchers (academic or industry) upon request to Dr. Van Etten at vanetten@hs.uci.edu.
References

Retention but significant reduction of BcrAbl transcript in hematopoietic stem cells in chronic myelogenous leukemia after imatinib therapyInternational Journal of Hematology 88:471–475.https://doi.org/10.1007/s1218500802211

MIP1α/CCL3mediated maintenance of leukemiainitiating cells in the initiation process of chronic myeloid leukemiaThe Journal of Experimental Medicine 210:2661–2673.https://doi.org/10.1084/jem.20130112

"Emergency" granulopoiesis in GCSFdeficient mice in response to Candida albicans infectionBlood 95:3725–3733.

Understanding normal and pathological hematopoietic stem cell biology using mathematical modellingCurrent Stem Cell Reports 7:109–120.https://doi.org/10.1007/s40778021001919

Macrophage inflammatory protein1 alpha receptors are present on cells enriched for CD34 expression from patients with chronic myeloid leukemiaBlood 86:4270–4277.

Natural course and biology of CMLAnnals of Hematology 94 Suppl 2:S107–S121.https://doi.org/10.1007/s002770152325z

Human chronic myeloid leukemia stem cells are insensitive to imatinib despite inhibition of BcrAbl activityThe Journal of Clinical Investigation 121:396–409.https://doi.org/10.1172/JCI35721

Feedback mechanisms control coexistence in a stem cell model of acute myeloid leukaemiaJournal of Theoretical Biology 401:43–53.https://doi.org/10.1016/j.jtbi.2016.04.002

Evolutionary dynamics of chronic myeloid leukemiaGenes & Cancer 1:309–315.https://doi.org/10.1177/1947601910371122

Predicting time to relapse in acute myeloid leukemia through stochastic modeling of minimal residual disease based on clonality dataComputational and Systems Oncology 1:e1026.https://doi.org/10.1002/cso2.1026

Retinoic acid signaling in myelopoiesisCurrent Opinion in Hematology 15:37–41.https://doi.org/10.1097/MOH.0b013e3282f20a9c

LongTerm followup of the French stop imatinib (STIM1) study in patients with chronic myeloid leukemiaJournal of Clinical Oncology 35:298–305.https://doi.org/10.1200/JCO.2016.68.2914

Multicenter independent assessment of outcomes in chronic myeloid leukemia patients treated with imatinibJournal of the National Cancer Institute 103:553–561.https://doi.org/10.1093/jnci/djr060

LongTerm outcomes of imatinib treatment for chronic myeloid leukemiaThe New England Journal of Medicine 376:917–927.https://doi.org/10.1056/NEJMoa1609324

Integration of mathematical model predictions into routine workflows to support clinical decision making in haematologyBMC Medical Informatics and Decision Making 20:28.https://doi.org/10.1186/s129110201039x

Mathematical modeling of genesis and treatment of chronic myeloid leukemiaCells Tissues Organs 188:236–247.https://doi.org/10.1159/000118786

GranulocyteMacrophage progenitors as candidate leukemic stem cells in blastcrisis CMLThe New England Journal of Medicine 351:657–667.https://doi.org/10.1056/NEJMoa040258

SoftwareCml_Treatment, version swh:1:rev:91df5a6dd841d97a2dbca6bc42a976749cfe52e3Software Heritage.

Differential regulation of myeloid leukemias by the bone marrow microenvironmentNature Medicine 19:1513–1517.https://doi.org/10.1038/nm.3364

Improving cancer treatments via dynamical biophysical modelsPhysics of Life Reviews 39:1–48.https://doi.org/10.1016/j.plrev.2021.10.001

Establishment of a reproducible model of chronicphase chronic myeloid leukemia in NOD/SCID mice using bloodderived mononuclear or CD34+ cellsBlood 91:630–640.

Optimal experimental design for mathematical models of haematopoiesisJournal of the Royal Society, Interface 18:20200729.https://doi.org/10.1098/rsif.2020.0729

Dynamical modelling of haematopoiesis: an integrated view over the system in homeostasis and under perturbationJournal of The Royal Society Interface 10:20120817.https://doi.org/10.1098/rsif.2012.0817

Regulation of hematopoiesis and hematological disease by TGFβ family signaling moleculesCold Spring Harbor Perspectives in Biology 9:a027987.https://doi.org/10.1101/cshperspect.a027987

Critical differences in hematopoiesis and lymphoid development between humans and miceJournal of Clinical Immunology 33:711–715.https://doi.org/10.1007/s1087501298443

Mathematical modelling of the hematopoietic stem cellniche system: clonal dominance based on stem cell fitnessJournal of Theoretical Biology 518:110620.https://doi.org/10.1016/j.jtbi.2021.110620

Blood cell Dynamics: half of a century of modelingMathematical Modelling of Natural Phenomena 10:182–205.https://doi.org/10.1051/mmnp/201611106

Strategies for calibrating models of biologyBriefings in Bioinformatics 21:24–35.https://doi.org/10.1093/bib/bby092

HematopoiesisCold Spring Harbor Perspectives in Biology 4:a008250.https://doi.org/10.1101/cshperspect.a008250

Designing combination therapies using multiple optimal controlsJournal of Theoretical Biology 497:110277.https://doi.org/10.1016/j.jtbi.2020.110277

Normal and leukemic SCIDrepopulating cells (Src) coexist in the bone marrow and peripheral blood from CML patients in chronic phase, whereas leukemic Src are detected in blast crisisBlood 87:1539–1548.

Lymphohematopoietic engraftment in minimally myeloablated hostsBlood 91:3681–3687.

Clonal selection and therapy resistance in acute leukaemias: mathematical modelling explains different proliferation patterns at diagnosis and relapseJournal of the Royal Society, Interface 11:20140079.https://doi.org/10.1098/rsif.2014.0079

High level engraftment of NOD/SCID mice by primitive normal and leukemic hematopoietic cells from patients with chronic myeloid leukemia in chronic phaseBlood 91:2406–2414.

ConferenceGlobal Stability Analysis of Strictly Positive Steady State for a Surviving Hematopoietic Stem Cells Models2018 Annual American Control Conference (ACC.https://doi.org/10.23919/ACC.2018.8430856
Article and author information
Author details
Funding
National Institutes of Health (1U54CA21737801A1)
 Arthur D Lander
 Richard A Van Etten
 John Lowengrub
National Institutes of Health (P30CA062203)
 Arthur D Lander
National Institutes of Health (R01 CA090576)
 Richard A Van Etten
National Science Foundation (DMS1763272)
 Arthur D Lander
 John Lowengrub
National Science Foundation (DMS1936833)
 John Lowengrub
National Science Foundation (DMS1953410)
 John Lowengrub
National Science Foundation (GRFP 16588)
 Abdon Iniguez
Simons Foundation (594598QN)
 Arthur D Lander
 John Lowengrub
National Institute of General Medical Sciences (GM136624)
 Jonathan Rodriguez
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
The authors thank the reviewers for their careful analysis of the manuscript and their suggestions, which significantly improved it. ADL, JSL, and RAV acknowledge the National Institutes of Health for partial support through grant nos. 1U54CA21737801A1 for a National Center in Cancer Systems Biology at UC Irvine and P30CA062203 for the Chao Family Comprehensive Cancer Center at UC Irvine. In addition, ADL and JSL acknowledge support from DMS1763272 and the Simons Foundation (594598QN) for an NSFSimons Center for Multiscale Cell Fate Research, RAV acknowledges support from R01 CA090576. JSL also acknowledges NSF grants DMS1936833 and DMS1953410. JR and AI each acknowledge support from NSF graduate research fellowships.
Ethics
This study was performed in strict accordance with the recommendations in the Guide for the Care and Use of Laboratory Animals of the National Institutes of Health. All of the animals were handled according to approved institutional animal care and use committee (IACUC) protocol (AUP19159) of the University of California, Irvine.
Version history
 Received: October 12, 2022
 Preprint posted: October 14, 2022 (view preprint)
 Accepted: April 26, 2023
 Accepted Manuscript published: April 28, 2023 (version 1)
 Accepted Manuscript updated: May 3, 2023 (version 2)
 Accepted Manuscript updated: May 3, 2023 (version 3)
 Accepted Manuscript updated: May 4, 2023 (version 4)
 Version of Record published: May 25, 2023 (version 5)
Copyright
© 2023, Rodriguez, Iniguez et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 417
 Page views

 90
 Downloads

 1
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Cancer Biology
 Structural Biology and Molecular Biophysics
Abelson tyrosine kinase (Abl) is regulated by the arrangement of its regulatory core, consisting sequentially of the SH3, SH2, and kinase (KD) domains, where an assembled or disassembled core corresponds to low or high kinase activity, respectively. It was recently established that binding of type II ATP site inhibitors, such as imatinib, generates a force from the KD Nlobe onto the SH3 domain and in consequence disassembles the core. Here, we demonstrate that the Cterminal αIhelix exerts an additional force toward the SH2 domain, which correlates both with kinase activity and type II inhibitorinduced disassembly. The αIhelix mutation E528K, which is responsible for the ABL1 malformation syndrome, strongly activates Abl by breaking a salt bridge with the KD Clobe and thereby increasing the force onto the SH2 domain. In contrast, the allosteric inhibitor asciminib strongly reduces Abl’s activity by fixating the αIhelix and reducing the force onto the SH2 domain. These observations are explained by a simple mechanical model of Abl activation involving forces from the KD Nlobe and the αIhelix onto the KD/SH2SH3 interface.

 Biochemistry and Chemical Biology
 Cancer Biology
Poly(ADPribose)ylation or PARylation by PAR polymerase 1 (PARP1) and dePARylation by poly(ADPribose) glycohydrolase (PARG) are equally important for the dynamic regulation of DNA damage response. PARG, the most active dePARylation enzyme, is recruited to sites of DNA damage via pADPrdependent and PCNAdependent mechanisms. Targeting dePARylation is considered an alternative strategy to overcome PARP inhibitor resistance. However, precisely how dePARylation functions in normal unperturbed cells remains elusive. To address this challenge, we conducted multiple CRISPR screens and revealed that dePARylation of S phase pADPr by PARG is essential for cell viability. Loss of dePARylation activity initially induced Sphasespecific pADPr signaling, which resulted from unligated Okazaki fragments and eventually led to uncontrolled pADPr accumulation and PARP1/2dependent cytotoxicity. Moreover, we demonstrated that proteins involved in Okazaki fragment ligation and/or base excision repair regulate pADPr signaling and cell death induced by PARG inhibition. In addition, we determined that PARG expression is critical for cellular sensitivity to PARG inhibition. Additionally, we revealed that PARG is essential for cell survival by suppressing pADPr. Collectively, our data not only identify an essential role for PARG in normal proliferating cells but also provide a potential biomarker for the further development of PARG inhibitors in cancer therapy.