Research Article

Predictive nonlinear modeling of malignant myelopoiesis and tyrosine kinase inhibitor therapy

Graduate Program in Mathematical, Computational and Systems Biology, University of California, Irvine, United States
Center for Complex Biological Systems, University of California, Irvine, United States
Department of Medicine, University of California, Irvine, United States
Department of Developmental and Cell Biology, University of California, Irvine, United States
Chao Family Comprehensive Cancer Center, University of California, Irvine, United States
Department of Biomedical Engineering, University of California, Irvine, United States
Department of Mathematics, University of California, Irvine, United States

Apr 28, 2023

https://doi.org/10.7554/eLife.84149

Open access
Copyright information

Version of Record

Accepted for publication after peer review and revision.

Download
Cite
Share
CommentOpen annotations (there are currently 0 annotations on this page).

Version of Record published: May 25, 2023 (This version)
Accepted Manuscript updated: May 4, 2023 (Go to version)
Accepted Manuscript updated: May 3, 2023 (Go to version)
Accepted Manuscript updated: May 3, 2023 (Go to version)
Accepted Manuscript published: April 28, 2023 (Go to version)
Accepted: April 26, 2023
Preprint posted: October 14, 2022 (Go to version)
Received: October 12, 2022

Figures
Tables
Additional files

25 figures, 6 tables and 1 additional file

Figures

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 (TD_m) and lymphoid (TD_l) cells. Modulation of the HSC and MPP self-renewal fractions (p₀ and p₁), division rates (η₁ and η₂), and fate switching probability (q₁) 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 TD_l, and dotted for TD_m). (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).

Figure 2

Download asset Open asset

Fluorescence-activated cell sorting (FACS) analysis of mouse Lin–Sca-1+c-Kit+ (LSK) bone marrow stem/progenitor cells and the proposed branched lineage hematopoiesis model.

(A) Gating schema for phenotyping hematopoietic stem cells (HSC, defined as LSK CD34–CD48–) and multipotential progenitors (MPP, defined as LSK CD34+ CD48+), and BrdU incorporation in their respective compartments. (B) Distributions of HSC (blue), MPP (orange), and other (purple) compartments on days 1, 3, and 7 in the bone marrow (BM) of control (CTRL) B6 mice and mice that received 50 cGy radiation. (C) Frequency of HSC and MPP proliferation in CTRL (gray bars) and irradiated (blue bars) mice measured by BrdU incorporation on days 1, 3, and 7. Data are shown as mean ± SEM. *p<0.05. (D) Representative histograms depicting the frequency of myeloid cells as measured by CD11b expression in mice 24 hr after intravenous administration of isotype control (iso) or RB68C5 (50 μg) antibody. (E) Left panel: bar graph showing the frequency of CD11b+ cells in BM of mice that were treated with isotype control antibody (Iso; orange bar, n = 3 ) or RB68C5 antibody (50 μg; blue bars, n = 3). Right panel: HSC (blue), MPP (orange), and other cell type (purple) frequencies from mice that received isotype or RB6-8C5 antibody. Data are shown as mean ± SEM. *p<0.05. (F) Proposed feedforward-feedback model of hematopoiesis with associated feedback strengths denoted with γ₁–γ₅. The negative feedback loops shown in blue correspond to those suggested by previous experimental data (Reynaud et al., 2011; Staversky et al., 2018), while the negative feedback and feedforward loops in red are supported by our cell depletion experiments in (**A–E**).

Figure 3 with 2 supplements

Download asset Open asset

Qualitative behavior of feedforward-feedback model and parameter sensitivity.

(A) The colored regions (i–iv) represent areas of design space in which there are distinct qualitative behaviors as a function of the feedback gains γ₁ and γ₃ of the hematopoietic stem cell (HSC) and multipotential progenitor (MPP) self-renewal fractions, respectively. White dots denote specific parameter combinations. (B) The dynamics for each cell compartment within each of the four design space regions (i–iv). Solid lines represent ordinary differential equation (ODE) solutions using the specific parameter combinations (black dots in A) while the lightly colored regions represent the range of ODE solutions resulting from perturbations in γ₁ and γ₃ in a range within 0.9–1.1 times their original values. The blue and black curves correspond to the HSC and MPPs, respectively, the green and turquoise curves correspond to the myeloid and lymphoid cells. (C) The return to equilibrium following partial depletion of mature myeloid cells (10%, 50%, 90%) using the parameter combination (white dot) in region iii.

Figure 3—figure supplement 1

Download asset Open asset

The dynamics from the steady state with perturbations of each parameter, as labeled, ranging from 90% to 110% of the original parameter value from Appendix 1—table 4.

The solid line represents the median of the perturbations, while the shaded region represents 95% of the perturbation dynamics.

Figure 3—figure supplement 2

Download asset Open asset

Effective parameters for proliferation, self-renewal, and branching as the solutions to the model for the normal hematopoietic system approach steady state.

Figure 4 with 6 supplements

Download asset Open asset

Extension of the model of hematopoiesis to chronic myeloid leukemia (CML).

(A) Schematic of two branched lineages consisting of normal and CML cell compartments. The two lineages share the same feedback architecture. The difference between the two lineages is the leukemic hematopoietic stem cell (HSC) self-renewal is less affected by negative feedback, denoted by p₀^L (see text). (B) Dynamics of hematopoiesis upon introduction of CML cells. We begin with having normal hematopoiesis at equilibrium. At time 0, 10⁴ leukemic stem cells (HSC^L, S^L) cells are introduced to the system and subsequently expand over time at the expense of the normal cells, which decrease. (C) Sensitivity analyses of the outcomes of CML hematopoiesis with values corresponding to the proportion of parameter sets where less than 50% of terminal cells are leukemic. (D) The fitness of the leukemic stem cells relative to the normal stem cells, as measured by the ratio of their characteristic self-renewal fractions ( ${\bar{p}}_{0}^{L}$ / ${\bar{p}}_{0}$ ) determines whether CML will progress and leukemic cells will take over the system after CML stem cells are introduced.

Figure 4—figure supplement 1

Download asset Open asset

The dynamics from steady state upon introduction of leukemic stem cells with perturbations of each parameter, as labeled, ranging from 90% to 110% of the original parameter value in Appendix 1—table 5.

The solid line represents the median of the perturbations, while the shaded region represents 95% of the perturbation dynamics.

Figure 4—figure supplement 2

Download asset Open asset

Figure 4—figure supplement 3

Download asset Open asset

Figure 4—figure supplement 4

Download asset Open asset

Figure 4—figure supplement 5

Download asset Open asset

Figure 4—figure supplement 6

Download asset Open asset

Variations in the other leukemic parameters and their associated gamma for all eligible parameter sets.

Plots display data using the same criteria as Figure 4C, but results show nearly complete insensitivity to changes in either value for the other leukemic parameters.

Figure 5 with 3 supplements

Download asset Open asset

Validation of model through simulated transplant.

Results of a transplant experiment from Reynaud et al., 2011. Schematic (A) depicting the experimental pipeline and results (B), adapted from Figure 2A-C in Reynaud et al., 2011. When HSC^L are transplanted into sublethally irradiated mice, chronic myeloid leukemia (CML)-like leukemia is induced and the myeloid cells expand. When leukemic multipotent progenitor cells (MPP^L, P^L) are transplanted, they do not stably engraft and transiently produce a larger fraction of differentiated lymphoid cells. (C) Simulated time evolutions of donor-derived HSC^L, MPP^L, and terminally differentiated lymphoid (TD_l), and myeloid (TD_m) cells when HSC^L (solid) or MPP^L (dashed) cells are transplanted. (D) Bar chart showing model predictions of the percentages of donor-derived myeloid and lymphoid cells after 35 d when HSC^L or MPP^L are transplanted, which is consistent with the experimental data in (B; see text).

Figure 5—figure supplement 1

Download asset Open asset

Heat map depicting the outcomes of transplant experiments in the presence of decrements of 50–70% HSC^L and 30–50% MPP^L from their equilibrium values (see text for details).

The heat map quantifies the number of parameter sets that are consistent with the experimental observations (Reynaud et al., 2011) that there is a simple majority of myeloid or lymphoid cells after 35 d when HSC^L or MPP^L are transplanted, respectively.

Figure 5—figure supplement 2

Download asset Open asset

The dynamics of each parameter set that does not match the experimentally observed behavior of the transplant experiments from Reynaud et al., 2011.

These parameter sets are removed from the pool of eligible parameter sets if in the worst-case scenario for each parameter set they are not consistent with the criteria outlined in Figure 5—figure supplement 1.

Figure 5—figure supplement 3

Download asset Open asset

Distributions of the remaining 478 parameters after removal determined through the depletion sweep.

Figure 6 with 4 supplements

Download asset Open asset

The response of chronic myeloid leukemia (CML) to tyrosine kinase inhibitor (TKI) therapy.

(**A–C**) Simulated cell dynamics of normal and leukemic cells in response to TKI therapy that is started at different time points in CML development (A, early times; B, intermediate times; C, late times). See text. (D) Simulated molecular response curves corresponding to the application of TKIs for each of simulations in (**A–C**). The simulated molecular response from (A) (blue) compares well with clinical data (symbols) measuring treatment responses to two different TKIs (imatinib, dasatinib) averaged across a cohort of patients (Glauche et al., 2018). The simulated molecular responses from (B) and (C) (orange solid and dashed curves) are indicative of primary resistance. (**E–G**) Experiments in chimeric mice (see text) that show that the size of the leukemic stem cell clone correlates with decreased response to TKI therapy. Peripheral blood (PB) leukocyte counts (E), percentage of PB granulocytes (F), and PB *BCR-ABL1*⁺ (leukemic) granulocyte chimerism (G) are shown in cohorts of mice treated with dasatinib. Blue symbols are chimeras bearing >90% *BCR-ABL1*⁺ HSC^L, red-orange symbols are chimeras bearing 46–75% *BCR-ABL1*⁺ HSC^L.

Figure 6—figure supplement 1

Download asset Open asset

Effective parameters for proliferation, self-renewal, and branching after leukemic stem cells are added to the normal system at steady state, and after tyrosine kinase inhibitor (TKI) therapy begins for the cases shown in Figure 6 in the main text.

The vertical gray line depicts the start of treatment (time $t = 0$ ). Black represents treatment at early times (shortly after chronic myeloid leukemia [CML] develops), blue denotes treatment at intermediate times, and green denotes treatment a long time after CML develops.

Figure 6—figure supplement 2

Download asset Open asset

The dynamics from steady state upon introduction of leukemic stem cells with perturbations of each leukemic parameter, as labeled, ranging from 90% to 110% of the original parameter value (Appendix 1—table 4 for normal cells, Appendix 1—table 5 for leukemic cells), except for $p_{0, m a x}^{L}$ , which is never varied above the $p_{0, m a x}$ value.

The solid line represents the unperturbed state, while the shaded region represents the entire range of the perturbation dynamics.

Figure 6—figure supplement 3

Download asset Open asset

Figure 6—figure supplement 4

Download asset Open asset

Figure 7 with 2 supplements

Download asset Open asset

Leukemic stem cell load alone does not predict response to tyrosine kinase inhibitor (TKI) therapy.

(A) Scatter plot of the distribution of simulated *BCR-ABL1* transcripts at 18 mo after start of treatment as a function of the HSC^L proportion at the start of TKI therapy for each of the 478 parameter sets. The time to reach MR3 (*BCR-ABL1* transcripts less than 0.1%) is indicated by the color. (B) The proportion of parameter sets that achieve MR3 within 50 mo (responders, blue) and those that do not (nonresponders, orange) shown as a joint distribution of the parameters $p_{0, m a x}$ and ${\bar{p_{0}}}^{L} / p_{0, m a x}$ using the minimum fitness threshold $\frac{{\bar{p_{0}}}^{L}}{p_{0, m a x}} \geq 0.7$ , which reveals the maximal self-renewal fraction of the normal HSC, $p_{0, m a x}$ (y-axis marginal distribution), distinguishes response across parameter combinations. (C) Dynamics of the effective leukemic stem cell self-renewal fraction for the parameter set used in Figure 6A–D (blue) and an arbitrary representative non-responsive parameter set (orange) during chronic myeloid leukemia (CML) development and before initiation of therapy (t < 0), and after application of TKI therapy (t > 0). (D) Early time dynamics (e.g., t = 0–3 mo; left of the vertical line) of the transcript levels reveal that it is difficult to distinguish responders (blue) from nonresponders (orange). At later times (e.g., t = 3–6 mo; right of the vertical line), the two populations are easier to distinguish. (E) Receiver operating curves (ROC) obtained from the 478 parameter sets using our new prognostic criterion based on the relative changes in transcript levels (solid) and the transcript halving time (dashed) for the first 3 mo (blue) and the second 3 mo (orange) after therapy. The prognostic thresholds (symbols) are identified by optimizing true and false positives rates. Early molecular response (EMR) at 3 mo (10% transcript levels) and 6 mo (1% transcript levels) are shown by the blue square and orange diamond, respectively. Inset: expanded view of the ‘elbow’ region of the ROC curves to display differences between the prognostic tests. Accuracy is improved using the 3–6 mo time window, and our new prognostic criterion outperforms the EMR and halving time prognostics in this time window. (F) The accuracy of the prognostic criteria applied to CML patient data (n = 7) treated using the same TKI dosing for the 6-month period after therapy is started. The results are consistent with the synthetic data in (E). (G) The prognostic criteria applied to patient data (n = 7) in which therapy could be changed but those changes were maintained for 6 mo (see text). Although the data are limited, the results are consistent with those in (E) and (F) suggesting increased accuracy using the 3–6-month window, and that the prognostic criterion based on relative change may yield more accurate predictions than EMR and halving time in the 3–6-month time frame.

Figure 7—figure supplement 1

Download asset Open asset

Parameter distributions by response to tyrosine kinase inhibitor (TKI) treatment using synthetic data.

Blue denotes a parameter set achieves MR3 within 50 mo (responder), and yellow denotes a parameter set that does not (non-responder).

Figure 7—figure supplement 2

Download asset Open asset

Comparison of the proportion of cases that successfully respond to therapy as a function of the ratio of characteristic values of leukemic and normal stem cells: ${\bar{p}}_{0}^{L} / {\bar{p}}_{0}$ .

The relative fitness is the critical determinant of whether chronic myeloid leukemia (CML) cells respond to tyrosine kinase inhibitor (TKI) treatment.

Figure 8 with 2 supplements

Download asset Open asset

Combining tyrosine kinase inhibitor (TKI) therapy with differentiation promoters enhances response to treatment.

(A) The proportion of the 478 parameter sets that achieve MR3 under combined TKI and differentiation therapy depends nonmonotonically on the strength Δ of the differentiation therapy, with the peak response (86.8%) occurring at Δ = 0.24. (B) The maximum stem cell self-renewal fraction for a single ${\bar{p}}_{0}^{L} / p_{0, m a x}$ in Figure 7B (marginal y) with hatching indicating the effects of the combination of TKI and differentiation therapy with Δ = 0.24. Blue hatching indicates nonresponders (who did not achieve MR3) that become responders (achieve MR3) while orange hatching indicates responders that become nonresponders upon combined treatment. Differentiation promoters allow nonresponders to TKI therapy with large self-renewal fractions to reach MR3. The opposite outcome, loss of MR3 in a TKI responder, primarily occurs only at the smallest self-renewal fractions. (C) Time evolution of *BCR-ABL1* transcripts during combination therapy, with Δ = 0.24 (black) and Δ = 0.5 (blue), using the parameter set from Figure 5B that does not achieve MR3 using TKI monotherapy (red). (**D, E**) The time derivatives of the number of normal (D) and leukemic (E) stem cells during combination therapy. The differentiation promoter attenuates the rapid increases in the rates of change at early times after therapy starts in both normal and leukemic cells, but the attenuation is much larger in the leukemic cells. This results in the growth of normal cells, while leukemic cells experience restricted growth or outright depletion depending upon the differentiation therapy strength. (F) Simplified diagram representing the key interactions between the cells and the impact on outcomes of TKI and combination therapy. Green: chronic myeloid leukemia (CML) hematopoiesis depicting the loss of normal stem cells and progenitors and the increase in leukemic stem cells and progenitors. Red fill: TKI treatment failure. The TKI-induced death of leukemic progenitors relieves negative feedback and increases stem cell self-renewal, resulting in increases in both normal and leukemic stem cells, and eventually their progeny (panel 2). The increases are larger for the leukemic cells because their self-renewal fraction is bigger. Increases in the leukemic progenitor compartment (panel 3) drive down the self-renewal fraction of normal stem cells proportionally more than for the leukemic stem cells. The increases in HSC^L also drive down proliferation rates, which makes the leukemic cells less responsive to TKI treatment. Altogether, this makes the leukemic cells more fit than the normal cells and results in therapy failure. Blue fill: treatment by combined TKI and pro-differentiation therapy reduces stem cell self-renewal relative to TKI monotherapies, equalizes the normal and leukemic self-renewal fractions, which limits leukemic stem cell growth and limits decreases in proliferation rates, making the HSC^L and MPP^L more susceptible to TKI-induced death (panel 2). This allows repopulation of the bone marrow by normal stem cells and progenitors to occur (panel 3).

Figure 8—figure supplement 1

Download asset Open asset

Changes in the parameter distributions from Figure 7—figure supplement 1 when combined differentiation and tyrosine kinase inhibitor (TKI) therapy is administered.

Blue denotes a parameter set achieves MR3 within 50 mo, and orange denotes parameter sets that do not. Orange hatching denotes a parameter set that achieves MR3 within 50 mo with TKI monotherapy, but does not achieve MR3 within 50 mo under combination therapy. Blue hatching denotes a parameter set that does achieve MR3 within 50 mo with combination therapy but does not with TKI therapy alone.

Figure 8—figure supplement 2

Download asset Open asset

Variations in $\frac{p_{0, m a x}^{L}}{p_{0, m a x}}$ and $\frac{γ_{1}^{L}}{γ_{1}}$ for all parameter sets with combination therapy show improvement of response with combination therapy compared to Appendix 1—figure 11 (bottom right).

Plots display data using the same criteria as Appendix 1—figure 11 (bottom right).

Appendix 1—figure 1

Download asset Open asset

Flow chart for design space analysis (DSA).

Given an ordinary differential equation (ODE), we can obtain the design space by obtaining all S-systems, steady states, and evaluated dominance conditions. S-systems that do not satisfy the dominance condition are not included in the design space.

Appendix 1—figure 2

Download asset Open asset

Lineage schematic depicting a stem cell and terminally differentiated cell with negative feedback onto stem cell self-renewal probability $p_{0}$ .

Appendix 1—figure 3

Download asset Open asset

Design space for four-cell lineage model.

The design space is showing a slice of parameter space varying the max self-renewal probability ( $\bar{p_{0}}$ ) and feedback gain for stem cells ( $γ$ ). If we are to sample parameters sets (shown by points), we observe oscillatory behavior in the full system, as shown by the time evolution plots.

Appendix 1—figure 4

Download asset Open asset

Pairwise parameter distributions for all 1493 parameter sets found from the gridsearch.

The overall distributions for each parameter are shown along the diagonal.

Appendix 1—figure 5

Download asset Open asset

Pairwise parameter distributions for the 563 parameter sets from Appendix 1—figure 4 that have feedforward gain $γ_{5} > 0.01$ .

The overall distributions for each parameter are shown along the diagonal.

Appendix 1—figure 6

Download asset Open asset

A spaghetti plot showing the dynamics of the effective proliferation, self-renewal, and branching parameters for 50 parameter sets under tyrosine kinase inhibitor (TKI) treatment started at early times.

Appendix 1—figure 7

Download asset Open asset

Comparison between our prognostic and alternate prognostics, as labeled, at the first and second 3 mo after the start of therapy.

Appendix 1—figure 8

Download asset Open asset

The ability of prognostic criteria to predict MR3 by 18 mo is evaluated using anonymized patient data (n = 11) in which patients received the same therapy for 6 mo either from the start of therapy or after a change in therapy.

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. The Early molecular response (EMR) prognostic criterion is not used when the patients have had a therapy change.

Appendix 1—figure 9

Download asset Open asset

The predictive ability of each prognostic criterion from Appendix 1—figure 8 but grouped based upon patient (n = 11) outcome (blue, responder; and yellow, nonresponder) and the time frame for prediction.

Appendix 1—figure 10

Download asset Open asset

Aggregating the accuracy of prognostic criteria predictions with patient data (n = 11) from Appendix 1—figure 8 where a general time of 0–3 mo contains both 0–3 mo after start of therapy and after therapy change.

The 3–6-month data is aggregated similarly.

Appendix 1—figure 11

Download asset Open asset

Variations in $\frac{p_{0, m a x}^{L}}{p_{0, m a x}}$ and $\frac{γ_{1}^{L}}{γ_{1}}$ for the individual parameter set from Figure 6 (top) and all eligible parameter sets (bottom).

Left heat maps display the proportion of parameter sets that maintain dominance of normal cells in the system in the lighter regions and the proportion that possess leukemic cell dominance in the darker regions. The right heat maps indicate proportion of cases where parameter sets achieve MR3 within 50 mo with lighter values associated with a higher proportion of response. Biologically relevant parameter combinations exist within darker regions on the bottom left and regions within orange to purple on the bottom-right plot.

Appendix 1—figure 12

Download asset Open asset

Combinations of leukemic values with similar overall response rates to tyrosine kinase inhibitor (TKI) therapy (see text) yield similar distributions of response.

Appendix 1—figure 13

Download asset Open asset

Simulated distributions of response to tyrosine kinase inhibitor (TKI) therapy as a function of initial proportion of leukemic stem cells is qualitatively similar across the parameter combinations from Appendix 1—figure 12 (from left to right): $\frac{p_{0, m a x}^{L}}{p_{0, m a x}} = 0.8$ and $\frac{γ_{1}^{L}}{γ_{1}} = 0.1$ , $\frac{p_{0, m a x}^{L}}{p_{0, m a x}} = 0.8$ and $\frac{γ_{1}^{L}}{γ_{1}} = 0.2$ , $\frac{p_{0, m a x}^{L}}{p_{0, m a x}} = 0.9$ and $\frac{γ_{1}^{L}}{γ_{1}} = 0.3$ , and $\frac{p_{0, m a x}^{L}}{p_{0, m a x}} = 1$ and $\frac{γ_{1}^{L}}{γ_{1}} = 0.5$ .

Appendix 1—figure 14

Download asset Open asset

Comparison of trajectories of response as functions of ${\bar{p}}_{0}^{L} / p_{0, m a x}$ and $p_{0, m a x}$ .

For each value of $\frac{p_{0, m a x}^{L}}{p_{0, m a x}}$ , there is one tree starting from $p_{0, m a x} = 0.5$ , then for each $\frac{γ_{1}^{L}}{γ_{1}}$ up to $\frac{γ_{1}^{L}}{γ_{1}} = 0.6$ the tree splits into six branches. Dot color denotes whether a parameter set responds (blue) or does not respond (orange).

Appendix 1—figure 15

Download asset Open asset

Bivariate histogram for response to tyrosine kinase inhibitor (TKI) therapy where the relative fitness of $H S C^{L}$ is ${\bar{p}}_{0}^{L} / {\bar{p}}_{0} \geq 0.50$ .

Individual variable marginal distributions are shown along the sub-axes.

Appendix 1—figure 16

Download asset Open asset

Bivariate histogram for response to tyrosine kinase inhibitor (TKI) therapy where the relative fitness of $H S C^{L}$ is ${\bar{p}}_{0}^{L} / {\bar{p}}_{0} \geq 0.60$ .

Individual variable marginal distributions are shown along the sub-axes.

Appendix 1—figure 17

Download asset Open asset

Prognostic value sweep applied to combinations of $\frac{p_{0, m a x}^{L}}{p_{0, m a x}}$ and $\frac{γ_{1}^{L}}{γ_{1}}$ with leukemic predominance > 60%.

Optimal values from Figure 7E are no longer the optimal across all combinations, but the 3–6-month time frame still outcompetes 0–3-month time frame.

Tables

Appendix 1—table 1

A sample of S-systems from Equations 14–18.

	S-system 1	S-system 2	S-system 3	S-system 4
$x_{S}^{'}$	$2 \bar{p_{0}} x_{n e w}^{- 1} η_{1} x_{S} - η_{1} x_{S}$	$2 \bar{p_{0}} x_{n e w}^{- 1} η_{1} x_{S} - η_{1} x_{S}$	$2 \bar{p_{0}} x_{n e w}^{- 1} η_{1} x_{S} - η_{1} x_{S}$	$2 \bar{p_{0}} x_{n e w}^{- 1} η_{1} x_{S} - η_{1} x_{S}$
$x_{P}^{'}$	$2 η_{1} x_{S} - 2 \bar{p_{0}} x_{n e w}^{- 1} η_{1} x_{S}$	$2 η_{1} x_{S} - η_{2} x_{S}$	$2 p_{1} η_{2} x_{S} - 2 \bar{p_{0}} x_{n e w}^{- 1} η_{1} x_{S}$	$2 p_{1} η_{2} x_{S} - η_{2} x_{S}$
$x_{T D_{l}}^{'}$	$2 q_{1} η_{2} x_{P} - d_{l} x_{T D_{l}}$	$2 q_{1} η_{2} x_{P} - d_{l} x_{T D_{l}}$	$2 q_{1} η_{2} x_{P} - d_{l} x_{T D_{l}}$	$2 q_{1} η_{2} x_{P} - d_{l} x_{T D_{l}}$
$x_{T D_{m}}^{'}$	$2 η_{2} x_{P} - d_{m} x_{T D_{m}}$	$2 η_{2} x_{P} - d_{m} x_{T D_{m}}$	$2 η_{2} x_{P} - d_{m} x_{T D_{m}}$	$2 η_{2} x_{P} - d_{m} x_{T D_{m}}$
0	$γ x_{P} - x_{n e w}$	$γ x_{P} - x_{n e w}$	$γ x_{P} - x_{n e w}$	$γ x_{P} - x_{n e w}$
Boundary 1	$2 η_{1} x_{S} > 2 p_{1} η_{2} x_{S}$	$2 η_{1} x_{S} > 2 p_{1} η_{2} x_{S}$	$2 η_{1} x_{S} < 2 p_{1} η_{2} x_{S}$	$2 η_{1} x_{S} < 2 p_{1} η_{2} x_{S}$
Boundary 2	$2 \bar{p_{0}} x_{n e w}^{- 1} η_{1} x_{S} > η_{2} x_{S}$	$2 \bar{p_{0}} x_{n e w}^{- 1} η_{1} x_{S} < η_{2} x_{S}$	$2 \bar{p_{0}} x_{n e w}^{- 1} η_{1} x_{S} > η_{2} x_{S}$	$2 \bar{p_{0}} x_{n e w}^{- 1} η_{1} x_{S} < η_{2} x_{S}$

Appendix 1—table 2

Parameter values used to model the normal hematopoietic system.

Parameter	Description
$p_{0, m a x}$	Maximal self-renewal fraction of HSC
$p_{1, m a x}$	Maximal self-renewal fraction of MPP
$q_{1, m a x}$	Maximum branching fraction from MPP to $T D_{L}$
$η_{1, m a x}$	Maximum HSC proliferation rate
$η_{2, m a x}$	Maximum MPP proliferation rate
$γ_{1}$	Feedback gain on HSC self-renewal fraction
$γ_{2}$	Feedback gain on HSC proliferation rate (from HSC)
$γ_{3}$	Feedback gain on MPP self-renewal fraction
$γ_{4}$	Feedback gain on MPP branching fraction
$γ_{5}$	Feedforward gain on MPP proliferation rate
$d_{L}$	Death rate of $T D_{L}$
$d_{m}$	Death rate of $T D_{M}$

HSC: hematopoietic stem cell; MPP: multipotential progenitor.

Appendix 1—table 3

Additional parameter values used to model the chronic myeloid leukemia (CML) cell population dynamics and effect of tyrosine kinase inhibitor (TKI) therapy.

Parameter	Description
$p_{0, m a x}^{L}$	Maximal self-renewal fraction of HSC $^{L}$
$p_{1, m a x}^{L}$	Maximal self-renewal fraction of MPP $^{L}$
$q_{1, m a x}^{L}$	Maximum branching fraction from MPP $^{L}$ to $T D_{L} ^{L}$
$η_{1, m a x}^{L}$	Maximum HSC $^{L}$ proliferation rate
$η_{2, m a x}^{L}$	Maximum MPP $^{L}$ proliferation rate
$γ_{1}^{L}$	Feedback gain on HSC $^{L}$ self-renewal fraction
$γ_{2}^{L}$	Feedback gain on HSC $^{L}$ proliferation rate (from HSC $^{L}$ )
$γ_{3}^{L}$	Feedback gain on MPP $^{L}$ self-renewal fraction
$γ_{4}^{L}$	Feedback gain on MPP $^{L}$ branching fraction
$γ_{5}^{L}$	Feedforward gain on MPP $^{L}$ proliferation rate
$d_{L}^{L}$	Death rate of $T D_{L} ^{L}$
$d_{m}^{L}$	Death rate of $T D_{M} ^{L}$
$T K I_{H S C}$	Death rate of HSC $^{L}$ due to TKI therapy
$T K I_{M P P}$	Death rate of MPP $^{L}$ due to TKI therapy

HSC: hematopoietic stem cell; MPP: multipotential progenitor.

Appendix 1—table 4

Parameter values used to model the normal hematopoietic system in Figures 3—7 in the main text.

Parameter	Value
$p_{0, m a x}$	0.756641
$p_{1, m a x}$	0.357913
$q_{1, m a x}$	0.032241
$η_{1, m a x}$	0.197639
$η_{2, m a x}$	0.47088
$γ_{1}$	0.513281
$γ_{2}$	0.197639
$γ_{3}$	0.47088
$γ_{4}$	0.513281
$γ_{5}$	0.543987
$d_{L}$	0.000165
$d_{m}$	0.000428

Appendix 1—table 5

Additional parameter values used to model the chronic myeloid leukemia (CML) cell population dynamics in Figures 4A ,B—7 in the main text.

Parameter	Value
$p_{0, m a x}^{L}$	0.756641
$p_{1, m a x}^{L}$	0.357913
$q_{1, m a x}^{L}$	0.032241
$η_{1, m a x}^{L}$	0.197639
$η_{2, m a x}^{L}$	0.47088
$γ_{1}^{L}$	0.2566405
$γ_{2}^{L}$	0.197639
$γ_{3}^{L}$	0.47088
$γ_{4}^{L}$	0.2566405
$γ_{5}^{L}$	0.543987
$d_{L}^{L}$	0.000165
$d_{m}^{L}$	0.000428
$T K I_{H S C}$	0.201311
$T K I_{M P P}$	0.024757

Appendix 1—table 6

Secondary parameter values used as a representative case to model nonresponsive chronic myeloid leukemia (CML) cell population dynamics in Figures 7C, D, 8 in the main text.

Parameter	Value
$p_{0, m a x}$	0.838481
$p_{1, m a x}$	0.009776
$q_{1, m a x}$	0.418627
$η_{1, m a x}$	0.226025
$η_{2, m a x}$	0.247591
$γ_{1}$	0.676962
$γ_{2}$	0.000219
$γ_{3}$	0.000326
$γ_{4}$	0.367564
$γ_{5}$	0.168132
$d_{L}$	0.071273
$d_{m}$	0.155649
$p_{0, m a x}^{L}$	0.838481
$p_{1, m a x}^{L}$	0.009776
$q_{1, m a x}^{L}$	0.418627
$η_{1, m a x}^{L}$	0.226025
$η_{2, m a x}^{L}$	0.247591
$γ_{1}^{L}$	0.338481
$γ_{2}^{L}$	0.000219
$γ_{3}^{L}$	0.000326
$γ_{4}^{L}$	0.367564
$γ_{5}^{L}$	0.168132
$d_{L}^{L}$	0.071273
$d_{m}^{L}$	0.155649
$T K I_{H S C}$	0.25
$T K I_{M P P}$	30

Additional files

MDAR checklist: https://cdn.elifesciences.org/articles/84149/elife-84149-mdarchecklist1-v5.docx
Download elife-84149-mdarchecklist1-v5.docx

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Mendeley

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

Jonathan Rodriguez
Abdon Iniguez
Nilamani Jena
Prasanthi Tata
Zhong-Ying Liu
Arthur D Lander
John Lowengrub
Richard A Van Etten

(2023)

Predictive nonlinear modeling of malignant myelopoiesis and tyrosine kinase inhibitor therapy

eLife 12:e84149.

https://doi.org/10.7554/eLife.84149

Share this article

Cite this article

Branched lineage model of normal hematopoiesis with feedback regulation.

Fluorescence-activated cell sorting (FACS) analysis of mouse Lin–Sca-1+c-Kit+ (LSK) bone marrow stem/progenitor cells and the proposed branched lineage hematopoiesis model.

Qualitative behavior of feedforward-feedback model and parameter sensitivity.

The dynamics from the steady state with perturbations of each parameter, as labeled, ranging from 90% to 110% of the original parameter value from Appendix 1—table 4.

Effective parameters for proliferation, self-renewal, and branching as the solutions to the model for the normal hematopoietic system approach steady state.

Extension of the model of hematopoiesis to chronic myeloid leukemia (CML).

The dynamics from steady state upon introduction of leukemic stem cells with perturbations of each parameter, as labeled, ranging from 90% to 110% of the original parameter value in Appendix 1—table 5.

The dynamics from steady state upon introduction of leukemic stem cells with perturbations of each parameter, as labeled, ranging from 90% to 110% of the original parameter value in Appendix 1—table 5.

The dynamics from steady state upon introduction of leukemic stem cells with perturbations of each parameter, as labeled, ranging from 90% to 110% of the original parameter value in Appendix 1—table 5.

The dynamics from steady state upon introduction of leukemic stem cells with perturbations of each parameter, as labeled, ranging from 90% to 110% of the original parameter value in Appendix 1—table 5.

The dynamics from steady state upon introduction of leukemic stem cells with perturbations of each parameter, as labeled, ranging from 90% to 110% of the original parameter value in Appendix 1—table 5.

Variations in the other leukemic parameters and their associated gamma for all eligible parameter sets.

Validation of model through simulated transplant.

Heat map depicting the outcomes of transplant experiments in the presence of decrements of 50–70% HSCL and 30–50% MPPL from their equilibrium values (see text for details).

The dynamics of each parameter set that does not match the experimentally observed behavior of the transplant experiments from Reynaud et al., 2011.

Distributions of the remaining 478 parameters after removal determined through the depletion sweep.

The response of chronic myeloid leukemia (CML) to tyrosine kinase inhibitor (TKI) therapy.

Effective parameters for proliferation, self-renewal, and branching after leukemic stem cells are added to the normal system at steady state, and after tyrosine kinase inhibitor (TKI) therapy begins for the cases shown in Figure 6 in the main text.

The dynamics from steady state upon introduction of leukemic stem cells with perturbations of each leukemic parameter, as labeled, ranging from 90% to 110% of the original parameter value (Appendix 1—table 4 for normal cells, Appendix 1—table 5 for leukemic cells).

The dynamics from steady state upon introduction of leukemic stem cells with perturbations of each leukemic parameter, as labeled, ranging from 90% to 110% of the original parameter value (Appendix 1—table 4 for normal cells, Appendix 1—table 5 for leukemic cells).

Leukemic stem cell load alone does not predict response to tyrosine kinase inhibitor (TKI) therapy.

Parameter distributions by response to tyrosine kinase inhibitor (TKI) treatment using synthetic data.

Comparison of the proportion of cases that successfully respond to therapy as a function of the ratio of characteristic values of leukemic and normal stem cells: p¯0L/p¯0.

Combining tyrosine kinase inhibitor (TKI) therapy with differentiation promoters enhances response to treatment.

Changes in the parameter distributions from Figure 7—figure supplement 1 when combined differentiation and tyrosine kinase inhibitor (TKI) therapy is administered.

Variations in p0,maxLp0,max and γ1Lγ1 for all parameter sets with combination therapy show improvement of response with combination therapy compared to Appendix 1—figure 11 (bottom right).

Flow chart for design space analysis (DSA).

Lineage schematic depicting a stem cell and terminally differentiated cell with negative feedback onto stem cell self-renewal probability p0 .

Design space for four-cell lineage model.

Pairwise parameter distributions for all 1493 parameter sets found from the gridsearch.

Pairwise parameter distributions for the 563 parameter sets from Appendix 1—figure 4 that have feedforward gain γ5>0.01.

A spaghetti plot showing the dynamics of the effective proliferation, self-renewal, and branching parameters for 50 parameter sets under tyrosine kinase inhibitor (TKI) treatment started at early times.

Comparison between our prognostic and alternate prognostics, as labeled, at the first and second 3 mo after the start of therapy.

The ability of prognostic criteria to predict MR3 by 18 mo is evaluated using anonymized patient data (n = 11) in which patients received the same therapy for 6 mo either from the start of therapy or after a change in therapy.

The predictive ability of each prognostic criterion from Appendix 1—figure 8 but grouped based upon patient (n = 11) outcome (blue, responder; and yellow, nonresponder) and the time frame for prediction.

Aggregating the accuracy of prognostic criteria predictions with patient data (n = 11) from Appendix 1—figure 8 where a general time of 0–3 mo contains both 0–3 mo after start of therapy and after therapy change.

Variations in p0,maxLp0,max and γ1Lγ1 for the individual parameter set from Figure 6 (top) and all eligible parameter sets (bottom).

Combinations of leukemic values with similar overall response rates to tyrosine kinase inhibitor (TKI) therapy (see text) yield similar distributions of response.

Comparison of trajectories of response as functions of p¯0L/p0,max and p0,max .

Bivariate histogram for response to tyrosine kinase inhibitor (TKI) therapy where the relative fitness of HSCL is p¯0L/p¯0≥0.50.

Bivariate histogram for response to tyrosine kinase inhibitor (TKI) therapy where the relative fitness of HSCL is p¯0L/p¯0≥0.60.

Prognostic value sweep applied to combinations of p0,maxLp0,max and γ1Lγ1 with leukemic predominance > 60%.

A sample of S-systems from Equations 14–18.

Parameter values used to model the normal hematopoietic system.

Additional parameter values used to model the chronic myeloid leukemia (CML) cell population dynamics and effect of tyrosine kinase inhibitor (TKI) therapy.

Parameter values used to model the normal hematopoietic system in Figures 3—7 in the main text.

Additional parameter values used to model the chronic myeloid leukemia (CML) cell population dynamics in Figures 4A ,B—7 in the main text.

Secondary parameter values used as a representative case to model nonresponsive chronic myeloid leukemia (CML) cell population dynamics in Figures 7C, D, 8 in the main text.

MDAR checklist

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)

Heat map depicting the outcomes of transplant experiments in the presence of decrements of 50–70% HSC^L and 30–50% MPP^L from their equilibrium values (see text for details).

Comparison of the proportion of cases that successfully respond to therapy as a function of the ratio of characteristic values of leukemic and normal stem cells: ${\bar{p}}_{0}^{L} / {\bar{p}}_{0}$ .

Variations in $\frac{p_{0, m a x}^{L}}{p_{0, m a x}}$ and $\frac{γ_{1}^{L}}{γ_{1}}$ for all parameter sets with combination therapy show improvement of response with combination therapy compared to Appendix 1—figure 11 (bottom right).

Lineage schematic depicting a stem cell and terminally differentiated cell with negative feedback onto stem cell self-renewal probability $p_{0}$ .

Pairwise parameter distributions for the 563 parameter sets from Appendix 1—figure 4 that have feedforward gain $γ_{5} > 0.01$ .

Variations in $\frac{p_{0, m a x}^{L}}{p_{0, m a x}}$ and $\frac{γ_{1}^{L}}{γ_{1}}$ for the individual parameter set from Figure 6 (top) and all eligible parameter sets (bottom).

Comparison of trajectories of response as functions of ${\bar{p}}_{0}^{L} / p_{0, m a x}$ and $p_{0, m a x}$ .

Bivariate histogram for response to tyrosine kinase inhibitor (TKI) therapy where the relative fitness of $H S C^{L}$ is ${\bar{p}}_{0}^{L} / {\bar{p}}_{0} \geq 0.50$ .

Bivariate histogram for response to tyrosine kinase inhibitor (TKI) therapy where the relative fitness of $H S C^{L}$ is ${\bar{p}}_{0}^{L} / {\bar{p}}_{0} \geq 0.60$ .

Prognostic value sweep applied to combinations of $\frac{p_{0, m a x}^{L}}{p_{0, m a x}}$ and $\frac{γ_{1}^{L}}{γ_{1}}$ with leukemic predominance > 60%.