Abstract
Recent experimental evidence suggests that acute myeloid leukaemias may originate from multiple clones of malignant cells. Nevertheless, it is not known how the observed clones may differ with respect to cell properties, such as proliferation and selfrenewal. There are scarcely any data on how these cell properties change due to chemotherapy and relapse. We propose a new mathematical model to investigate the impact of cell properties on the multiclonal composition of leukaemias. Model results imply that enhanced selfrenewal may be a key mechanism in the clonal selection process. Simulations suggest that fast proliferating and highly selfrenewing cells dominate at primary diagnosis, while relapse following therapyinduced remission is triggered mostly by highly selfrenewing but slowly proliferating cells. Comparison of simulation results to patient data demonstrates that the proposed model is consistent with clinically observed dynamics based on a clonal selection process.
1. Introduction
Leukaemia is a clonal disease of the haematopoietic system leading to extensive expansion of malignant cells that are nonfunctional and cause impairment of blood cell formation. Recent experimental evidence indicates that the malignant cell population might be composed of multiple clones [1], maintained by cells with stemlike properties [2,3]. A clone consists of genetically identical stem and nonstem cells. Relapse of the disease after therapy is a common problem of leukaemias [1].
To understand better the origins of acute leukaemia relapses, a genetic interdependence between clones at diagnosis and relapse has been investigated using gene sequencing and other techniques. In most cases of acute lymphoblastic leukaemia (ALL), the clones dominating relapse were already present at diagnosis but were undetectable by routine methods [4–6]. Owing to quiescence, very slow cycling or other intrinsic mechanisms [5,6], these clones survive chemotherapy and eventually expand [5,6]. This implies that the main mechanism of relapse in ALL is based on a selection of existing clones and not an acquisition of therapyspecific mutations [5]. Similar mechanisms have been described for acute myeloid leukaemia (AML), where clones at relapse are genetically closely related to clones at primary diagnosis [1,7] and did not have to acquire additional mutations during the course of disease [8,9].
Based on these findings, the evolution of malignant neoplasms can be interpreted as a selection process [10–12] of cells with properties that enable them to survive treatment and to expand efficiently. Cells with different mutations may have different growth properties [1]. Chemotherapy significantly alters growth conditions of cells and therefore, it may have a strong impact on the selection process. If cells dominating at diagnosis are sensitive to therapy, minor clones with intrinsic resistance [5,6,13] may expand more efficiently once the competing clones are eliminated by the treatment. The latter could explain manifestation of different cell clones at diagnosis and at relapse without a need for additional mutations in between.
The mechanism of the underlying selection process and its impacts on the disease dynamics and on the response of cancer cells to chemotherapy are not understood. Gene sequencing studies allow the genetic relationship between different clones to be deciphered, nevertheless the impact of many detected mutations on cell behaviour remains unclear [1], and often passenger mutations cannot be distinguished from relevant genetic changes [4]. Many authors, e.g. [14,15], have provided evidence for the heterogeneity of leukaemia stem cells (LSCs) attempting to identify LSC characteristics, for review see [16]. This heterogeneity is further supported by the results of gene sequencing studies [1,17,18]. The multifactorial nature of the underlying processes severely limits the intuitive interpretation of experimental data. Mathematical modelling is a powerful technique to close this gap and to provide quantitative insights into cell kinetics, fate determination and development of cell populations. It allows a systematic study of processes not yet accessible by experimental procedures. Mathematical models have been widely applied to analyse the regulatory mechanisms controlling the haematopoietic system and its diseases: for review see [19–22] and references therein.
The aim of this work is to investigate the impact of cell growth properties on the clonal selection process in acute leukaemias before and after treatment. We introduce mathematical models of the dynamics of leukaemia, which are extended versions of the models proposed earlier by our group [23–25]. The novel ingredients of the models in this work are: (i) heterogeneity and multiclonal structure of LSCs, (ii) different plausible feedback mechanisms and (iii) effects of chemotherapy.
As the mechanism of interaction between healthy and leukaemic cell lines is not well identified, we propose two models (figure 1). In the first one, we assume that leukaemic cells depend on haematopoietic growth factors and interact with haematopoietic cells via competition for these factors. The second model is based on the assumption that autonomous leukaemic clones compete with haematopoietic cells for niches in bone marrow, which leads to increased cell death because of overcrowding. The latter is supported by experimental findings showing signalindependent activation of important cell functions [26–28] and by an increased cell degradation observed in leukaemic patients [29–31]. Such interactions have not been considered in previous models.
The models proposed in this paper do not account for new mutations. Motivated by the experimental findings described above [5,6,8,9], we rather aim to understand which aspects of the dynamics of leukaemias can be explained by a selection process alone. It is interesting that expansion of a clone at relapse that could not be detected at diagnosis owing to limited sensitivity of methods can be misinterpreted as the occurrence of mutations [5]. This scenario seems to be relevant in the case of acute leukaemias with a shortduration treatment administration. Many acute leukaemias are genetically relatively stable in comparison to other cancers [32,33]. For this reason, on average many replications are necessary to acquire a new mutation. Consequently, it is less probable that cells acquire mutations during short treatment and, therefore, intrinsic resistance to therapy may be important, as suggested by available evidence [5]. The latter does not hold true for longterm drug administration, such as imatinib treatment in the case of chronic leukaemias. For this reason, our work focuses on the acute leukaemias. For completeness of this work and to check how mutations might influence the model dynamics as it concerns results presented in this paper, we have developed a version of the model with mutations. The simulations of the model confirm our conclusions for the model with mutations. The model and simulations are presented in appendix A.
Using mathematical models, we aim to identify which cell properties are compatible with intrinsic resistance to therapy and efficient expansion after treatment, and to compare them with the cell properties selected for before treatment. We perform computer simulations describing evolution of a multiclonal population of leukaemic cells during disease development and the contribution of different clones to the entire cancer cell population at different time points. The heterogeneity of the system is given by a certain number of leukaemic clones already present at the beginning of our observation. The models provide information on the influence of cell properties on the growth dynamics of the different clones in the presence or in the absence of chemotherapy. This allows us to understand how the cell properties selected for before treatment differ from those selected for during and after treatment and how treatment could be optimized to reduce relapses. Finally, we compare qualitatively model simulations to patients' data from clinical routine to show that the proposed models are consistent with clinical observations concerning the response to therapy and the time intervals between relapses of the disease. Details of mathematical formulation and parametrization of the proposed models are presented in appendix A.
2. Material and methods
2.1. Mathematical models
2.1.1. Model assumptions
The models used in this study are based on the models of healthy haematopoiesis proposed and analysed in [23,34–36] and extended to account for evolution of a single leukaemic clone in [25].
Based on the classical understanding of haematopoiesis [37], we assume that the system consists of an ordered sequence of different maturation states, socalled compartments. To describe time evolution of cell populations, we apply ordinary differential equations. The enormous amount of cells forming the haematopoietic system justifies this approach [37,38]. Evolution of small cell population in the posttherapy period is modelled by cutting off the initial data which are below a minimal threshold, as it was, for example, proposed in [39].
We model time dynamics of one healthy cell lineage and an arbitrary number of leukaemic clones. In the description of cell differentiation within each cell line, we choose a twocompartment version of the multicompartment system established in [23]. The model focuses on the maintenance of primitive cells and differentiation from undifferentiated, proliferating cells to differentiated, postmitotic cells. In the case of healthy haematopoiesis, the proliferating cells are haematopoietic stem cells (HSCs), haematopoietic progenitor cells (HPCs) and precursor cells, the postmitotic cells are mature cells, e.g. white cells. The twocompartment architecture is based on a simplified description of the multistages differentiation process. Nevertheless, as shown in [24,35,36] models consisting of two compartments capture the desired dynamics of the multicompartmental cell population. This reduces the complexity of the differentiation process to focus on mechanisms and effects of competition between different cell lines.
Each proliferating cell type is characterized by the following cell properties:

— proliferation rate, describing how often a cell divides per unit of time;

— fraction of selfrenewal, describing the fraction of daughter cells returning to the compartment occupied by the mother cells that gave rise to them. Based on our earlier work and on compatibility with clinical data [23], we assume that the fraction of selfrenewal of haematopoietic cells is regulated by feedback signalling;

— death rate, describing what fraction of cells dies per unit of time. For simplicity, we assume that under healthy conditions proliferating cells do not die and postmitotic mature blood cells die at a constant rate. We assume the same for leukaemic cells in Model 1, while in Model 2 (see below), we consider, additionally, cell densitydependent death rates for all bone marrow cell types if the marrow space is overcrowded. The considered marrow cell types include immature haematopoietic cells and mitotic and postmitotic leukaemic cells. Overcrowding is defined when marrow cell counts exceed the steadystate marrow cell count by two to three times. In this case, the death rate of postmitotic leukaemic cells consists of their intrinsic death rate and the death triggered by spatial competition.
Production of healthy blood cells is regulated by a negative feedback [40–42], mediated by cytokines, such as GCSF or EPO [37,42,43]. If there is a shortage of blood cells of a certain type, the concentration of signalling molecules increases and stimulates expansion of precursor cells. This effect is modelled using a negative feedback loop as proposed in [23]. Analysis and simulation of the model of healthy haematopoiesis, validated based on the clinical observations after stem cell transplantations [23,44,45], indicate that the regulation of the selfrenewal is a more efficient mechanism than the regulation of the proliferation rates. Similar conclusions have been drawn using the models of multistage cell lineages applied to regeneration and maintenance of the mouse olfactory epithelium [46,47]. Therefore, in the remainder of this paper we assume that the regulatory mechanism is based on the feedback inhibition of selfrenewal depending on the level of mature cells.
2.1.2. Model of the healthy cell line
We denote by p^{c} the proliferation rate of mitotic haematopoietic cells and by a^{c} the corresponding fraction of selfrenewal. The death rate of mature blood cells is denoted by We denote the concentration of healthy cell types at time t by c_{1}(t), c_{2}(t), corresponding to mitotic and mature cells, respectively. The flux to mitosis at time t equals p^{c}(t)c_{1}(t). During mitosis, a mother cell is replaced by two daughter cells. The outflux from mitosis at time t equals 2p^{c}(t)c_{1}(t), of which the fraction 2a^{c}(t)p^{c}(t)c_{1}(t) stays in compartment 1 (process referred to as selfrenewal). The fraction moves to compartment 2 (process referred to as differentiation).
We denote the value of the feedback signal at time t by s(t), which takes values between zero and one. Selfrenewal of a certain cell type at time t is assumed to be given as a maximal possible selfrenewal of this cell type multiplied by s(t). Following [23,25], we chose which can be derived from cytokine kinetics [23]. The constant k^{c} depends on the rate of extrahaematopoietic cytokine degradation by liver or kidney and on the rate of cytokine degradation by haematopoietic cells. The latter depends on the densities of cytokine receptors on haematopoietic cells [45].
We obtain the following system of ordinary differential equations, where corresponds to the maximal possible selfrenewal of HSCs. 2.1 2.2 2.3
The two different models proposed in this paper differ with respect to the interaction of leukaemic and haematopoietic cells. We consider two cases. In Model 1, leukaemic cells depend fully on haematopoietic cytokines, whereas in Model 2 they are totally independent of environmental signalling. In this sense, Models 1 and 2 can be understood as the two opposite extremes of a continuum. In reality, both mechanisms, competition for environmental signals and direct inhibition or death of haematopoietic cells, may contribute to impaired haematopoietic function [48]. A schematic of the model is given in figure 1.
2.1.3. Model 1
We assume that leukaemic cells depend on the same feedback signal as their healthy counterparts and that the postmitotic leukaemic cells (blasts) decrease the supply of the factor. It describes a competition between healthy and leukaemic cells for survival signals, which results in downregulation of selfrenewal. A schematic of the model is given in figure 1.
To write the corresponding equations, we denote the number of leukaemic clones by n. As for the haematopoietic cells we consider mitotic and postmitotic cell compartments for each leukaemic clone. Let denote the proliferation rate of mitotic cells in leukaemic clone i and the corresponding maximal fraction of selfrenewal. By we denote the clearance rate of postmitotic cells of clone i. Denote by the level of mitotic cells of clone i and by the level of postmitotic cells at time t. These assumptions result in the following system of ordinary differential equations: 2.4 2.5 2.6 2.7 2.8 2.9 2.10 2.11
The expression for s(t) is a special case of where we assume that for all i. This simplification corresponds to the observation that the density of cytokine receptors is similar on cells of all leukaemic clones. For the major cytokine of the myeloid line, GCSF [41], this is true for many patients [49]. As there is evidence that in some patients receptor densities may differ between different leukaemic clones [49], we have repeated all simulations with a randomly chosen value for each clone, ranging from 30 to 100% of k^{c}. This heterogeneity had no significant impact on the model results. As in many cases the receptor density on leukaemic cells is of the same order of magnitude as that on haematopoietic cells [49,50], we assume also for the simulation of patient examples.
2.1.4. Model 2
There is evidence that in some leukaemias malignant cells show constitutive activation of certain signalling cascades and thus may become independent of external signals [26–28]. We consider this scenario in Model 2. In contrast to Model 1, we assume that leukaemic cells are independent of haematopoietic cytokines, whereas the haematopoietic cell types depend on the nonlinear feedback described earlier. Interaction between the healthy and cancer cell lines is modelled through a competition for space resulting in an increased cellular degradation, for example, owing to overcrowded bone marrow space. This is consistent with the observation of an increase of markers for cell death such as LDH [29–31]. Several mechanisms underlying this spatial competition have been proposed: (i) physical stress owing to overcrowding leads to extinction of cells (e.g. [51]; recently challenged by [52]), (ii) competition for a limited niche surface expressing certain receptors (contact molecules) necessary for survival of the cells [53,54] and apoptosis if no contacts to these molecules can be established [55].
We model the space competition by introducing a death rate that increases with the number of cells in bone marrow and acts on all cell types residing in bone marrow, i.e. mitotic and postmitotic leukaemic cells as well as mitotic haematopoietic cells. For simplicity, we assume in Model 2 that all leukaemic cells stay in bone marrow, as the number of leukaemic cells exiting bone marrow is highly variable among individuals and only partially dependent on the leukaemia subtype [56–58] and as it is not well understood which mechanisms are responsible for marrow egress and high interindividual variability. The presented results are robust with respect to this assumption: we repeated all simulations for the cases that 10, 50 or 90% of the most mature leukaemic blasts exit bone marrow. This has an impact on the time dynamics of marrow blast count but does not influence the cell properties that are selected for.
Let d(x) be an increasing function with This function describes the death rates of bone marrow cells in dependence of bone marrow cell counts x. We assume that under healthy conditions there exists no cell death owing to overcrowding. Enhanced cell death can be observed only if total bone marrow cellularity increases beyond the threshold level. This assumption is in line with bone marrow histology [59]. Therefore, we assume that d(x) = 0 for where is the steadystate count of mitotic healthy cells.
Assuming that the haematopoietic cell lineage is regulated as described above, we obtain the following system of differential equations: 2.12 2.13 2.14 2.15 2.16 2.17 2.18 2.19 2.20Here denotes the mitotic cells of clone i, their fraction of selfrenewal and their proliferation rate. The level of postmitotic cells of clone i is denoted as . In the absence of marrow overcrowding, these cells die at rate .
2.1.5. Chemotherapy
We focus on classical cytotoxic therapy acting on fast dividing cells, which is introduced to the models by adding a death rate proportional to the proliferation rate. The assumption is motivated by the fact that many of the classical therapeutic agents used for the treatment of leukaemias act on cells in the phase of division or DNA replication [60]. Therefore, the rate of induced cell death is proportional to the number of cycling cells. We assume that the linear factor, denoted by k_{chemo}, is identical for all mitotic cells. Under chemotherapy, the equation for mitotic haematopoietic cells in Model 1 takes the form 2.21
Similarly, we obtain for mitotic cells of leukaemic clone i 2.22
Chemotherapy in Model 2 is introduced analogously.
2.2. Simulations
We perform numerical simulations of the models to investigate which leukaemic cell properties lead to survival advantage during evolution of leukaemogenesis and recurrence under chemotherapy. As explained before, the models do not account for additional mutations taking place during the therapy. Instead, we investigate evolution of a certain number of leukaemic clones present at a starting time point. We assume that in healthy individuals the haematopoietic cells are in a dynamic equilibrium, i.e. production of each cell type equals its clearance. Initial conditions for the computer simulations are equilibrium cell counts in the haematopoietic cell lineage and a small cell number for different leukaemic clones. We assume that the initial number of leukaemic clones in each patient is 50. This number is arbitrarily chosen. All presented simulations were repeated for different numbers of leukaemic clones (between 3 and 100), which led to comparable results (see the electronic supplementary material, figures S1 and S2). We assume that primary diagnosis and diagnosis of relapse occur when healthy blood cell counts are decreased by 50% of their steadystate value. We perform simulations for 50 patients, i.e. 50 different sets of initial data and model parameters, with 50 leukaemic clones per patient. The growth properties of the leukaemic clones are chosen randomly within certain ranges. The choice of model parameters is described in appendix A. The simulations follow the following algorithm:
(i) We start from healthy equilibrium in the haematopoietic lineage and one mitotic cell per kilogram of body weight for each leukaemic clone and run simulations until the number of healthy mature blood cells decreases by 50%. We investigate properties of the clones with the highest contribution to the total leukaemic cell mass. The clones under consideration are those which together constitute 80% of the total leukaemic cell mass. In the following, we denote these clones as ‘significantly contributing clones'. This procedure is taken to reflect the sensitivity of the detection methods. In more than 90% of the patients, two to four clones sum up to more than 95% of the total leukaemic cell mass. Taking a threshold between 80 and 95% to define ‘significantly contributing clones' has little influence on the result. Furthermore, more than 97% of the clones that are considered as insignificant by this method consist of less than 1% of the leukaemic cell mass. This number is in agreement with the detection efficiency reported in the literature [61].
(ii) Next, we simulate chemotherapy. For simplification, we consider seven applications of cytotoxic drugs (one per day during 7 following days, corresponding to standard inductions). Simulations show that the number of drug applications has no influence on the presented qualitative results. As proposed in the literature [39], we assume that a cell population has become extinct if it consists of less than one cell. Initial conditions for the posttherapy period are obtained from cell counts after therapy where counts of extinct populations are set to zero. We continue simulations until mature blood cell counts decrease by 50%, and then assess the cell properties of the clones contributing to relapse.
Calibration of the haematopoietic part of the model to clinical data and parameters for simulation of two patient examples can be found in appendix A. As in clinical routine only few key mutations are monitored, we choose patient examples with different key mutations detected at diagnosis and at relapses. Such data are relatively rare, therefore we focus on two patients. Simulations are performed using standard ODEsolvers from the Matlabsoftware package (v. 7.8, The MathWorks, Inc., Natick, MA, USA) which are based on Runge–Kutta schemes.
3. Results
3.1. Clonality at diagnosis
We solve the models numerically to obtain insight into the contribution of different leukaemic clones to the total leukaemic cell mass. Simulations indicate that at the diagnosis rarely more than three to four clones significantly contribute to the total leukaemic cell mass. In most cases, more than 40–50% of the total leukaemic cell mass originates from a single leukaemic clone. This finding is identical for both considered models.
3.2. Properties of clones at diagnosis
Simulations indicate that the clones significantly contributing to the leukaemic cell mass have high proliferation rates and high selfrenewal potential (high fraction of symmetric selfrenewing divisions). Such configuration of parameters leads to an efficient cell expansion. The properties of clones contributing significantly to leukaemic cell mass at diagnosis are depicted in figure 2. This finding is identical for both considered models.
3.3. Clonality at relapse
The clonality at relapse is comparable to the clonality at diagnosis. Rarely more than three clones significantly contribute to the total leukaemic cell mass. This finding is the same for both considered models.
3.4. Properties of clones at relapse
The properties of the leukaemic clones responsible for relapse depend on the efficiency of chemotherapy. We run computer simulations for varied efficiency of chemotherapy, namely different death rates imposed on mitotic cell compartments. In the case of inefficient chemotherapy, i.e. killing rates of mitotic cells being relatively small, the clones present at diagnosis are also responsible for relapse. These clones have high proliferation rates and high selfrenewal potential. In the case of more efficient chemotherapy, i.e. killing rates of mitotic cells being higher, the clones responsible for primary presentation differ from the clones responsible for relapse. Compared to the clones leading to primary presentation, the clones responsible for relapse have low proliferation rates but high selfrenewal potential. The properties of clones contributing significantly to leukaemic cell mass at diagnosis and at relapse are depicted in figure 3. Both models lead to similar results.
The result that slow cycling is an important selective mechanism and is compatible with the finding that cells in minimal residual disease samples are highly quiescent [6]. It is further supported by the fact that addition of anthracyclines, which act independent of cell cycle [62], leads to improved outcome of relapse therapies in ALL [63].
3.5. Treatment of relapse
If the same treatment strategy as in the case of primary treatment is applied to a relapsed patient, remission time is significantly shorter (figure 4). Second relapse is mostly triggered by the same clones as primary relapse. With repeated chemotherapy, clonal composition changes in favour of the clones with minimal proliferation (Clone 5 in figure 4). This finding is in agreement with data from clinical practice in ALL suggesting that the clones selected for at relapse possess inherently reduced sensitivity to treatment [5] and may be also responsible for second relapse [5]. The dynamics of leukaemic cells in our model are in good agreement with data from clinical practice: chemotherapy is able to reduce leukaemic cell load after relapses [4], nevertheless, this reduction does not lead to durable remission [63]. This reflects the worse prognosis of relapsed patients [13,63,64]. The increasing fraction of cells with reduced drug sensitivity predicted by the simulations explains the experimental finding that cells present at relapse are more resistant to chemotherapy than cells present at initial diagnosis [13,64]. It also shows that repetition of the same induction therapy leads to worse results in relapse compared with primary manifestation [63]. The selection of slowly cycling cells predicted by our model seems to be an important mechanism in AML. It was demonstrated that induction of cell cycling enhances chemosensitivity of leukaemic cells [65] and improves patient outcome after therapy [66]. Our model suggests that repeated chemotherapy can lead to the selection of clones that are not competitive in the natural environment, i.e. which can be outcompeted by clones sensitive to chemotherapy after cessation of the treatment.
3.6. Shortterm expansion efficiency does not correlate with longterm selfmaintenance
If leukaemic cell behaviour depends on haematopoietic cytokines (Model 1), the current signalling environment influences the expansion of leukaemic clones. In this scenario, it is possible that fast proliferating cells with low selfrenewal potential dominate the leukaemic cell mass during an initial phase. If, with increasing leukaemic cell mass, selfrenewal becomes downregulated, e.g. owing to occupation of bone marrow niche, eventually the cell clone with the highest affinity to selfrenewal survives, although its proliferation might be slow. An example of time evolution during an early phase is depicted in figure 5.
3.7. Late relapses can originate from clones that were already present at diagnosis
Simulations of Model 1 indicate that late relapses, e.g. relapses after more than 3 years, can originate from clones that were already present at diagnosis but did not significantly contribute to the leukaemic cell mass at that time. These relapses are triggered by very slow proliferating cells which survive chemotherapy and then slowly grow. At primary diagnosis, fast proliferating clones dominate. The slowly proliferating clones are then selected by chemotherapy. This finding is able to explain relapses without additional mutations occurring after primary diagnosis. Thus, temporary risk factor exposure (e.g. chemicals or radiation) can also be responsible for very late relapses and presentations.
3.8. Comparison of simulations to patient data
To check whether the proposed modelling framework is consistent with the observed dynamics of leukaemia, we calibrate the model to data of two patients with multiple relapses. The selected two patients showed different AMLtypical mutations. Properties of leukaemic cells and their impairment owing to chemotherapy cannot be measured directly and the effects of specific mutations on cell dynamics are not well understood. The available data include time periods between induction/consolidation chemotherapy and relapse as well as the percentages of leukaemic blasts in the bone marrow at diagnosis, followups and relapse. In addition, emergence and subsequent elimination of leukaemia driving mutations (FLT3, MLLPTD) in the bone marrow cells were precisely monitored using molecular biology methods [68–70]. We verify whether, and under which assumptions concerning the cell behaviour, the proposed model is compatible with clinical observations. This can serve as a qualitative ‘proof of principle’ and leads to hypotheses concerning changes in cell properties induced by the respective mutations. We assume that each mutation is associated with one leukaemic cell clone. We interpret differences at diagnosis and at relapse as the result of a clonal selection process owing to chemotherapy and cell properties. For this study, we apply Model 2, as simulations over a large range of parameters showed that remissions shorter than 150 days are only compatible with Model 2.
Simulations of the evolution of leukaemic clones in the two patients are depicted in figures 6 and 7. The results show that bone marrow blast fraction can be well described by the model. In Patient 1, FLT3ITD mutation of a length of 39 bp is detected at diagnosis. This mutation becomes extinct and the relapse is triggered by two different FLT3ITD mutations (42 and 63 bp). This behaviour is reproduced in the model simulation. At diagnosis, leukaemic cell mass is mainly derived from one clone while at relapse two different clones contribute to leukaemic cell mass.
In Patient 2, both FLT3ITD mutation and MLLPTD mutation were detected at diagnosis. The MLLPTD mutation practically did not contribute to relapse. The model reflects this scenario. At diagnosis, two different clones contribute to leukaemic cell mass, one of which becomes extinct and is not detected at the relapse. In this patient, the clone responsible for relapse behaves similar to the HSC lineage. Thus, classical cytotoxic treatment would not lead to its eradication. This is an indication for application of new antileukaemic drugs, if feasible, or for bone marrow transplantation.
4. Discussion
We have examined the impact of cell properties on clonal evolution in acute leukaemias during the course of disease. We have considered two different mathematical models, representing different modes of interactions between normal haematopoietic and leukaemic cells. In Model 1, leukaemic cells depend on haematopoietic cytokines, niches or other environmental factors. In Model 2, the leukaemic cells are independent of these aforementioned determinants and the only interaction between benign and malignant cells is owing to a competition for bone marrow space.
Model simulations suggest that clones with a high proliferation rate and a high selfrenewal are favoured at primary diagnosis. The results indicate that the number of clones significantly contributing to the leukaemic cell mass is relatively small, even if a large number of clones with different leukaemia driving mutations might coexist in the bone marrow. For example, in our simulations it was reduced from 50 to 2–5. This result is in agreement with data from recent gene sequencing studies and explains these data. In these studies [1,15], at most four contributing clones were detected in the case of AML and at most 10 in the case of ALL. In many patients, this number was even smaller. Our study implies that clonal selection owing to different growth characteristics is an efficient mechanism to reduce the number of clones contributing to leukaemic cell burden. Clones not contributing to primary disease manifestations might rest in a slowly proliferating or quiescent state and expand at relapse. Chemotherapy exerts a strong selective pressure on leukaemic clones, and thus has a considerable impact on the clonal composition during relapse.
In the case of insufficient chemotherapy, the relapse can be triggered by the same clones as the primary disease. In the case of more intensive therapy regimens, relapses are mostly triggered by different clones than primary disease. This has also been concluded from experimental studies [1]. Our models suggest that chemotherapy selects for slowly proliferating clones with high selfrenewal properties. Depending on efficiency of the therapy, it is also possible that clones with high proliferation and high selfrenewal potential are responsible for relapse.
In this study, we have focused on classical cytotoxic chemotherapy, mostly acting on mitotic cells. This explains the selection of slowly proliferating clones, among which those with high selfrenewal potential have a competitive advantage, as shown in earlier studies [23–25]. High proliferation rates constitute a disadvantageous factor under cytotoxic treatment, as fast proliferating cells are responsive to even moderately intensive therapy regimens. Relapses due to such clones are only possible if LSCs at the same time have a high selfrenewal potential, which is an advantageous factor for expansion and survival. Otherwise, they would be outcompeted by slowly proliferating cells with high selfrenewal. Fast proliferating cells with low selfrenewal have never been observed at relapse in our simulations. Their emergence at relapse could only be explained by additional mutations acquired after initial treatment. The selection of slowly proliferating cells may explain emergence of resistance in relapses. In such case, applying an identical therapeutic regimen to primary presentation and relapse has limited effects in the absence of new mutations.
The principle of clonal competition in leukaemia evolution and the fact that resistant subclones might be responsible for relapse have been discussed for a long time [5]. Using mathematical modelling, we have provided for the first time evidence that selfrenewal potential is a major force behind this mechanism and that cells responsible for a relapse show high selfrenewal in nearly all cases. This finding is new and cannot be concluded from biological data so far.
In appendix A, we study a model that includes occurrence of new mutations in addition to the selection process. In this scenario, the number of clones detectable at diagnosis and at relapse and their respective properties are practically identical to the scenario without mutations. This finding underlines that clonal selection has an important impact on the evolution of leukaemic cell properties.
The exact nature of interaction between leukaemic and haematopoietic cells is not well understood. Moreover, it is well known that leukaemias show high interindividual heterogeneity concerning symptoms and survival [67]. Therefore, it is possible that different mechanisms may be relevant in different cases. Simulation results suggest that the evolving cell properties are robust with respect to the assumptions on the exact mode of interaction between haematopoietic and leukaemic cells and are similar in different scenarios and different patients. Common features of both models are: (i) relapses can be explained by cells that were already present at diagnosis. (ii) Before therapy clonal evolution selects for cells with high proliferation rate and high selfrenewal. (iii) Cytotoxic treatment selects for cells with slow proliferation and high selfrenewal. Thus, it is possible to draw conclusions on leukaemic cell properties, even if their interaction with healthy haematopoiesis is not known in detail.
Nevertheless, the two proposed models exhibit some different dynamical properties, namely: (i) complete remissions lasting less than 150 days are only possible in Model 2. (ii) In Model 2, it is possible that leukaemic and nonleukaemic cells coexist at ratios compatible with sufficient haematopoiesis for long periods. (iii) In Model 1, clones can temporarily expand and then be outcompeted. In Model 2, the clone with the fastest expansion is dominant at all times until treatment. (iv) In Model 2, the leukaemic cell load can be reduced to a new steady state under chronic application of cytostatic drugs. In Model 1, expansion of leukaemic cells can be reduced in speed but eventually healthy haematopoiesis will be outcompeted. This may have application in the treatment of fast relapsing patients, as fast relapse can only be explained by Model 2.
Up to now, it cannot be decided which model is more realistic. For each of the models, there exist supportive findings. Model 1 is supported by observations on expression of growth factor receptors by leukaemic cells similar to those by haematopoietic cells [49,50], expansion of leukaemic cells in the presence of cytokines in some patients [71] and dependence of leukaemic cell selfrenewal and proliferation on chemokines needed for haematopoietic cell maintenance [72]. The facts supporting Model 2 are enhanced cell death in marrow samples [73] and increased markers for cell death/cell lysis in serum [29,74], independence of leukaemic cells from important environmental signalling cues in the presence of some mutations [26] and necessity of physical contacts to marrow stromal cells needed for cell survival [53–55]).
Our models support the hypothesis that processes of clonal selection are important mechanisms of leukaemia relapse, which can be responsible for expansion of different cell clones without a need for new mutations. A testable prediction of our models is that more sensitive methods should reveal larger numbers of different clones that exist but do not significantly contribute to the leukaemic cell mass. Another prediction is that cells present at relapse show mutations responsible for high selfrenewal.
Calibration of the models to patient data shows that the proposed framework is compatible with the observed clinical course in the considered two datasets. The predicted selection of slowly proliferating cells with high selfrenewal ability is consistent with clinical observations. Our results may have relevance for personalized medicine. Deep sequencing techniques might provide information on the genetic interdependence of the clones present at diagnosis and relapse [1]. Our model suggests that insufficient therapy may lead to the presence of the same clones at diagnosis and relapse. If the clones present at diagnosis and relapse are not identical but related, i.e. they share common somatic mutations [1], relapse may be due to a selection process. In this case, it is probable that the clones present at relapse show a slow proliferation and a high selfrenewal. One possible implication might be the application of cellcycle independent drugs, such as those used in targeted therapies.
Funding statement
This work was supported by the Collaborative Research Center, SFB 873 ‘Maintenance and Differentiation of Stem Cells in Development and Disease’. A.M.C. was supported by ERC Starting Grant Biostruct and EmmyNoetherProgramme of German Research Council (DFG).
Acknowledgements
The authors would like to thank Prof. Marek Kimmel for many helpful advices during preparation of the manuscript.
Appendix A
A.1. Calibration of the haematopoietic cell lineage
In the absence of leukaemic clones, both considered models reduce to the same model of the haematopoietic system. In steady state, this model has the following form: A1 A2 A3
Assume we know and . It holds A4 A5
Knowing a^{c}, we can calculate such that the steadystate population size c̄_{2} is satisfied. We calibrate the model to the data on production of neutrophil granulocytes, which constitute the majority of mature white blood cells (50–70%). Lymphocytopoiesis is a complicated process involving lymphatic organs and not only the bone marrow [75]. Therefore, we restrict ourselves to the myeloid line. It holds for the steadystate count of neutrophils [76], We interpret as the total amount of mitotic neutrophil precursors in bone marrow. Based on the data from [59], we assume that about 20% of bone marrow cells are mitotic precursors of neutrophils (the interindividual variations are considerable). The total bone marrow cellularity is about 10^{10} cells per kilogram of body weight [77]. Therefore, we take A6
Assuming an average blood volume of 6 l and an average body weight of 70 kg, we calculate a mature neutrophil count of A7
Neutrophils have a halflife in the bloodstream, T_{1/2}, of about 7 h [78]. From this, we calculate A8
We obtain A9i.e. about once per 1.5 days. We know from bone marrow transplantation [79] that patients need about 15 days to reconstitute to 5 × 10^{8} neutrophils per litre of blood (4 × 10^{7} kg^{–1}) after infusion of 5 × 10^{6} immature cells per kilogram of body weight. For A10this constraint is met.
A.2. Model parameters
For the haematopoietic branch, we chose parameters obtained from the calibration in §A.1. For simplicity, we assume k^{c} = k^{l} for the feedback mechanism in Model 1. We set the clearance rate of blasts (in the absence of effects of overcrowding) to . This is based on the apoptotic indices (fraction of dying cells) reported in the literature which are ≈0.19 ± 0.16 (19 ± 16%) [80,81]. Choosing blast clearance between 0.1 and 0.5 changes the speed of leukaemic cell accumulation but, as revealed by additional simulations, not the cell properties that are selected.
We chose In histological images of healthy adult bone marrow, a large part of the bone marrow cavity consists of fat and connective tissue and is free of haematopoietic cells. To reflect this fact, we set , where is the steadystate count of mitotic healthy cells. In the simulations, d_{const.} was set to 10^{−}^{10}. This choice implies that if bone marrow cell counts are three times higher than in the steady state, the additional death rate owing to overcrowding is of the order of magnitude of mature cell clearance. The results remain unchanged qualitatively, even if values of d_{const.} vary within different orders of magnitude.
We apply chemotherapy on 7 following days for 2 h. Different treatment intervals lead to comparable results. In the depicted simulations k_{chemo} has been set to values between 20 and 30 for less efficient therapy and to values between 40 and 60 for efficient chemotherapy. For different choices similar results are obtained. The higher k_{chemo}, the stronger the selection for high selfrenewal and slow proliferation and the lower the probability of relapse. The lower k_{chemo}, the higher the probability that clones contributing to primary presentation are among clones contributing to relapse. For k_{chemo} between 40 and 45, the obtained results are similar to the results obtained in experiments [1]. The parameters of the leukaemic clones are chosen randomly from uniform distributions, assuming that cells divide at most twice per day and that selfrenewal is between zero and one.
A.3. Calibration to patient examples
Both patients were treated within clinical trials at the University Hospital of Heidelberg after obtaining their written consent. Details on the patients' characteristics and therapeutic regimens can be found in the electronic supplementary material, table S1. Model parameters can be found in the electronic supplementary material, tables S2 and S3. For both patients, the presence of specific key mutations was assessed in clinical routine. We chose cases where mutations get lost because of treatment and new mutation are detected at relapse. We interpret this as the result of clonal evolution. As the two patients harbour different mutations, their leukaemic cells can have different properties. Chemotherapy is modelled by increasing death rates for mitotic cell types during the duration of each cycle. For simplicity, we did not model kinetics of single chemotherapeutic agents. Instead, the therapyinduced death rates are assumed to remain constant from the first to the last day of each treatment cycle. In pharmacology, the exposition to a drug is measured using the ‘area under the curve’ (AUC). This is the integral of concentration (or drug effect) over time [82]. The AUC in our case is k_{chemo} · Δt, where Δt is the period of drug action. The AUC over 1 day of therapy is similar for the single patient examples and the simulations in figure 3. Only myeloablative treatment before transplantation has a higher AUC. The presented results are based on Model 2. Model 1 is not compatible with remissions lasting less than 150 days. For simplicity, we count all leukaemic cell types as blasts.
A.4. Model with mutations
In the following, we describe an extension of Model 1 which includes mutations. There is an evidence that a preleukaemic HSC compartment serves as a reservoir of accumulated mutations [33]. This hypothesis is supported by the finding that some of the mutations recurrently observed in leukaemia already exist in the HSC compartment of a majority of leukaemia patients [33]. These preleukaemic HSCs seem to behave similar to normal HSC [33]. The hypothesis is that a relatively small number of additional hits may transform these preleukaemic HSCs into LSCs. Nevertheless, details of the underlying dynamics are not well understood [33].
We make the following assumptions:

— LSC with new properties can be generated either from preleukaemic HSC or from LSC owing to acquisition of mutations. This is in line with the current knowledge [7,32].

— For simplicity, we assume that the influx of new LSC from the preleukaemic compartment is constant in time. This assumption is made due to simplicity, because at the moment the dynamics of the preleukaemic compartment is not well understood [33]. We neglect mutations leading from normal HSC, i.e. nonpreleukaemic HSC to LSC.

— We assume that most mutations in LSC are acquired during replication of the genome and neglect other possible origins. In line with [7], we neglect mutations leading to dedifferentiation of nonLSC leukaemic cells.
Let be the level of LSC of clone i at time t. The flux to mitosis is then Out of mitosis, we obtain , where is the fraction of LSC selfrenewal of clone i at time t. We assume that the fraction ν of these cells is mutated, ν takes into account replication errors in relevant genes and is assumed to be constant. The influx α_{i}(t) of mutated LSCs due to new mutations occurring in clone i at time t is therefore .
We obtain the following set of equations describing dynamics of clone i: A similar system of equations has been obtained by Traulsen et al. [83]. As is considered to be postmitotic, we do not distinguish between cells that acquired a mutation during the divisions and those that did not.
The influx α(t) of mutated cells at time t is given by where γ is the constant influx from the preleukaemic compartment and N(t) the number of leukaemic clones present at time t.
We accept the rate α(t) as the rate of an inhomogeneous Poisson process. Poisson processes describe rare events [84,85], therefore they are a suitable tool to model mutations. It is known from probability theory that, if τ_{1} is a jump time of the inhomogeneous Poisson with rate λ(t), then the next jump time τ_{2} can be generated by solving the equation for τ_{2}. Here u is a uniformly distributed random variable [86]. We further know that if u is uniformly distributed in then −log(1 − u) is exponentially distributed with parameter 1 [85].
We simulate the system with mutations as follows: at time t_{0} = 0 we draw an exponentially distributed random number r_{1} with parameter 1. We simulate the system until the timepoint t_{1} which fulfils At timepoint t_{1}, a mutation occurs that gives rise to a new LSC. This is modelled by adding to the system a new LSC clone, consisting of one LSC and no less primitive leukaemic cells. We assume that the mutation occurs in a random gene position therefore we assign random cell properties to the new clone, i.e. selfrenewal and proliferation chosen randomly from uniform distributions (proliferation rate between 0.01 and 0.9, selfrenewal between 0.5 and 1). This choice is made for the sake of simplicity as details about the impact of mutations on cell behaviour and the underlying probability distributions are not known [1]. We then draw another random number r_{2} and continue simulations until timepoint t_{2} fulfilling etc. We start simulations from the equilibrium of the haematopoietic system and one LSC with random properties.
The results obtained from these simulations are similar to the results from the model without mutations. At primary diagnosis, cells show high selfrenewal and high proliferation while at relapse cells show high selfrenewal and reduced proliferation (electronic supplementary material, figure S3). The proliferation rates differ significantly between diagnosis and relapse (p < 10^{−}^{6} in Kruskal–Wallis test), while selfrenewal does not differ significantly (p ≈ 0.7 in Kruskal–Wallis test).
 Received January 23, 2014.
 Accepted February 20, 2014.
© 2014 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/3.0/, which permits unrestricted use, provided the original author and source are credited.