Dynamics of competing heterogeneous clones in blood cancers explains multiple observations - a mathematical modeling approach

1.   Introduction

Hematopoietic stem cells reside in the bone marrow and are responsible for continuously replenishing the blood cells; red blood cells, white blood cells, or platelets. The hematopoietic stem cells possess the abilities of self-renewal and extensive proliferation which are essential for the maintenance of the formation of blood cells, hematopoiesis ^[1].

The proliferation rate for each of the $10^4-10^5$ normal hematopoietic stem cells is about once per every $25$th to $50$th week ^[2,3]. Cell numbers are amplified through a sequence of maturation stages, with estimated numbers of mitotic events being 31 ^[4]. As a result approximately one million new mature blood cells are generated per second corresponding to $10^{11}$ mature blood cells daily which is comparable to the number of stars in the Milky Way ^[5]. Evidently, a tight regulation of hematopoiesis is crucial and disturbances to the regulation of the stem cell cycle may have severe consequences ^[6]. Several hematopoietic diseases arise from mutations in the hematopoietic stem cells causing the regulation of hematopoiesis to be disturbed. Such distortion may cause disturbances or uncontrolled growth of dysfunctional blood cells resulting in diseases, e.g., blood cancer.

Approximately 175, 000 Americans are each year diagnosed with myeloid malignancies, e.g., acute myeloid leukemia (AML), myelodysplastic syndromes (MDS), chronic myeloid leukemia (CML) and the Philadelphia-negative myeloprofilerative neoplasms (MPNs). MPNs develop slowly on a time scale of 15-25 years whereas other leukemias, like AML and CML develop on a time scale of 1 year or even on a shorter time scale ^[7].

MPNs associated with the Philadelphia-negative mutations are driven by the JAK2V617F (hencefort JAK2), MPL, CALR or TET2 mutations ^[8]. MPNs mutations relate to the diseases; essential thrombocythemia (ET), polycythemia vera (PV) and primary myelofibrosis (PMF) ^[9]. The chronic MPNs are characterized by excessive myeloproliferation, which cause disturbance of hematopoiesis ^[10]. As a consequence, patients with MPNs have a high risk of thrombosis with cardiovascular complications and increased propensity to develop autoimmune and chronic inflammatory diseases ^[11,12,13]. Likewise, patients with chronic inflammatory disease have an elevated risk of an MPN disease and other cancers ^[14,15]. Moreover, MPNs may ultimately develop into AML which is a more rapidly progressing disease with a poor prognosis ^[16].

1.1.   Purpose of the present study

We present a proof of concept model referred to as the Multiple Clone Cancitis (MCC) model allowing for few as well as multiple clones in a flexible mathematical framework. The MCC model describes; the interaction of several types of competing stem cells in the bone marrow, whereof some may be malignant, their mature offspring in the blood circulation, the amount of dead cells and the activity of the immune system degrading the debris from the dead cells as well as the cytotoxic T-cell response of the adaptive immune system (see Figure 1 for a conceptual diagram of the MCC model).

The model is a generalization of the model presented by Andersen et al. in 2017 ^[17], analyzed in ^[18] and investigated using a reduction approach in ^[19], expanded to include the adaptive T-cell response explicitly by Ottesen et al. in 2019 ^[20], and used for prediction of a randomized trial in Ottesen 2020 ^[21]. In contrast to its predecessors, the presented MCC model allows for multiple competing clones and thereby enables reproduction of complex dynamics which its predecessors cannot account for.

We use the properties of the micro-environment, i.e., the cell to cell inhibition and immunoediting, to explain very different clinical observations in a diversity of hematopoietic disorders using the MCC model. In particular, we focus on the MCC models ability to reproduce the clinical observations;

● Observed fast oscillations in cell count of time scale 1 month for MPN patients.

● Consequences of CHIP leading to good, fair, or poor responders to treatment.

● Observed cell count during hydroxyurea treatment of MPNs.

● Elimination, equilibrium, and escape phases related to immunoediting.

● Developed resistance to treatment.

● Observed slow oscillations in cell count of time scale 2-3 months for CML patients.

Mathematical modeling of the above mentioned bullets have, to our knowledge, not been investigated before in a unified framework. By additionally including changes in the self-renewal rates the model also explains the evolution of clones over a life span of a human. The structure of the MCC model is chosen deliberately to ensure that the MCC model inherits the previous findings of the Cancitis model, e.g., observed non-oscillatory cell count progression and observed levels of cancer specific related cytokines, c-reactive proteins, and lactic dehydrogenase. Hence, the MCC model provides an agile model framework with the ability to unify all phenomena described.

1.2.   Previous studies

As early as in 1970 Kennedy reported cyclic oscillations (period 35-40 days) in leukocyte counts in CML during hydroxyurea therapy ^[22]. Chikkappa et al. confirmed the finding in 1976 (period 53-69 days) and in addition similar oscillations were also observed in platelets, neutrophils, eosinophils, basophils and reticulocytes which continued for a period after off anti-leukemic therapy ^[23]. In a similar investigation from 1999 Fortin and Mackey estimated the period to 37-83 days ^[24]. Later Hirrayama et al. found oscillations in platelets and white blood cells in a patient diagnosed with CML (period 8 weeks corresponding to 244 days) ^[25]. More recently, Clapp et al. in 2015 and 2018 suggested such oscillations to be related to an inhibitory immune response by the amount of mature cells in a modeling study ^[26,27]. In 2019 Knauer et al. ^[28] used a mathematical model to achieve oscillations with a 20 day period fitting data from a patient with cyclic neutropenia. Considering patients with JAK2 positive MPNs (JAK2$^+$ MPNs) faster oscillations (period 27 days on average) have also been observed clinically by several studies ^[29,30,31].

Competition or selection of clones as a driver for cancer evolution has been of long term interest. In this study we consider a clone as a group of cells having identical genome, i.e., DNA. Different clones may arise from mutations altering the genome and hereby constituting another clone. Such mutated cells may or may not be malignant and generally multiple sequential mutations are believed to take place before malignant cells arise. The frequency of somatic mutations in the hematopoietic stem cells increases many-folds with age; from an insignificant level for age less than 20 years to a frequency around 10% at age above 70 ^[32]. Based on a meta study, Heuser et al. 2016 concluded that clonal hematopoiesis of indeterminate potential (CHIP) is common and the frequency increases with age to around 10% at age 70-80 ^[33]. Park et al. investigated the role of CHIP and its role for stem cell transplantation ^[34]. Similar cytogenetic clonal evolution in relation to CML was detected by Cortes et al. in 2004 over a span of different disease burdens ^[35].

Mutations may accumulate during growth, resulting in heterogeneous malignant cells referred to as cancer heterogeneity ^[36]. Walter et al. reported that blood samples for AML patients may contain persistent antecedent founding of a clone containing 182 to 660 somatic mutations. Subclones harboring up to hundreds of mutations may emerge or outgrow the founding clone ^[37]. Ding et al. attributed relapse in AML treated patients to the existence of various clones which may lie latently, possible as a small amount of cells ^[38]. Cancer heterogeneity affect the response to therapy and can present a great challenge for cancer treatment ^[39]. Stiehl et al. investigated the role of heterogeneity of AML clones and resistance to therapy using a mathematical model ^[40].

Clinical experiences show that therapy may work well on a subetaoup having MPNs but not on others for unknown reasons. These subetaoups are denoted good and poor responders, respectively. A third subetaoup responds partially on the treatment, giving rise to coexistent equilibrium (or possibly an coexisting oscillatory) states where malignant clones can coexist with healthy clones for long time.

The three Es of cancer immunoediting describe elimination, equilibrium, and escape phases in cancer development ^[41]. In the presented work the three Es of immunoediting are related to CHIP as well as to the coexisting steady or oscillate states becoming resistant or gaining additional mutations which are less immuno-suppressible.

2.   Materials and methods

2.1.   The Multiple Clone Cancitis Model

To understand complex systems it is preferable first to understand simple related systems which can be done either by in vitro or in silico experiments. This approach holds for investigating cancers too, where a single isolated or two competing cell types are often studied.

Inspired by clinical observations, the MCC model has been formulated in attempt to gain an understanding of the more complex dynamics such as resistance and oscillations.

The mutated stem cells (MSCs) include malignant as well as non-malignant stem cells. The malignant stem cells share the characteristics of the wild type stem cells (WSCs) but may have acquired the ability of 'indefinite' self-renewal ^[42,43] and may ultimately suppress the production of healthy blood cells ^[6]. The corresponding mature blood cells are shortly referred to as hematopoietic mature cells consisting of both wild type mature cells and mutated mature cells whereof some may be cancerous mature cells. Distortion in the number of mature blood cells may affect many vital functions such as fighting off infections, preventing serious bleeding or causing thrombosis.

A conceptual diagram of the MCC model is depicted in Figure 1. The WSCs can self renew, differentiate, die or mutate to MSCs with altered properties such as the malignant stem cells, which behave similar to the WSCs. The hematopoietic stem cell niches provide negative feedback, restoring stem cell numbers upon perturbations. Mature blood cells are the outcome of multiple mitosis events. For brevity we omit intermediate maturation stages, normally represented by approximately 30 generations of proliferative progenitor cells, and let the differentiation output of hematopoietic stem cells enter the production of mature cells with a multiplication factor corresponding to the amplification of intermediate cell divisions. The apoptotic cells have to be cleared by immune cells which provide feedback to the hematopoietic stem cell production. Prior to a malignant mutation, the model is calibrated such that typical clinical values are reproduced, e.g., stem cell numbers, mature cell numbers, time for restoration upon perturbation and well established life time of blood cells. Inflammation - a hallmark of cancer - is known to provide further growth advantage to mutated cells, thus the stem cell niche and the inflammatory level are included to investigate this effect ^[44].

The MSCs may share the characteristics of malignant stem cells and thereby cause cancer heterogeneity. Cancer heterogeneity are often believed to arise through a sequential mutation process, since an agammaessive malignant cell type is likely to arise only from such a sequence of mutations. A conceptual diagram of sequential mutations is illustrated in Figure S.2 where the WSCs mutate into an MSC line which then sequentially mutates further $n$ times. Likewise, a conceptual model of parallel mutations is illustrated in Figure S.1, where WSCs mutate multiple times in parallel into distinct MSCs. Each cell after multiple mutations may be the result of the mutations arriving in any order, thus resulting in cells with identical DNA. However, the phenotype of the cells may differ causing the cells to function differently, i.e., having different RNA activities.

Note that in Figure 1, each hematopoietic cell line can either be the result of a sequential or parallel mutation. Thus a very large class of mutational landscapes is possible due to these phenomena. In this paper, we will mainly consider the two limits, either the WSCs mutate into different MSCs (parallel mutations), or the WSCs mutates into a single MSC line which initiate a sequence of mutations (sequential mutations). The last mutation scenario is considered to happen more frequently than the first, but in reality, it is expected that an intermediate scenario is the most realistic one.

In simulations we may randomly choose an existing cell to mutate and randomly choose the type of mutation in each sufficiently small time intervals. Doing a Monte Carlo simulation of such random process gives information about mean and variation but less mechanistic insight. However, exhausting simulation studies reveal that the points outlined in the present paper are robust with respect to which type of mutation we choose. Thus we mainly discuss the deterministic sequential MCC model to ease the explanatory insight gained.

2.2.   A quick guide to equations - the Multiple Clone Cancitis Model

The model equations arise from formulating the biological mechanisms in subsection 2.1 directly into nonlinear, ordinary, differential equations. A differential equation is formulated for each of the cell types; WSC ($y_0$), MSC ($y_{1}$, $y_{2}$, ..., $y_{n}$) and their mature blood cell counterparts ($z_0, z_1, \ldots, z_n$), dead cells ($a$) and the immune cells ($s$). The default parameters are listed in Table A1. Furthermore, a useful scaling of each of the variables is introduced to bring the model into dimensionless form which is a common approach for a mathematical analysis of such equations. The chosen scaling constants are listed in Table A2.

In line with the biological knowledge, the mature blood cells, immune cells and dead cells quickly equilibrate with the slower stem cell dynamics. This implies that the mature cells, dead cells and immune cells can be explicitly expressed in terms of WSC and the $n$ MSC types, thus the model is reduced significantly from $2(n+1)+2$ to $n+1$ equations. Note that $n$ can can be chosen as any positive integer and the model can therefore be used for any number of mutations. Let the dimensionless cell types and time scale be denoted by $Y_0, Y_1, \ldots, Y_n$ and $T$, respectively. The resulting dimensionless equations are (see Supplementary for derivation).

with the immune cells, $S$, expressed as

$J$ is a parameter describing the inflammatory stimuli not caused by the mutated cells, e.g., exogenous inflammatory stimuli. Increasing this value leads to an increased production of immune cells. $D_{i}$ describes the death rate of cell type $i$ and $K_{i}$ is a death rate of the second order elimination of mutant $i$ due to activated immune clearance of the mutants by cytotoxic T-cells. $B_i$ describes the balance between loss of cells from the stem cell pools and elimination from the dead cell pool due to immune cell clearance.

To maintain function, stem cells must interact with the stem cell niche. The stem cell niche uses negative feedback mechanisms to stabilize wild type stem cells under normal circumstances, e.g., in absence of cancer mutants. When several stem cell types are present, there is a competition for space and resources through the stem cell niche feedback (the $\Phi$-function),

where the weighting coefficient $C_{i, j}$ describes the inhibiting strength of cell type $j$ onto cell type $i$ for $i, j\in\{0, 1, 2, \ldots, n\}$.

A way to include the fact that cancer cells are less constrained by environmental signals than hematopoietic stem cells is to enforce $\Phi_{i} > \Phi_{0}$ for $i = 1, \ldots, n$, i.e., large weighting coefficients imply high niche feedback inhibition strengths. In the current work we are exploring the $C$-dependence on the overall dynamics of cell populations. Finally, $R_i$ denote the self renewal rates of the mutated cells relative to the wild type.

Stochastic or continuous mutation processes between stem cell types are incorporated through $\Xi_i$. Hence, the first mutant is initiated by a wild type mutation, while secondary mutations may happen sequentially (the mutant mutates further), in parallel (a new mutation from the wild type cell pool), or as a mixture of these, where any of the existing stem cells wild type or mutated ones may mutate. The stochastic process is a scaled Poisson process but in the deterministic limit $\Xi_i = \sum_{j = 0}^n ({M_{i, j} Y_j}-{M_{j, i} Y_i})$ where $M_{i, j}$ describes the rate that cell type $j$ mutates into cell type $i$. For the parallel case the expression for $\Xi_i$ reduces to $\Xi_0 = \sum_{j = 1}^n -M_{j, 0}Y_0$ and $\Xi_i = M_{i, 0}Y_0$ for $i = 1, 2, \ldots, n$, while for the sequential case the expressions are given by $\Xi_0 = -M_{1, 0}Y_0$ and $\Xi_i = M_{i, i-1}Y_{i-1}-M_{i+1, i}Y_{i}$ for $i = 1, 2, \ldots, n$ where $M_{n, n+1} = 0$, i.e., the last mutation cannot initiate further differentiation.

By presenting the equations in dimensionless form as in (2.1), the death rates as well as the self renewal rates for WSCs become normalized to 1, i.e., $D_0 = 1$ and $R_0 = 1$, respectively. Likewise, the self niche feedback becomes 1, i.e., $C_{i, i} = 1$. The advantage in putting the equations in dimensionless form is not only to ease the analysis, but also to make reliable simulations and promote insights. However, the results of the in silico investigations are presented using the dimensional variables ($y$ and $z$) instead of the corresponding dimensionless ones ($Y$ and $Z$), while the parameters will be listed in dimensionless form.

2.3.   In silico methodology

For the simulations we will investigate how clones with distinct micro-environment can give rise to complex dynamics observed in clinical studies such as resistance and oscillations. The properties of the micro-environment are modeled by the cell to cell inhibition ($C$) and immunoediting ($K$) and these are the parameters of main interest in the simulations while keeping the remaining parameters fixed at their default values listed in Table A3. In sections 3.7 and 3.1 we will also consider the effect of the self renewal rate $R$.

The default parameter values for the models listed in Table A3 are for the JAK2$^+$ MPNs ^[20] and for each experiment, only the parameter values deviating from the default value will be reported.

At this point we emphasize that the structure of the model is independent of the specific blood cancer that a given mutation may cause. However, the parameter values depend on the specific mutations causing the cancer to develop. Thus, the model structure is general but the parameter values are mutation-specific taking all involved clones into account.

3.   Results

Before turning our attention to the results of the MCC model, we briefly outline the results for the scenario with one or two clones which is described by the MCC models predecessors.

For the simple case of a single clone, the dynamical behavior of the model is simple - either the cell number approaches the carrying capacity or the cell type goes extinct. The case of two competing cell types, i.e., the wild type and one mutant, has been extensively studied in ^[18,20] illuminating cancer elimination, escape or equilibrium depending on the parameter values. No oscillatory dynamics has been observed in the studies, only trajectories approaching steady state values. For JAK2$^+$ MPNs the parameter values were estimated to reproduce observed cell counts, see ^[17,20]. Remarkably, the immune response reproduces observed average values for cohorts.

Now we will present the simulations for the MCC model which allow for the number of clones to exceed two and thereby enabling more complex dynamics than its predecessors can describe.

3.1.   Additional mutation explains elimination, equilibrium and escape

The first investigation is for the case with two distinct MSCs present in the bone marrow in addition to the WSC. The aim is to simulate a scenario where a less agammaessive MSC is reduced by treatment and thereby enabling a more agammaessive MSC to initiate a fatal growth, i.e. wild type cells get extinct.

The two MSCs differ in their susceptibility to the immune response and by their self-renewal rates. The second malignant cell line $y_2$ is more agammaessive as it has a lower susceptibility to the immune response and a higher self-renewal rate. However, $y_{2}$ is suppressed by the presence of $y_{1}$, corresponding to a scenario where $y_1$ could be an earlier mutation with advantageous growth conditions due to the lack of presence of additional mutations in the initial growth phase.

Figure 2a depicts the scenario where the MSCs reach a co-existing steady state with the WSCs. The malignant cells, $y_{2}$ cannot initialize a fatal growth as it is sufficiently suppressed by $y_{1}$. Note there is a non-fatal reduction in the WSCs which could lead to the onset of targeted treatment of $y_{1}$, i.e., treatment that does only affect $y_1$. This is exactly the scenario explored in Figure 2b, where the non-fatal reduction in $y_0$ triggers a treatment of the virtual patient which boosts the immune response to $y_{1}$ (increasing $K_1$) at year 40. The treatment results in a reduction of the amount of $y_{1}$, which initially gives an improved prognosis as the amount of WSCs ($y_0$) increases. However, the inhibiting effect for the more agammaessive malignant cells $y_{2}$ wears off as the amount of $y_{1}$ reduces, and a fatal malignant growth of $y_{2}$ begins.

3.2.   Resistance causes fatal growth

The second simulation scenario explores the concept of resistance. At time $t = 0$ there is a WSC line $y_0$ and single MSC line $y_1$. The two cell lines enter a co-existing steady state as the immune surveillance keeps the malignant cell type at a non-fatal level. Thus, a fatal growth can initiate if $y_1$ acquires a mutation making the mutated cells $y_2$ less susceptible to the immune response.

In Figure 3a the MSC line $y_{1}$ does not gain resistance and the patient may live with a coexisting steady state for the healthy and malignant cells. However, in Figure 3b the non-fatal malignant cell type $y_{1}$ gives rise to a resistant cell line $y_2$ which initializes a fatal growth. This behavior resemble that described by Shlush et al. in 2017 ^[45] and discussed by others in ^[46,47,48] where the basic idea is that some malignant stem cells may survive treatment either due to previously developed resistance or for other reasons thus causing a fatal relapse for some AML patients. A sustained remission is seen in less than 35% of the young adults and 15% of older patients ^[48].

3.3.   Micro-environment causes oscillations in a heterogeneous clonal landscape

The stem cell pool can be highly heterogeneous ^[3] and we will investigate such a scenario where circular inhibiting effects between the MSCs give rise to oscillations and possible complications of such oscillations. At least three clones additional to the WSCs are needed to account for various observed dynamic complexities. Oscillating cell numbers could possibly be explained by having three or more malignant stem cells that inhibit each other.

For the case $n = 3$ we could have that $y_{1}$ is highly inhibited by $y_{2}$ but only slightly inhibited by $y_{3}$. In addition $y_{2}$ is not highly affected by the presence of $y_{1}$ but is highly affected by the presence of $y_{3}$ and lastly $y_{3}$ is highly affected by $y_{1}$ but not much by $y_{2}$. This specific form of inhibiting effect is illustrated in Figure 4 and the corresponding trajectories are seen in Figure 5.

Figure 6a explores a similar scenario but where the WSCs are affected differently by the three types of malignant cells, i.e., the WSCs are most inhibited by the presence of $y_{3}$ and least inhibited by $y_{1}$. Hence, if the virtual patient receives targeted treatment then the outcome depends on the targeted clone.

3.4.   Outcome of single cell type - targeted treatment depends on heterogeneous composition

The outcome of applying treatment in Figure 6a depends severely on which of the malignant cells the intervention addresses. Targeted treatment of $y_{1}$, $y_{2}$ and $y_{3}$ is depicted in Figure 6b-6d, respectively.

Figure 6b shows treatment targeting $y_{1}$ which inhibits the WSCs the least of the existing MSCs. The suppression of $y_1$ may unintentionally improve the growth conditions for the more agammaessive $y_{3}$ and initiate a fatal growth. Figure 6c shows that treatment targeting $y_{2}$ may increase the amount of WSCs as the treatment of $y_{2}$ improves the growth conditions for the less agammaessive $y_{1}$ which inhibits the WSCs the least. Figure 6d shows that treatment targeting $y_{3}$ results in a somewhat neutral outcome.

3.5.   Hydroxyurea treatment results in oscillatory cell counts

The aim of the simulated study presented in Figure 7 is to model a scenario observed for patients treated with hydroxyurea, namely a response having a fast drop followed by a rebound with oscillations as seen in the clinic.

Initially a fatal growth of a malignant cell type $y_{1}$ starts increasing the total amount of blood cells rapidly. This increase in cell counts triggers symptoms such as thrombosis and a targeted treatment for $y_1$ is initiated. The targeted treatment successfully reduces the amount of $y_{1}$ and after the treatment, $y_{1}$ enters a cyclic state with the two less agammaessive malignant types $y_{2}$ and $y_{3}$. The treatment initiates a drop in the total amount of cells. After some time the total amount of cells increases and stabilizes with oscillations at a non-fatal level.

3.6.   Oscillations observed in clinical data - patient specific model calibration

As mentioned, the MCC model is not mutation specific, such as JAK2, and can model the faster developing CML by changing the parameter values governing the micro-environment as an example. Observed oscillations in white blood cells for CML and platelets for a JAK2$^+$ MPN, are depicted in Figure 8, after treatment onset with hydroxyurea. The data are from ^{[22,23,24,25,31]}, respectively. The oscillations differ in period for the patients and the period of the oscillations in Figure 8 are reported to be 35-40 days, 53-69 days, 37-83 days, and 244 days for CML, respectively, while it is 27 days for MPNs.

To compare model output with data the outputs have been normalized to achieve the means and the amplitudes of the data. The linear drift seen in Figure 8b is artificial imposed on the model output to make the periods comparable. However, the deterministic model struggles to reproduce the irregularities of the measured oscillations, e.g., those in Figure 8a. For the data fits we have chosen $M_{i, j} = 0$ corresponding to ignoring additional mutations during the relatively short time spans considered.

3.7.   Inhibiting and stimulating properties generate severe clonal evolution

A common understanding is that an agammaessive malignant cell type arises from a sequence of mutations where the fitness of the malignant type increases with the number of mutations. We will explore both a deterministic and stochastic scenario depicted in Figures 9a and 9b, respectively. In the deterministic scenario the fitness of the malignant daughter clone increases in each mutation by an increase in the self-renewal rate and an inhibiting effect on the mother cell line. The result is a step-wise reduction in the amount of WSCs as depicted in Figure 9a.

In the stochastic scenario we have two cell lines present $y_0$ and $y_1$ at time $t = 0$. During the simulation they are both allowed to mutate. Each time a mutation happens, the daughter cells will have the parameters of the mother cell but inhibiting effect ($C$), self-renewal rate ($R$) and T-cell specific death rate ($K$) are sampled from log-normal distributions with the parameter values of the mother cells as mean values. Hence, the daughter cells have the possibility of becoming either more agammaessive or less fit than the mother cells. The new clone will have the possibility to mutate itself and thereby enabling the simulation to have a mixture of the sequential and the parallel aspects. The outcomes of the deterministic and the stochastic approach are analogues despite the variability in the stochastic approach and they illustrate how the healthy cells are step-wise out-competed due to the clonal evolution.

4.   Sensitivity, robustness and correlations

A sensitivity analysis of the parameters used for Figure 5 is provided in Supplementary section S.2. The sensitivity was performed for the non-default parameter values $K$ and $C$ with a $10\%$ change. Perturbing the parameter values by $10\%$ did not change the qualitative behavior of the system, i.e., oscillations were still observed for the perturbed parameter values. However, the amplitude and frequency of the oscillations did change as a consequence of the perturbation of the parameters as illustrated in Figures S.3 and S.4 and quantified in Table 1.

Moreover, the sensitivity analysis revealed an inverse correlation between $K$ and $C$ such that increasing $K$ by $10\%$ and lowering $C$ by 10% has almost identical simulation outcome. These observations align with the model understanding as a decrease in $K$ and an increase in the cells inhibitory effect will promote growth, thus giving larger oscillations whereas an increase in $K$ and a decrease in the cells inhibitory effects will decrease the growth potential thus making the amplitude smaller.

5.   Discussions

The presented work builds on previous models by allowing for the emergence of multiple populations of wild-type stem cells to mutated stem cells, which interact with each other as well as indirectly via the stem cell niche and the innate and adaptive immune response. A simulation study is described, exploring model dynamics under a minimal number of different parameter combinations focusing on the niche feedback inhibitions and the strength of the adaptive immune response. The different in silico scenarios encompass competition between clones of haematopoietic stem cells, wild-type ones and malignant ones, including during treatment with hydroxyurea, in a heterogeneous mixture of multiple clones. Findings are related back to various clinical and experimental observations. The model predictions show agreement between model predictions and observations of white blood cell and platelet counts, JAK2 allele burden, and cytokine counts in patients with haematologic malignancies. This is in favor of the validity of the model.

5.1.   Elimination, equilibrium and escape

The most common telling is that if a driver mutation happens then the malignant cell type increases sigmoidally while the healthy cell type decreases in concert as illustrated in Figure 2a. The healthy cell type does not need to go extinct as the two cell types may coexist and they may even coexist at a non-fatal level for the virtual patient, see Figure 2a. The healthy cell type may also coexist with a heterogeneous cancer, i.e., multiple malignant clones, as illustrated in Figure 5.

If the healthy cell line has entered such a coexisting state a treatment of one malignant cell type may leave room for the most fit of the remaining clones to increase vitally. We consider three clones, a healthy cell type and two malignant cell types, the second more agammaessive than the first and examine the result of treatment which leads to an increase in the T-cell related death rate of the malignant clone. In case of no treatment the first malignant clone may raise to a level causing complications, e.g., risk of thrombosis, but without eradicating the normal cell type while the second malignant clone persists at an insignificant level (Figure 2a). If instead the first malignant clone is suppressed by treatment, the second and more agammaessive malignant clone may bloom if it is insensitive to the treatment, thereby suppressing the healthy cells (Figure 2b). Hence, the treatment initiates a fatal development for the patient, who may had preferred to live with the risk of thrombosis.

In earlier studies, it is found that coexistence of clones being possible only for equal parameter values for the self-renewal rate of stem cells ($R$) ^[28]. We disprove this since we have shown that co-existence can occur for different values of $R$, as illustrated in Figure 2a. Moreover, the simulation shows that the clone with highest self-renewal rate ($R$) does not need to become the dominating clone if an earlier arrived clone has a sufficient suppressing niche feedback effect on the first one (high $C$ value).

5.2.   Resistance

A similar outcome results if the clone becomes resistant or less sensitive to the immune surveillance over time, illustrated by Figure 3. Notice, the developed resistance or decreased sensitivity to the immune response will cause the agammaegated cell count to increase toward a fatal level, as illustrated by the Brachiosaurus shaped development shown in Figure 3. Thus, the model explains the three Es of cancer immunoediting, elimination, equilibrium, and escape phases. First cancer grows exponentially but then the immune system becomes activated and initiates a defense, leading to stagnation of the cancer burden at a plateau. After a while either some of the malignant cells become resistant to the immune response or an additional mutation occurs resulting in a less immune-sensitive clone emerge. In both cases the immuno-insensitive clone starts to grow exponentially which is interpreted as the cancer escaping the immune response ^[41].

5.3.   Oscillatory cell counts caused by micro-environment

Medical action to changes in levels of clones may also be subject to caution. Consider three clones in addition to the normal clone as is the common idea of CHIPs where the number of cells of each clone may oscillate as illustrated in Figure 5. Assume that the three malignant clones have the same parameter values (but a little higher self-renewal than the normal clone to suppress it slightly) except for their mutually inhibitory effects; malignant clone 2 possesses a raised inhibition of malignant clone 1, malignant clone 3 has a raised inhibition of malignant clone 2, and malignant clone 1 imposes a raised inhibition of malignant clone 3 compared to the rest of the inhibitions (circular inhibition, see conceptual diagram in Figure 4). The number of cells of each malignant clone may oscillate, as illustrated in Figure 5. If in addition, some of the malignant clones have increased inhibition onto the normal clone, then oscillations in the normal cell counts will be more pronounced.

Figure 5 shows such an example, where the third malignant clone has a slightly increased inhibition of the normal clone. Depending on the values of the parameters, the oscillations and their mean may be more or less pronounced as illustrated in Figure 6. The period of the normal clone oscillating and its amplitude comes from the dominating malignant clone period and amplitude. If the patients blood is screened for specific mutations using NGS some of the clones may have common mutations. Thus one, two, or three of the malignant clones in said example may be registered as the same mutation, e.g., the different clones may all have the JAK2 mutation and thus only differ due to apparently 'insignificant' differences. In Figure 6a the sum of the clones is shown as an example and we may imagine that all three mutations may be classified as the JAK2 mutation. In a clinical setting such measurement may be taken rarely, e.g., each third month, every sixth month, yearly or not at all. To reliably observe periodic oscillations at least 3-5 measurements during such period are necessary otherwise measurements will be considered as 'random' fluctuations.

For the three species Lotka-Volterra equation (generalized logistic growth) global stable steady states of one, two or three (co-existing) species may occur. However, periodic solutions ^[49,50,51] or even heteroclinic orbits may exist ^[52,53] and for five species chaotic solutions exist ^[54]. Thus, oscillation dynamic of the current model is not surprising. However, another possible explanation for the observed oscillations could be that it is due to time delay.

For simplicity consider a single-type malignant stem cell pool. From Murray ^[55] we have that the steady state of the delay logistic equation with delay $T$ is locally asymptotic stable if $0 \leq r_yT < \frac{\pi}{2}$ and unstable if $r_y T > \frac{\pi}{2}$. Murray proved that for $r_y T = 2.1$ ($> \frac{\pi}{2}$) a periodic solution with period $P \approx 4.54 T$ exist. If the hydroxyurea induced oscillation with period $P \approx 100$ days is observed, corresponding to a time delay $T \approx 25$ days, then $r_y < \frac{\pi}{50} \approx 0.064$ per day holds in the case of MPNs, CML and in AML. Hence, such time delay would hardly cause oscillatory behavior.

5.4.   Targeted treatment

Continuing with the complex scenario of four clones of which three are malignant and one of the malignant clones possess a slightly elevated inhibitory effect on the normal clone, we return to the treatment discussion. We consider a hypothetical treatment which reduces the self-renewal rate for one clone only, yielding enormous differences in treatment outcome, depending on which clone the treatment targets. In Figures 6b, c, and 6d the self-renewal of clone 1, 2, and 3 is reduced at year 30. If malignant clone 1 is affected by the treatment, the malignant clone 3 expands resulting in a fatal eradication of the normal clone. If the malignant clone 2 is suppressed by treatment then the malignant clone 1 increases slightly but the normal clone almost recover fully resulting in a bearable coexisting steady state. If malignant clone 3 is reduced by the treatment, the outcome is neutral as the oscillations stops but the allele burden is unchanged for the patient. This coexisting steady state does have less functional cells and relatively many malignant cells. By changing parameter values the level of cells in the steady state may be adjusted. The relatively high number of malignant cells increase the probability of getting more agammaessive mutations. For a given probability (it is a Poisson process) the numbers of such are proportional to the number of mother cells. If the number of malignant cells, for some reasons is elevated by a factor $k$ then the number of more agammaessive additionally mutation per time increases by the same factor ^[56]. A reduction in the number of normal cells may cause reduced transport of oxygen and carbon dioxide by red blood cells throughout the body.

We emphasize that the model accounts for the apparently paradoxical effect of treatment, namely, how eradication of a clone may stimulate the growth of another, more fatal, one.

5.5.   Calibration to data

In clinical practice roughly half of all MPN-patients receiving treatment are characterized as good responders and the rest as poor responders when ignoring those taken off the treatment due to adverse events. ^[56] This phenomenon of good versus poor responders is not well understood but the multiple clone approach gives a plausible hypothesis. The multiple discoveries of visible oscillations in cell numbers over time support the hypothesis: Cyclic oscillations with period 35-40 days in leukocyte counts in CML during hydroxyurea therapy have been reported ^[22], similar oscillations with a period of 53-69 days were observed in platelets, neutrophils, eosinophils, basophils and reticulocytes too and continuing for a period after anti-leukemic therapy in ^[23]. In ^[24] oscillations with a period of 37-83 days were estimated, whereas ^[25] found oscillations with period of 244 days in platelets and white blood cells in a patient diagnosed with CML. All of these periods can be explained by the model, noticing that the time-scale results from the physiological parameters, which differ from mutation to mutation and thus from diagnosis to diagnosis, and in ^[29,30,31] faster oscillations with period 27 days have been observed clinically for JAK2$^+$ MPNs. The MCC model accurately reproduces all these clinical observations as illustrated in Figure 8. Additional support for such oscillations comes from the fact that CHIPs have been observed in a number of studies ^{[32,33,37,38]}.

In the clinic we observe strong fluctuations in JAK2 allele burden data and leukocyte data from MPN-patients treated with hydroxyurea. These fluctuations are typically not pronounced until after one year after treatment onset. The typical transient is a steep decay in allele burden (the fraction of malignant mature cell count) and total cell count followed by relapse toward the levels before treatment. After such transient the JAK2 allele burden may partly decay, stagnate, or even increase but in a fluctuating way. The MCC model can reproduce such observations as illustrated in Figures 6 and 7. The model illustrates the current clinical experience and knowledge that treating a blood cancer may introduce oscillatory cell counts.

5.6.   Clonal evolution

The MCC model may also explain the paradigm of a mutational evolution landscape as a sequential cascade of mutations. First one not very harmful mutation (possible as a result of a sequence of mutations) capable of surviving occurs. Next, an additional mutation (or sequence of such) occurs resulting in a more competitive clone, which more or less out-competes the first clone. These mutated stem cells may mutate further resulting in more and more harmful clones. In each step the preceding clone is inhibited and the normal cell type become more and more inhibited by the resulting clone after each additional mutation as illustrated in Figure 9a, where for simplicity only the self-renewal is increased in each mutational step. In Figure 9b, a stochastic simulation of a scenario reveals the same overall outcome, namely that the WSCs are stepwise outcompeted by increasingly agammaessive mutations.

6.   Conclusion and future work

In conclusion, the MCC model contributes to a better understanding and explanation of several observations in patients with myeloid cancers, including the previously observed (before the interferon-era and treatment with tyrosine kinase inhibitors) slow oscillations in cell counts of time scale 2-3 months for some CML patients. Similarly, the MCC model may explain the increased risk of thrombosis in MPNs despite cytoreductive therapy by demonstrating fast oscillations in cell counts of time scale 1 month for some MPN patients. Importantly, in the context of "Tumor Immune Surveillance" the different phases of immunoediting, elimination, equilibrium and escape, are explained in the model as well.

The model also suggests an evolution of CHIP towards the development of overt MPNs and the development of MPNs in the biological continuum from the early cancer stages towards the advanced myelofibrosis stage and ultimately leukemic transformation. Indeed, all these scenarios, being likely driven by chronic inflammation and inflammatory cytokines with increasing genomic instability and formation of subclones over time are explained in the model. Thus, most observations are a consequence of the model independent of whether one assumes sequential, parallel, or stochastic mutations. Finally, the large class of phenomena explained by the model and its ability to reproduce clinical data serve as a strong support for the validity of the model.

As a perspective for future studies, targets for therapeutic control by combination of drugs might be added to the model and an optimal control strategy might then be proposed to qualitatively study the behavior of the population of clones, with an objective of eradication or containment.

A.   Default parameters

B.   Non-default simulation parameters