Author response: Predictive nonlinear modeling of malignant myelopoiesis and tyrosine kinase inhibitor therapy
Year: 2023
Type: peer-review
Abstract: Full text Figures and data Side by side Abstract Editor's evaluation Introduction Results Discussion Methods Appendix 1 Data availability References Decision letter Author response Article and author information Metrics 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 BCR-ABL1 tyrosine kinase. Tyrosine kinase inhibitors (TKIs) have proved effective in treating CML, but there is still a cohort of patients who do not respond to TKI therapy even in the absence of mutations in the BCR-ABL1 kinase domain that mediate drug resistance. To discover novel strategies to improve TKI therapy in CML, we developed a nonlinear mathematical model of CML hematopoiesis that incorporates feedback control and lineage branching. Cell–cell interactions were constrained using an automated model selection method together with previous observations and new in vivo data from a chimeric BCR-ABL1 transgenic mouse model of CML. The resulting quantitative model captures the dynamics of normal and CML cells at various stages of the disease and exhibits variable responses to TKI treatment, consistent with those of CML patients. The model predicts that an increase in the proportion of CML stem cells in the bone marrow would decrease the tendency of the disease to respond to TKI therapy, in concordance with clinical data and confirmed experimentally in mice. The model further suggests that, under our assumed similarities between normal and leukemic cells, a key predictor of refractory response to TKI treatment is an increased maximum probability of self-renewal of normal hematopoietic stem cells. We use these insights to develop a clinical prognostic criterion to predict the efficacy of TKI treatment and design strategies to improve treatment response. The model predicts that stimulating the differentiation of leukemic stem cells while applying TKI therapy can significantly improve treatment outcomes. 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.sa0 Decision letter Reviews on Sciety eLife's review process Introduction Chronic myeloid leukemia (CML) is a myeloproliferative neoplasm of the hematopoietic system, which normally produces billions of mature myeloid and erythroid cells on a daily basis, is tightly regulated, and accommodates massive increases in the production of individual cell types in response to physiological and pathological stresses. The hematopoietic system is organized hierarchically as a collection of progressively more differentiated cells starting from a hematopoietic stem cell (HSC) located in the bone marrow (BM) and ending with postmitotic terminally differentiated myeloid and lymphoid cells (Rieger and Schroeder, 2012; Liggett and Sankaran, 2020). CML is characterized by an overproduction of myeloid cells including mature granulocytes (neutrophils, basophils, and eosinophils) and their immediate precursors (metamyelocytes, myelocytes, and promyelocytes), and of myeloid progenitors (Jamieson et al., 2004) including multipotential progenitors (MPPs) and committed progenitors (common myeloid progenitors [CMP], granulocyte-macrophage progenitors [GMPs], and megakaryocyte-erythroid progenitors [MEPs]). Untreated, the disease has three distinct phases (Chereda and Melo, 2015). In the initial 'chronic' phase, the differentiation of myeloid progenitors is essentially normal, resulting in excessive levels of mature postmitotic neutrophils and their immediate precursors. In later stages of the disease (accelerated phase and blast crisis), differentiation is reduced and expansion of immature progenitors is observed. Additional clonal karyotypic abnormalities are typically only observed during the accelerated and blast crisis phases (Hehlmann et al., 2020). CML has one of the simplest cancer genomes. It is driven by a single genetic abnormality arising somatically in an HSC, the Philadelphia (Ph) chromosome, the result of a balanced translocation between chromosomes 9 and 22 that creates a fusion of the genes for BCR and ABL1. The product of the BCR-ABL1 fusion gene is a dysregulated cytoplasmic protein-tyrosine kinase, BCR-ABL1. CML thus represents a natural model of dysregulated granulocytopoiesis (Quintás-Cardama and Cortes, 2009). Cell biological studies have shown that Ph+ cells expressing markers of normal HSC are capable of engrafting immunodeficient mice (Sirard et al., 1996; Lewis et al., 1998), implying that these cells are leukemia-initiating or leukemic 'stem' cells (LSCs). More mature committed progenitors in CML, like normal progenitors, lack sustained self-renewal capacity and cannot stably engraft immunodeficient mice nor generate hematopoietic colonies in vitro upon serial replating (Huntly et al., 2004). The proportion of LSCs in the BM is highly variable across CML patients at diagnosis and can range from a few percent to nearly 100% (Petzer et al., 1996; Diaz-Blanco et al., 2007; Abe et al., 2008; Thielen et al., 2016), perhaps reflecting different periods of time patients spend in chronic phase before they are diagnosed, different rates of disease progression, or both. There is persuasive experimental evidence of significant feedback regulation of different cell compartments in the dynamics of myeloid cell production in both normal and CML hematopoiesis, including signaling between the normal and CML cells (Jiang et al., 1999; Devireddy et al., 2005; Vicente-Dueñas et al., 2009; Naka et al., 2010; Reynaud et al., 2011; Zhang et al., 2012; Krause et al., 2013; Walenda et al., 2014; Welner et al., 2015). For instance, experiments in a mouse model of CML provided evidence that IL-6, produced by leukemic neutrophils, blocked MPP differentiation toward a lymphoid fate, implying feedback from the myeloid lineage onto MPPs (Reynaud et al., 2011). Surprisingly, our knowledge of the details of feedback regulation in hematopoiesis is still incomplete, especially for granulopoiesis, where even late-stage feedback interactions are poorly understood. For example, two cytokines, granulocyte colony-stimulating factor (G-CSF) and granulocyte-macrophage colony-stimulating factor (GM-CSF), can pharmacologically increase neutrophil production, but mice lacking both cytokines maintain baseline neutrophil levels and can still increase neutrophil production in response to infection (Basu et al., 2000). In many cases, it is not known which cell types are providing and receiving the feedback, what signals are used, and what aspects of proliferative cell behavior they influence (i.e., proliferation rates, renewal probability, or progeny fate choice). In spite of this knowledge deficit, CML can be treated quite effectively using selective small-molecule tyrosine kinase inhibitors (TKIs) of the BCR-ABL1 kinase. TKIs such as imatinib, dasatinib, and nilotinib, which inhibit proliferation and increase apoptosis of Ph+ cells, have dramatically lowered CML death rates (Gambacorti-Passerini et al., 2011). The response to TKI therapy in CML is monitored primarily by determining the level of BCR-ABL1 mRNA transcripts in peripheral blood, normalized to a control RNA and expressed as a percentage on an International Scale (Arora and Press, 2017). BCR-ABL1 transcript levels, an approximation of the proportion of circulating malignant cells at any given time, generally decrease exponentially in patients responding to TKI therapy resulting in at least two distinct slopes when plotted semi-logarithmically—an initial rapid decline attributed to TKI-induced killing of more mature myeloid cells, and a subsequent slower decline postulated to represent lower death rates of more primitive leukemic stem/progenitor cells (Michor et al., 2005). Clinical resistance to TKI therapy in CML is a significant problem and is classified as acquired resistance (increasing BCR-ABL1 transcript levels following a substantial decrease) or primary resistance (lack of an adequate initial response). Many patients with acquired resistance have developed mutations in the BCR-ABL1 kinase domain that mediate pharmacological resistance to the TKI (Ernst and Hochhaus, 2012). By contrast, 10–15% of newly diagnosed CML patients fail to achieve an 'early molecular response,' defined as the level of BCR-ABL1 transcripts being less than 10% at 3 mo (Hanfstein et al., 2012; Marin et al., 2012). Clinical data indicate that switching TKIs may not benefit these patients (Yeung et al., 2012; Yeung et al., 2015), suggesting that this group is destined to do poorly regardless of the specific inhibitor used. BCR-ABL1 mutations are generally not present in this group of patients (Zhang et al., 2009; Pietarinen et al., 2016), and thus the mechanism(s) underlying this primary resistance is unclear. We hypothesized that these variable patient responses to TKI therapy arise from nonlinearity introduced by non-cell-autonomous interactions between normal and CML cells. To test this hypothesis, we developed a novel mathematical model of CML hematopoiesis and TKI treatment that incorporates lineage branching and interactions between normal and CML cells through feedback and feedforward regulation. Mathematical modeling of leukemia has a long history aimed at understanding disease progression and improving treatment response using single and combination targeted therapies and immunotherapy (Whichard et al., 2010; Pujo-Menjouet, 2015; Brunetti et al., 2021; Kuznetsov et al., 2021; Roeder and Glauche, 2021). Further, recent efforts have been made to integrate mathematical modeling in clinical decision-making to design personalized therapies (Hoffmann et al., 2020; Engelhardt and Michor, 2021). Many models of leukemia have utilized simplified lineage architectures and minimal feedback (Roeder et al., 2006; Komarova and Wodarz, 2007; Horn et al., 2008; Foo et al., 2009; Hähnel et al., 2020; Pedersen et al., 2021). While these models can be made to fit the multiphasic disease response data of CML to TKI treatment, the simplicity of the models can make these fitted parameters of limited clinical value. More physiologically accurate, nonlinear models that account for cell–cell signaling and lineage branching are expected to improve clinical relevance. Mathematical models that incorporate feedback signaling have been developed in normal (Engel et al., 2004; Marciniak-Czochra et al., 2009; Mahadik et al., 2019, Mon Père et al., 2021) and diseased (Wodarz, 2008; Sachs et al., 2011; Krinner et al., 2013; Stiehl et al., 2014; Stiehl et al., 2015; Crowell et al., 2016; Woywod et al., 2017; Jiao et al., 2018; Stiehl et al., 2018; Zenati et al., 2018; Park et al., 2019; Sharp et al., 2020) hematopoiesis. Because of the vast number of possible ways in which feedback models of normal hematopoiesis and leukemia can be configured, mathematical models tend to greatly simplify the lineage architectures and the feedback interactions among the cell types. For example, Manesso et al., 2013 developed a hierarchical ordinary differential equation (ODE) model of normal hematopoiesis containing multiple cell types and branch points (16 cell types and 4 branch points) in the lineage tree. Limiting the feedback loops to involve only local, negative regulation (e.g., regulation by self and immediate progenitor/progeny in the lineage tree) results in about 106 models, which enabled the use of a stochastic optimization algorithm to obtain parameters consistent with homeostasis and a requirement for a rapid return to equilibrium following system perturbations. In the context of leukemia, the model architectures are typically much simpler. Generally, models of leukemia introduce a parallel mutant lineage with the same structure as that used to model the normal hematopoietic cells but with different parameters. For example, Wodarz developed an unbranched lineage ODE model of normal and leukemia stem and differentiated cells in which feedback from the differentiated cells controlled whether the stem cells divided symmetrically or asymmetrically, and demonstrated this provides a mechanism for blast crisis in CML to occur without additional mutations (Wodarz, 2008). Krinner et al. incorporated positive and negative feedback regulation of differentiation and proliferation in an unbranched lineage model that combined a discrete agent-based model for the stem cell compartment with an ODE system for the progenitor and differentiated cells to provide a detailed view of the stem cell dynamics and to test the effect of therapies (Krinner et al., 2013). Stiehl et al., 2015 developed an unbranched lineage ODE model of normal and leukemic cells in which only negative feedback regulation of stem and progenitor cell self-renewal fractions was considered, and this was further limited to arise only from factors produced by the postmitotic, mature normal and leukemic cells. Later work extended this approach to investigate clonal selection and therapy resistance (Stiehl et al., 2014), the role of cytokines on leukemia progression (Stiehl et al., 2018), combination treatment strategies (Banck and Görlich, 2019), and niche competition (Stiehl et al., 2020). Clonal competition was also considered in an ODE feedback model of CML (Woywod et al., 2017) and a stochastic model with feedback (Dinh et al., 2021). Simpler unbranched lineage models of normal and leukemic cells in which only the normal cells respond to feedback but normal and leukemic cells compete for space in the BM have been used to investigate regimes of coexistence of normal and leukemic cells (Crowell et al., 2016; Jiao et al., 2018) and design combination therapies using optimal control algorithms (Sharp et al., 2020). Here, we develop a nonlinear ODE model of normal and CML hematopoiesis using a general approach that integrates an automated method, design space analysis (DSA; Fasani and Savageau, 2010), with data gleaned from previously published experiments, and from two new in vivo experiments presented here that separately decrement the number of stem cells and terminally differentiated myeloid cells in the BM of mice. This approach enables us to systematically select among plausible model architectures and signaling interactions without a priori knowledge of which cells are providing and receiving signaling stimuli. We start with a model for normal hematopoiesis that accounts for stem, multipotent progenitor cells, and two types of terminally differentiated cells representing the myeloid and lymphoid lineage branches. This approach allows us to reduce the potential model space from about 60,000 models to a single model class and reveals the existence of feedforward and feedback mechanisms. We then extend the model to incorporate CML hematopoiesis by introducing a parallel lineage of CML cells with the same model architecture but with different parameters. The model captures the dynamics of CML at various stages of the disease and exhibits variable response to TKI treatment consistent with that observed in clinical data. The model suggests biomarkers of primary resistance, identifies the underlying mechanisms governing the response to TKI therapy, and suggests new treatment strategies. 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 (TDm) and lymphoid (TDl). The HSC self-renew with fraction (e.g., probability) p0 or differentiate with fraction 1-p0. That is, the fraction of HSC that remain as HSC after division is p0. The MPPs self-renew with fraction p1 and differentiate into either lymphoid or myeloid cells with fractions q1 and 1-p1-q1, respectively. The HSC and MPPs divide with rates η1 and η2 and the myeloid and lymphoid cells die at rates dm and dl, respectively. The ODEs that govern the dynamics of the cells are given in 'Methods.' Assuming that there is either positive or negative regulation of the self-renewal and differentiation probabilities and division rates of any cell type from any other cell type results in 59,049 models, counting each combination of regulated signaling as a separate model. Figure 1 Download asset Open asset Branched lineage model of normal hematopoiesis with feedback regulation. (A) Branched lineage model consisting of hematopoietic stem cells (HSC; S), multipotent progenitor cells (MPP; P), and postmitotic, terminally differentiated myeloid (TDm) and lymphoid (TDl) cells. Modulation of the HSC and MPP self-renewal fractions (p0 and p1), division rates (η1 and η2), and fate switching probability (q1) through feedback can arise from any cell type. The different line styles correspond to regulation by a particular cell type (dashed for S, solid for P, dot-dashed for TDl, and dotted for TDm). (B) Using Design Space Analysis, four candidate model classes are identified that differ in how HSCs are regulated. (C) Using biological data from the literature as discussed in the text, we reduced the model space by hypothesizing that factors secreted by terminally differentiated myeloid cells direct the fate of MPPs (e.g., IL-6) and those by MPPs suppress HSC self-renewal (e.g., CCL3). To select the most physiologically accurate models, we first filtered the models using an automated approach (DSA) developed by Savageau and co-workers (Savageau et al., 2009; Fasani and Savageau, 2010; Lomnitz and Savageau, 2016) that enables models to be distinguished based on their range of qualitatively distinct behaviors, without relying on knowledge of specific values of the parameters. This method relies on identifying boundaries in parameter space that separate qualitative behaviors of a particular model, which is much more efficient than searching for model behaviors directly. The boundaries can be approximated from a sequence of inequalities that identify regions where one term on the right-hand side of each ODE (e.g., the rate of change) dominates all others in the sources (positive terms) and another dominates the sinks (negative terms). This is known as a dominant subsystem (S-system) of the model. The number of S-systems in each model depends on the number of combinations of positive and negative terms in the rates of change. If the equilibria of the S-systems, which are determined analytically, are not self-consistent (e.g., consistent with the assumed dominance of terms reflected in the inequalities) or the equilibria are not stable, then the S-system is rejected. If all the S-systems of a particular model are rejected, then that model is removed from further consideration. Models with at least one self-consistent and stable S-system are viable candidates for further analysis. DSA can be easily automated to make the analysis of very large numbers of equations feasible. Details are provided in 'Methods' and 'Appendix 1' (Section 1). The result of this procedure is the elimination of all but the four model classes shown in Figure 1B, which negative regulation of the stem cell self-renewal fraction but differ by where this regulation The models the classes at least one S-system and have qualitative The between models in a class in whether or not there is or regulation on the of the parameters from any of the cell types. This the number of possible models to work has feedback mechanisms in both normal and malignant hematopoiesis. is produced by differentiated myeloid cells and to MPPs toward a myeloid fate (Reynaud et al., 2011; Welner et al., 2015). negative feedback known as fate have been shown to provide an effective for control of cell proliferation and of in et al., 2015). The known as produced in BM by myeloid progenitors et al., 2016), to inhibit the proliferation and self-renewal of normal HSC et al., et al., 2018), but CML HSC are to et al., et al., 2013). In hypothesizing these we at a single model class as shown in Figure In this there are model which differ only in how the HSC and MPP cell division rates and the MPP self-renewal fraction are regulated. results that is a candidate feedback factor expressed in the myeloid compartment (TDm) with the to the fraction q1 of MPPs that differentiate into lymphoid cells. is a candidate factor negative feedback from the MPP onto HSC To further the models, we cell biological experiments in mice to information about cell–cell interactions by separately the stem cell and myeloid cell of HSC increases HSC and MPP proliferation in mice were treated with previously shown to be to HSC in the BM et al., The BM stem/progenitor compartment was by in and on and using the in Figure time points and the number of mice at each time were by a hierarchical for optimal experimental design of mathematical models of hematopoiesis et al., 2021). In the suggests time points after was with time points the time points provide more information about division rates, while the time points provide more information about the feedback parameters. after treatment, we observed an decrease in the of the HSC compartment in the mice by increase in proliferation rates for both HSC and MPPs There was significant in MPP and the system to equilibrium by results that the HSC negative feedback on their division rate and the division of MPPs through a negative feedforward on Figure Download asset Open asset cell sorting analysis of mouse bone marrow stem/progenitor cells and the lineage hematopoiesis model. (A) for hematopoietic stem cells defined as and multipotential progenitors defined as and in their (B) of HSC MPP and other compartments on and in the bone marrow (BM) of control mice and mice that (C) of HSC and MPP proliferation in and mice by on and Data are shown as the of myeloid cells as by in mice after of control or the of cells in BM of mice that were treated with control 3 or HSC MPP and other cell type from mice that or Data are shown as model of hematopoiesis with feedback with The negative feedback loops shown in correspond to those by previous experimental data (Reynaud et al., 2011; et al., 2018), while the negative feedback and feedforward loops in are by our cell experiments in of mature myeloid cells increases the MPP mice were treated with the and their BM was 1 after treatment This treatment in an decrease in mature BM myeloid cells, as by and and was by a increase in the of the MPP compartment and a decrease in the HSC compartment results that there is a negative feedback from the myeloid cells onto the MPP self-renewal fraction all these results into we at the model shown in Figure The negative feedback loops shown in correspond to those by previous experimental data, while the negative feedback and feedforward loops in are by the cell experiments presented 'Methods' for a detailed of the these data were from we that cell–cell signaling occurs in for model of hematopoiesis To parameters for the model in Figure a algorithm was The ODE model is given in 'Methods' and Appendix 1 (Section The model parameters and death rates, self-renewal and branching were using a for each Appendix 1 (Section and Appendix and 3 for details and a parameter parameter values were the model was for long If a parameter in values consistent with the range of values previously for a hematopoiesis model et al., that parameter was of possible parameter a of parameter were We further the candidate parameter by only those with large feedforward on the MPP division rate in to on the novel feedforward This reduced the number of parameter to and their are shown in Appendix of these parameter can be of as representing the of a patient by different individual for example, to or that result in a hematopoietic The different parameter thus model a range of across individual CML patients. The values of the parameters used are given in Appendix of hematopoiesis model DSA can be used to qualitative model behaviors and how the model is to of key parameters. Here, we on the feedback and on the HSC and MPP self-renewal Appendix for and for for other Figure 1). in Figure DSA identifies four regions in the and which the dynamics are by different Using a parameter in each design space by as a we a parameter in which we and in a range the of their In Figure the of each of the cell is starting from an initial in which there are only a number of The different correspond to the parameter in the four regions of the design space the dynamics are shown for the ODE The solid results from the parameter and the the range of behaviors when the parameters are The and correspond to the HSC and while the and correspond to the terminally differentiated myeloid and lymphoid cells. While the system to equilibrium for all parameter the approach to equilibrium is The dynamics in regions and are while those in regions and are not (e.g., the equilibria in regions and are stable while those in regions and are stable Further, the the the the approach toward The cell numbers and in each design space are different as In regions and the
Source: eLife
Institution: University of California, Irvine
Language: English
Subfield: Hematology
Field: Medicine
Domain: Health Sciences
Open Access status: gold
Estimated APC paid: $2,500
OpenAlex ID: W4378215295