Clonal selection and therapy resistance in acute leukaemias: mathematical modelling explains different proliferation patterns at diagnosis and relapse

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 self-renewal. 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 multi-clonal composition of leukaemias. Model results imply that enhanced self-renewal may be a key mechanism in the clonal selection process. Simulations suggest that fast proliferating and highly self-renewing cells dominate at primary diagnosis, while relapse following therapy-induced remission is triggered mostly by highly self-renewing 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.


Introduction
Leukaemia is a clonal disease of the haematopoietic system leading to extensive expansion of malignant cells that are non-functional 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 stem-like properties [2,3]. A clone consists of genetically identical stem and non-stem 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][5][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 therapy-specific 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][11][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][20][21][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][24][25]. The novel ingredients of the models in this work are: (i) heterogeneity and multi-clonal 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 signal-independent activation of important cell functions [26][27][28] and by an increased cell degradation observed in leukaemic patients [29][30][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 short-duration 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 long-term 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 multi-clonal 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.

Material and methods
2.1. Mathematical models 2.1.

Model assumptions
The models used in this study are based on the models of healthy haematopoiesis proposed and analysed in [23,[34][35][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, so-called 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 post-therapy 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 two-compartment version of the multi-compartment system established in [23]. The model focuses on the maintenance of primitive cells and differentiation from undifferentiated, proliferating cells to differentiated, post-mitotic cells. In the case of healthy haematopoiesis, the proliferating cells are haematopoietic stem cells (HSCs), haematopoietic progenitor cells (HPCs) and precursor cells, the post-mitotic cells are mature cells, e.g. white cells. The two-compartment architecture is based on a simplified description of the multi-stages differentiation process. Nevertheless, as shown in [24,35,36] models consisting of two compartments capture the desired dynamics of the multi-compartmental 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 self-renewal, 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 self-renewal 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 post-mitotic 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 density-dependent 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 post-mitotic leukaemic cells.
Overcrowding is defined when marrow cell counts exceed the steady-state marrow cell count by two to three times. In this case, the death rate of post-mitotic 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][41][42], mediated by cytokines, such as G-CSF 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 self-renewal is a more efficient mechanism than the regulation of the proliferation rates. Similar conclusions have been drawn using the models of multi-stage 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 self-renewal depending on the level of mature cells.

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 self-renewal. The death rate of mature blood cells is denoted by d c 2 : 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 self-renewal). The fraction 2(1 À a c 1 (t)) p c (t)c 1 (t) 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. Self-renewal of a certain cell type at time t is assumed to be given as a maximal possible self-renewal of this cell type multiplied by s(t). Following [23,25], we chose s(t) ¼ 1/(1 þ k c c 2 (t)), which can be derived from cytokine kinetics [23]. The constant k c depends on the rate of extra-haematopoietic 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 a c max corresponds to the maximal possible self-renewal of HSCs. : (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.

Model 1
We assume that leukaemic cells depend on the same feedback signal as their healthy counterparts and that the post-mitotic 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 self-renewal. 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 post-mitotic cell compartments for each leukaemic clone. Let p l i denote the proliferation rate of mitotic cells in leukaemic clone i and a l i max the corresponding maximal fraction of self-renewal. By d l i 2 . 0 we denote the clearance rate of post-mitotic cells of clone i. Denote by l i 1 (t) the level of mitotic cells of clone i and by l i 2 (t) the level of post-mitotic cells at time t. These assumptions result in the following system of ordinary differential equations: : (2:11) The expression for s(t) is a special case ofs , where we assume that k l i ¼ k l 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, G-CSF [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 k l i 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 k l ¼ k c for the simulation of patient examples.

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][27][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][30][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 post-mitotic 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][57][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 lim x!1 d(x) ¼ 1: 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 x c 1 , where c 1 is the steady-state count of mitotic healthy cells.
Assuming that the haematopoietic cell lineage is regulated as described above, we obtain the following system of differential equations:

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.

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 steady-state 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 post-therapy 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 ODE-solvers from the MATLAB-software package (v. 7.8, The MathWorks, Inc., Natick, MA, USA) which are based on Runge-Kutta schemes.

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.

Properties of clones at diagnosis
Simulations indicate that the clones significantly contributing to the leukaemic cell mass have high proliferation rates and high self-renewal 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.

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.

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 self-renewal 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 self-renewal 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].

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 chemo-sensitivity 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.

Short-term expansion efficiency does not correlate with long-term self-maintenance
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 self-renewal potential dominate the leukaemic cell mass during an initial phase. If, with increasing leukaemic cell mass, self-renewal becomes downregulated, e.g. owing to occupation of bone marrow niche, eventually the cell clone with the highest affinity to self-renewal survives, although its proliferation might be slow. An example of time evolution during an early phase is depicted in figure 5.

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   Figure 5. First phase of leukaemic clone evolution: at the beginning fast proliferating clones with low self-renewal can dominate. They are later outcompeted by clones with high self-renewal, which is an advantage under high competition for niche spaces, needed for self-renewal. If there exist clones with high self-renewal and high proliferation, they will dominate during this first phase of leukaemic evolution. Each line type corresponds to one leukaemic clone. Blasts are immature cells used for diagnosis of leukaemias. In the course of the disease, blasts accumulate and outcompete haematopoiesis. Blast counts greater than 5% are considered as pathological [67]. The simulations are based on Model 1. rsif.royalsocietypublishing.org J. R. Soc. Interface 11: 20140079 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.

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 AML-typical 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, follow-ups and relapse. In addition, emergence and subsequent elimination of leukaemia driving mutations (FLT3, MLL-PTD) in the bone marrow cells were precisely monitored using molecular biology methods [68][69][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, FLT3-ITD mutation of a length of 39 bp is detected at diagnosis. This mutation becomes extinct and the relapse is triggered by two different FLT3-ITD 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 FLT3-ITD mutation and MLL-PTD mutation were detected at diagnosis. The MLL-PTD 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 anti-leukaemic drugs, if feasible, or for bone marrow transplantation.

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 self-renewal 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 self-renewal properties. Depending on efficiency of the therapy, it is also possible that clones with high proliferation and high self-renewal 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 self-renewal 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 self-renewal potential, which is an advantageous factor for expansion and survival. Otherwise, they would be outcompeted by slowly proliferating cells with high self-renewal. 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 rsif.royalsocietypublishing.org J. R. Soc. Interface 11: 20140079 the first time evidence that self-renewal potential is a major force behind this mechanism and that cells responsible for a relapse show high self-renewal 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 self-renewal. (iii) Cytotoxic treatment selects for cells with slow proliferation and high self-renewal. 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 non-leukaemic cells coexist at ratios compatible with sufficient haematopoiesis for long periods. 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 self-renewal 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][54][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 self-renewal.
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 self-renewal 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 rsif.royalsocietypublishing.org J. R. Soc. Interface 11: 20140079 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 self-renewal. One possible implication might be the application of cell-cycle independent drugs, such as those used in targeted therapies.
i.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 a c % 0:87 (A 10) this 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 d l i 2 ¼ 0:1. 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 d(x) ¼ d const: Á max {0, x À x max }. 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 x max % 2 c 1 , where c 1 is the steady-state count of mitotic healthy cells. In the simulations, d const. was set to 10 210 . 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 self-renewal 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 self-renewal 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 rsif.royalsocietypublishing.org J. R. Soc. Interface 11: 20140079 of each cycle. For simplicity, we did not model kinetics of single chemotherapeutic agents. Instead, the therapy-induced 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 . Dt, where Dt 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. non-preleukaemic 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 non-LSC leukaemic cells.
Let l i 1 (t) be the level of LSC of clone i at time t. The flux to mitosis is then l i 1 (t) p l i (t): Out of mitosis, we obtain 2a l i (t) p l i (t)l i 1 (t), where a l i (t) is the fraction of LSC self-renewal of clone i at time t. We assume that the fraction n of these cells is mutated, n takes into account replication errors in relevant genes and is assumed to be constant. The influx a i (t) of mutated LSCs due to new mutations occurring in clone i at time t is therefore 2a l i (t) p l i (t)l i 1 (t)n.
We obtain the following set of equations describing dynamics of clone i: d dt l i 1 (t) ¼ 2a l i (t) p l i (t)l i 1 (t)(1 À n) À p l i (t)l i 1 , d dt l i 2 (t) ¼ 2(1 À a l i (t)) p l i (t)l i 1 (t) À d l i 2 l i 2 (t) and a i (t) ¼ 2a l i (t) p l i (t)l i 1 (t)n: A similar system of equations has been obtained by Traulsen et al. [83]. As l i 2 is considered to be post-mitotic, we do not distinguish between cells that acquired a mutation during the divisions and those that did not.
The influx a(t) of mutated cells at time t is given by a(t) ¼ g þ P N(t) i¼1 a i (t), where g is the constant influx from the preleukaemic compartment and N(t) the number of leukaemic clones present at time t.
We accept the rate a(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 t 1 is a jump time of the inhomogeneous Poisson with rate l(t), then the next jump time t 2 can be generated by solving the equation Ð t2 t1 l(t)dt ¼ Àlog(1 À u) for t 2 . Here u is a uniformly distributed random variable u [ [0, 1] [86]. We further know that if u is uniformly distributed in u [ [0, 1], then 2log(1 2 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 Ð t1 t0 l(t)dt ¼ r 1 : 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. self-renewal and proliferation chosen randomly from uniform distributions (proliferation rate between 0.01 and 0.9, self-renewal 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 Ð t2 t1 l(t)dt ¼ r 2 , 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 self-renewal and high proliferation while at relapse cells show high self-renewal and reduced proliferation (electronic supplementary material, figure S3). The proliferation rates differ significantly between diagnosis and relapse ( p , 10 26 in Kruskal -Wallis test), while self-renewal does not differ significantly ( p % 0.7 in Kruskal -Wallis test).