Mathematical models of leukaemia and its treatment: a review

Leukaemia accounts for around 3% of all cancer types diagnosed in adults, and is the most common type of cancer in children of paediatric age (typically ranging from 0 to 14 years). There is increasing interest in the use of mathematical models in oncology to draw inferences and make predictions, providing a complementary picture to experimental biomedical models. In this paper we recapitulate the state of the art of mathematical modelling of leukaemia growth dynamics, in time and response to treatment. We intend to describe the mathematical methodologies, the biological aspects taken into account in the modelling, and the conclusions of each study. This review is intended to provide researchers in the field with solid background material, in order to achieve further breakthroughs in the promising field of mathematical biology.


Introduction
Leukaemia is a cancerous disease in which blood cells display abnormal proliferation and invade other tissues, affecting the production and function of other blood cells. This malignancy is the most common cancer type in children from birth to 14 years of age and accounts for around 3-4% of all cancers diagnosed in developed countries [70]. Almost half a million  [11]. Generally, blood cancer overall 5-year survival in adolescents and young adults is about 50%. Although survival in children is higher and improving, blood cancer is still the major cause of cancer death in paediatric patients.
Most types of blood cancer start in the bone marrow, which is where blood is produced. In most blood cancers, the normal development process, starting from stem cells and leading to a hierarchy of more differentiated cells, is perturbed by the uncontrolled abnormal growth of specific types of blood cell.
There are three major types of blood cancer. Leukaemias are caused by the rapid production of lymphocytes and myelocytes (both white blood cells) in the bone marrow, which can also circulate in the blood stream. On the other side, lymphomas are a type of blood cancer affecting the lymph nodes and comprising abnormal lymphocytes. This type of white blood cells are involved in the adaptive immunity, producing specific antibodies during host defence. These cells multiply and collect in lymph nodes and other tissues and impair the lymphatic system's functionality to remove unnecessary fluids from the body and fight infections. Finally, myeloma is a cancer of the plasma cells also in the bone marrow, which produce disease-and infection-fighting antibodies.
It is well understood how blood cells differentiate from stem cells into more specialised cells, as represented schematically in Fig. 1. At the top of the hierarchy governing normal haematopoiesis there are the haematopoietic stem cells (HSCs) [91,121]. Multipotent haematopoietic stem cells can give rise to either lymphoid or myeloid progenitors. Lymphoid progenitors can generate either lymphoblasts, which will become B or T lymphocytes, or Natural Killer cells. These are all part of the immune system. Myeloid progenitors can also lead to a broad variety of cells, including erythrocytes, thrombocytes, or other cells of the non-specific immune system.
Although the classical understanding of haematopoiesis has considered cell types to be discrete compartments, current knowledge of the process [56] considers the evolution of cell types as a continuum process (see Fig. 2). This is because haematopoietic cells acquire lineage  [56] features through a continuous process involving the expression of different characteristic molecules [118].
In this framework, the type of cell that becomes cancerous determines the specific type of blood cancer. For instance, leukaemia can be either myeloid (or myelogenous), or lymphoid (or lymphoblastic, or lymphocytic). Also, leukaemias can be distinguished by the maturation stage of the transformed cells. Acute leukaemias affect blast cells (specifically immature, young lymphocytes and myelocytes), and grow very fast. Chronic leukaemias cause an accumulation of mature cells, leading to slowly growing cancers. Thus, there are four different classes of leukaemias: Acute Lymphoblastic leukaemia (ALL), Chronic Lymphocytic leukaemia (CLL), Acute Myelogenous leukaemia (AML) and Chronic Myelogenous leukaemia (CML) [31]. However, it is not completely clear whether the hierarchical organisation is preserved in blood cancers. Myeloid leukaemias seem to be hierarchically organised, whereas acute lymphoblastic leukaemias are not [9,42,90].
The origins of the mutations are essential to understanding self-renewal and differentiation fractions for cancer cells [84,90]. These probabilities could be explained by some of the basic hallmarks of cancer [38], such as sustaining proliferative signalling, resisting cell death, immortality, evading growth suppressors, and metastasis. The detection of these hallmarks is essential to tailoring treatments, which depend on classifying each patient within risk groups [16]. Currently, patients are assigned to a risk group depending on several factors, including the cell's morphology, the results of molecular or biochemical analysis, and the so-called flow cytometry techniques [36,64]. This is done by taking samples of the bone marrow (where the haematopoiesis process occurs) are obtained and characterised in terms of immunophenotypic patterns [117], which can be standardised [116].
Mathematical modelling may offer a new perspective in Oncology, specifically in blood cancers, with a huge potential to develop new strategies to characterise tumours and personalise treatments [3,13,86]. The cancer hallmarks relevant to each tumour type, coming from a specific cell lineage and a certain maturation stage, can be accounted for in mathematical models, usually as parameters to be estimated or as equations which model the dynamics of blood cancer development.
Blood cancers, and specifically leukaemias, have been one of the first types of cancers that has been thoroughly studied by applied mathematicians. It should be noted that there are many mathematically-grounded studies published in this field in high impact medical and general-science journals.
Leukaemias are a 'global' disease of the bone marrow, and as such spatial effects are usually ignored. They can be modelled mathematically, in an initial approach, using ordinary differential equations. More complex models have used partial differential equations, but to describe the evolution of some kind of trait or subpopulation, rather than spatial variables. Furthermore, blood cell counts are an easy way to gather information about the evolution of the disease. Putting the data together has led to substantial interest in the disease from modellers and clinicians managing the disease.
Only a small fraction of the data available during routine clinical procedures is used for diagnosis, and incorporated into the models developed so far. This review focuses on the role of mathematical models based on differential equations. However, mathematical techniques for (big-)data analysis [100] also have huge potential for providing answers to specific questions of relevance to leukaemia, whether alone or in combination with other mathematical methods. For instance, some studies have pointed out their potential use in avoiding expert manual gating of the data to identify leukaemic clones [1], analysing mass cytometry data [123], or predicting treatment response [26].
Our review is intended to expand the available literature on blood cancers [20,31] to incorporate more studies and greater detail, by focusing on leukaemia. Our plan in this paper is as follows. Firstly, we summarise mathematical models based on differential equations describing the growth of myeloid leukaemias. This focus reflects the fact that myeloid leukaemias are the commonest among adults. The only models that exist for lymphoblastic leukaemia concern treatment. We then review mathematical models for different types of leukaemia treatment. Finally, we discuss the results and summarise our conclusions.

Mathematical models of myeloid leukaemias
Myeloid leukaemias arise from alterations of cells of the myeloid lineage, and is considered a clonal disorder of the haematopoietic stem cells (HSCs). The condition may lead to an increase in myeloid cell, erythroid cell or platelet counts, not only in peripheral blood but also in the bone marrow. As described above, the two general types are chronic myeloid leukaemia (CML) and acute myeloid leukaemia (AML), depending on the maturation stage of the cells. In CML cells mature during the chronic phase, while in AML blast cells fail to mature, generating large amounts of blasts, i.e. immature cells [60,103].

Stem-cell based models of myeloid leukaemias
Stem-cell based models for myeloid leukaemia are based on mathematical models of the normal blood generation process, called haematopoiesis. The role of stem cells in cancer was recently reviewed in [111] in terms of mathematical models which can characterise cell behaviour in normal cell development. For blood cells, an important haematopoiesis model was proposed by Marciniak et al. [65]. The main assumption of this model was that the   (1). Cells were grouped into n different maturation stages, c i with i = 1, . . . , n. As cells mature, and the j index increases, their proliferation rates p c j increase, whereas the self-renewal fractions a c j decrease process of differentiation, i.e., the ability of a cell to change from one type to another, was described in several discrete maturation stages, beginning with stem cells as the first stage of maturation.
As cells mature, their proliferation rate increases, while the self-renewal fraction lowers, where self-renewal was understood as the probability of having the same properties and fates as their parent cell. This process is summarised in Fig. 3. The model includes different cell subpopulations with n different maturation stages and feedback signalling to regulate haematopoiesis.
The mathematical model describing the dynamics comprises a set of ODEs for the several compartments of normal cells (c i ) and another set of ODEs for the leukaemic cells (l j ) where c i = c i (t) denotes the density (or number) of healthy cells in each maturation stage i = 1, . . . , n, p c i are the proliferation rates of healthy haematopoietic cells in mitosis, a c i,max are the self-renewal fractions, and d c i the death rates for every cell maturation stage. The notation is analogous for the leukaemic cells, l j = l j (t) for j = 1, . . . , m, and the constants p l j , a l j,max and d l j for j = 1, . . . , m. Feedback signalling was described in that study using the cytokine effect function s(t). Cytokines are small proteins which assist in regulating fraction chemical signalling in cells. Cytokine concentration is modelled by the equation where k c and k l are the signalling regulation strength, for both normal and leukaemic cells, respectively. These parameters are sensitive to the number of mature healthy and leukaemic cells, c n (t) and l n (t). This signalling was assumed to control the dynamics of cell proliferation and differentiation in the mathematical model. Figure 4 shows an example of evolution towards the homoeostatic equilibrium of the healthy haematopoietic cell compartments for n = 6. The model of Eq. (1) in [106] was built on the basis of the haematopoiesis model of [65]. The main conclusion of the mathematical study of [106] was that both self-renewal fractions and proliferation rates could be indicators of poor prognosis. Similar models were also studied in [109], where some mathematical properties, including linear stability analysis, and necessary and sufficient conditions for the expansion of malignant cell clones, were studied for related models.
A similar model by the same group [105] described the differentiation process as a twostage process, but considered instead the multi-clonal nature of leukaemia, the feedback processes and the role of treatment. The study performed numerical simulations for 'in-silico' virtual patients, and obtained estimated parameters from the tumour growth data of two real patients. The researchers concluded that self-renewal might be a key mechanism in the clonal selection process. It was also stressed that late relapses could originate from clones that were already present at diagnosis, a question that has been the subject of discussion in the biomedical research literature. Stem cell self-renewal has been reviewed in terms of their impact on the dynamics of cell populations in [110], concluding that a high self-renewal fraction can lead to faster cancer growth.
A similar model, [108], accounted for genetic instability through the inclusion of the possibility of mutations, an essential hallmark in cancer evolutionary dynamics. Through comparison of patient data and simulations, the authors highlighted the fact that the selfrenewal potential of the first emerging leukaemic clone would have a major impact on the emergence of clonal heterogeneity so that it might serve as a biomarker of patient prognosis. A recent study of the group [59] on acute leukaemias formalised the clonal selection dynamics via integro-differential equations. They concluded that clonal selection was driven by the self-renewal fraction of Leukaemic Stem Cells (LSCs), constructing numerical solutions based on patient data parameters from the existing literature. These simulations showed that high self-renewal for LSC clones was a marker of stability in the presence of interclonal heterogeneity.
The model set out in Eq. (1) was further used in [120] to study feedback signals from myelodysplasic syndrome (MDS) clones and their effect on normal haematopoiesis. The model was fitted using serum samples from 57 MDS patients and five healthy controls. On the basis of the numerical simulations, the authors reached the conclusion that a high selfrenewal fraction of MDS-initiating cells may be critical for the development of the disease. It was conjectured that remission could be achieved if this parameter could be lowered.
Considering the dependence of leukaemic cell to cytokines, the model (1) is compared in [107] to a mathematical model including cytokine-independent leukaemic cell proliferation. In it, leukaemic cells are not controlled by cell signalling as in Eq. (1g), but instead a death rate is included that increases with the number of cells in the bone marrow, and acts on all cell types residing in bone marrow. This allows the authors to explain unexpected responses in some patients, such as blast crises or remission without chemotherapy. This was done by assigning patient data to two different groups that differ with respect to overall survival: those with cytokine-dependent or cytokine-independent leukaemic cell populations.
In [5] a mathematical model extension to the system (1) was developed to study chemotherapy dynamics for AML. The authors explained the cytarabine and anthracyclinelike chemotherapy effect by including specific terms in the respective equations for LSCs (l 1 ) and mature leukemic cells (l 2 ). The effect of cytarabine was included with a negative term −k cyt p l l 1 only for the mitotic phase, while the anthracycline effect was considered for the mitotic phase (−k anthra p l l 1 ) and mature cell phase (−k anthra l 2 ). This is, These equations were analogously constructed for HSCs (c 1 ) and mature healthy cells (c 2 ). The term d(x) = 10 −10 · max 0, x − 4 · 10 9 cells/kg represented an additional death rate describing the fraction of bone marrow cells dying because of overcrowding, which was dependent on the cell number x = c 1 + l 1 + l 2 . Finally, signalling s was chosen only to be dependent on the healthy mature cell number as s(t) = 1 1+k c c 2 . In particular, the authors compare the following: (A) a single-induction course with 7 days cytarabine and 3 day of anthracycline-like treatment; (B) a 7 + 3 course and a bone marrow evaluation that leads, in case of insufficient leukemic cell reduction, to the provision of a second chemotherapy course. Depending on the leukemia growth, the authors simulated both therapy schemes. The main conclusion was that therapy (A) was better for achieving a first complete remission while therapy (B) was better thought for lower therapy intensity and therefore to obtain less side effects.
Due to experimental and ethical limitations, the role of the stem cell niche in human leukemia is mostly unknown and has been particularly addressed for hematological malignancies. In [122], the authors collected bone marrow aspirates derived from 61 AML patients and 11 healthy donors, and so constructed a mathematical model including competition of healthy and leukemic cells for niche spaces and dislodgement of healthy cells from the niche by leukemic cells. Cells were subsequently divided in leukemic (l j ) and non-leukemic (c j ). Leukemic stem cells (LSCs) were denoted as l 1 , and non-leukemic hematopoietic stem cells (nl-HSCs) were denoted as c 1 . Furthermore, they considered progenitor (c 2,i , l 2,i ) and precursor cells c 3,i , l 3,i (limited either by a number of n 1 or n 2 divisions, respectively) and mature cells (c 4 , l 4 ). The subindex i denoted the number of divisions since entrance into their corresponding state. For stem cells, after mitosis, the probability that a randomly chosen niche space which the progeny tries to occupy is empty was denoted as p e , while being occupied by a LSC was p l . This last being dislodged had a probability function p r (η c /η l ), dependent on the fraction η c /η l of niche affinity η c of nl-HSCs over niche affinity η l of LSCs. The constructed model was based on these parameters to account for the flux to differentiation from the nl-HSCs as well as LSCs. Progenitor and precursors cells followed the structure of the model from [65], with cell death d only included for the mature cell compartments and signalling s. The results from the model suggested that competition of non-leukemic and leukemic stem cells for niche spaces is responsible for the decline of non-leukemic hematopoietic stem cells before relapse, indicating leukemic stem cell persistence. The authors, who suggested that competition of HSC and LSC for niche spaces also takes place in humans, also extended this work in [112], considering the same data and contributing to AML risk stratification. They developed this mathematical model of the human stem cell niche in AML, which was able to predict different risk groups, based on whether the proliferation rates from HSC and LSC are similar or not and whether the proliferation rate for LSC is high enough. These assumptions were compared to the number of diagnostic HSC (CD34+CD38−ALDH+) count, which was proposed as a prognostic tool in [122] if this number is in a threshold range. As these counts were similar for patients with different proliferation parameters, they concluded that a simple algorithm based on a cut-off value for HSC will not be sufficient for risk stratification.
More generally and not specifically for AML, a cooperative model of HSC homeostasis in which stem and niche cells mutually interact was proposed in [7]: for S p = S p (t)) the number of proliferative HSCs, N u = N u (t) the number of unoccupied niche cells and C = C(t) the number of HSC-niche cell pairs. In this model, λ denotes the division rate, and δ and δ the loss rates of proliferate and niche-bounded HSCs. Parameter μ > 0 denoted the net death rate of unoccupied niche cells, while ν is the rate of proliferation of HSC-bound niche cells. Finally, constants k 1 and k 2 described association and dissociation of the corresponding cells types. These parameters were obtained from available experimental data. The authors calculated a quasi-steady state approximation of the model and studied substantial perturbations of the initial conditions. They were then able to describe homeostatic recovery of the HSC compartment after irradiation versus apparent lack of recovery after HSC ablation. Theoretically, they also proposed that the outflux of differentiated cells from HSCs can be regulated by the affinity of HSCs for niche cells. Finally, in [4], migration of hematopoietic stem cells between the blood and bone marrow within a host were described via an ODE model. Cells were considered to be in bone marrow niches up to a fixed number N of niches. The number of cells in peripheral blood counts and in the bone marrow niches at a given time were denoted as s i and n i , respectively. The subindex i denote the number of host/wildtype cells for i = 1, while i = 2 denote mutant / donor cells. They considered the following ODEs model: where, for i = 1, 2, parameters d i were detachment rates, a i attachment rates, β i reproduction rates and δ i death rates. The model, parametrised using existing empirical findings from mice from literature, captures scenarios of clonal dominance and stem cell transplantation. The authors found that chimerism can be improved by injecting the host with multiple small doses, and accounted for the detectable levels of expansion for neutral and advantageous clones coming from one HSC.

Cell-cycle-based mathematical models of myeloid leukaemias
In some CML patients, symptoms may recur [71]. This is why periodicity is specifically studied for this disease. Thus, several authors considered the cell cycle in order to explain periodicity. The cell cycle is the process regulating cell division. It is a multi-stage process including, firstly, mitosis (M), the process of nuclear division; and a stage called interphase, the interlude between two M phases. In the interphase, three different substages occur: the G 1 phase, in which the cell prepares DNA synthesis; the S phase, where DNA replicates; and the G 2 phase, where the cell prepares for mitosis. Any cell, before going the S phase, can enter a resting state called G 0 , where the cell becomes quiescent and remains in a non-proliferating stage. This process is summarised in Fig. 5.
Many mathematical models have considered different aspects of the cell cycle [124]. However, many of those models, arising from the so-called systems biology approach, are quite complex. Due to the periodic nature of the cell cycle in proliferating cell populations, several mathematical models have tried to account for this cycling behaviour in a simplified form. Periodicity and other dynamic behaviours of haematological diseases are reviewed in [33]. The G 1 phase prepares for DNA replication and synthesis in the S phase. G 2 prepares cells for the mitosis phase M. While in G 1 , cells can become quiescent, entering a resting phase G 0 Specifically, the model in [61] described the dynamics of blood multipotential stem cells. Their approach was to write equations for a population N (t) of cells in the resting phase G 0 and another P(t) of proliferating cells, described by where N τ = N (t − τ ), τ being the cell cycle time. The function β(N ) = β 0 / (1 + (N /N * ) n ) is the mitotic re-entry rate, i.e. the rate of cell entry into proliferation, where β 0 , N * , n are parameters. The parameter δ is the total differentiation fraction from the G 0 phase, and γ is the fraction of irreversible cell loss from all portions of the proliferating-phase stem-cell population. Taking values for these parameters from the literature, the authors concluded that the origin of aplastic anaemia and periodic haematopoiesis could be related to irreversible cell loss from the blood multipotential stem compartment. The authors in [62] studied Eq. (5), to describe the existence and stability of long-period oscillations of stem cell populations in periodic chronic myelogenous leukaemia. This was made possible by studying a contractive return map, such that a fixed point of the return map gave a stable periodic solution of the model equation. This was computed in such a way that there was no analytic formula for the periodic solution in the limiting case n → ∞.
Other work based on the (5) model, such as [23], gives estimates of the model parameters for a typical normal human, and explored the changes in some of these parameters necessary to account for the quantitative data on leukocyte, platelet and reticulocyte cycling in 11 patients with Periodic Chronic Myelogenous leukaemia (PCML). Their results indicated that the critical model parameter changes required to simulate the PCML patient data were an increase in the amplification in the leukocyte line, an increase in the differentiation fraction from the stem cell compartment into the leukocyte line, and the rate of apoptosis in the stem cell compartment. In a companion study [22], they found that the parameter changes Fig. 6 Schematic view of the cell compartments in [97]. A, G 0 , M, R and B represent different cell compartments, α, β, γ are the corresponding coupling rates between them, and λ is the irreversible blood cell loss that mimic untreated cyclical neutropenia correspond to a decreased amplification (increased apoptosis) within the proliferating neutrophil precursor compartment, and a decrease in the maximal rate of re-entry into the proliferative phase of the stem cell compartment. The case of granulocyte colony stimulating factor treatment was also studied. Safarishahrbijari and Gaffari [101] used the equations for red blood cells and platelets from [23] and for leukocytes from [33] to identify parameters in PCML. The inclusion of new parameters resulted in a better fit of clinical data and from the data extracted from both platelet and leukocyte models.
Pujo-Menjouet and Mackey [89], performed a local stability analysis of the model (5) and found the conditions for Hopf bifurcation to occur. Periodic oscillations were studied depending on five haematopoietic stem cell parameters: the mitotic rate sensitivity, the maximal rate of cell entry into proliferation from the resting G 0 phase, the differentiation and apoptosis rate and the time to entry into mitosis. Extensions of this work [37], have proven that, under periodic treatment, there is a periodic solution with the same period. This could be related to the observed oscillatory behaviour of blood cells' counts under treatment in CML.
A different type of models to describe myeloblastic leukaemias have been constructed on the basis of the work of Rubinow and Lebowitz [97]. The model itself was based on granulocytopoiesis, also studied by these authors in [96]. In this first work, qualitative analysis was performed, supporting evidence for alterations which presumably occurs in cyclic neutropenia. For both models they considered four compartments for healthy cells as shown schematically in Fig. 6: the active A and G 0 cell compartments, representing the proliferative pools, and the maturation M and reserve R cell compartments, which finally ended in the blood pool B. For the leukaemic cells, only active and G 0 cells were considered, of which only a certain fraction were released into the blood, with no further maturation stages.
For this model, and in terms of myeloid leukaemia, the presence of a leukaemic population destabilises the homoeostatic state of the normal population, which is stable in the absence of leukaemic cells. In [98], the authors found differences between normal and leukaemic cell populations but including treatment into the model from Fig. 6: firstly, the recovery rate was higher for normal cells, as compared to leukaemic cells from the action of cytotoxic treatment. Secondly, the S-phase duration was different for the two populations. This led the authors to the conclusion that, for patients with a "slow" growing leukaemic cell population, remission could be achieved with one or two courses of treatment, whereas for those with a "fast" growing leukaemic cell population, a similar aggressive treatment achieved remission only at the cost of great toxic effects on the normal cell population.
These mathematical models described both the processes of normal blood and myelogenous leukaemia development. Fokas et al. [32] did the same, but for CML. The authors showed how CML cells could ultimately outnumber normal cells. They used the model to study the relationship between proliferation and maturation and proposed a solution to the apparent contradiction between decreased proliferation and increased production, by assuming that a greater fraction of CML cells is produced by division rather than by maturation.
Another mathematical model [35] included details on cyclins D, E and B, a family of proteins that help to control the cell cycle. Their production has a direct influence on the transition of a cell in the G 0 , G 1 and G 2 phases, respectively. Flow-cytometry data profiles for three leukaemia cell lines were analysed in this study (K-562, MEC-1, and MOLT-4, from AML, CLL and ALL patients, respectively). For the S phase, DNA replication was considered, as it is key before a cell can produce new daughter cells. The authors assumed that G = G(C E , t) was the number of cells in G 0 /G 1 at time t with a cyclin E content C E . Similarly, they denoted by S = S(DN A, t) and M = M(C B , t) the number of cells in the S and G 2 /M phases that had DNA content (represented as the variable DN A) and cyclin B content C B at time t, respectively. These assumptions taken together led to the model where r G→S , r M→G are the transition fractions from G 2 /M to G 0 /G 1 and from G 0 /G 1 to the S phase, respectively. Good agreement was found between experimental results and the model simulations. The authors claimed that the model could help in the identification of heterogeneous leukaemia clones at diagnosis and post-treatment, and that it could have the potential to predict future outcomes in response to induction and consolidation chemotherapy as well as relapse kinetics. This is due to the fact that the model allows the prediction of backward and forward culture evolution given known cell line-specific cell-cycle kinetics and initial conditions.

Other data-based mathematical models of myeloid leukaemia
Myeloid leukaemia models are the most studied in the literature. [102], for example, describes acute myeloid leukaemia (AML) using a multi-lineage multi-compartment model of the haematopoietic system and feedback via cytokines and chemokines. Analysis of the model suggested that self-renewal probabilities, mitotic rates and cytokine growth factors produced in peripheral blood determined leukocyte homoeostasis. The mitosis rate of cancer was found to be the parameter with the strongest prognostic value. A comparison of three mathematical models that describe CML progression and aetiology was undertaken in [63]. The authors sought to identify which models could provide the best description of disease dynamics and their underlying mechanisms. The first considered the following dynamic system where x 0 were HSCs, x 1 healthy progenitors, x 2 differentiated cells, y 0 LSCs and y 1 differentiated leukaemic cells, with parameters a, b, c, d, e as the corresponding self-renewal, production and death rates. Also, k was the carrying capacity and z = x 0 + x 1 + x 2 + y 1 + y 2 the total number of cells. The second model, in [68], was a shorter version of the (16) model, to be presented in detail later. The third was [34], which allowed competition between HSC and LSCs. This latter model was based on the following ODEs: The healthy cells x i and leukaemic cells y i were considered at different stages i = 0, . . . , 3 of differentiation and a compartment of quiescent cells was also added for each type, x q and y q . Parameters α x and α y denote the rates in which immature cells enter this quiescent compartment, while β x and β y are the rates in which they return to the maturation stage i = 0. Parameters r x and r y denote the stem cell division rates, while and the corresponding parameters a, b, c for also denote this for more mature stages. Cell death is denoted with parameter d. Finally, the functions φ x = 1/(1 + p x (x 0 + y 0 )) and φ x = 1/(1 + p y (x 0 + y 0 )) allow that both leukaemic and normal cells remain at a constant abundance once they reach the steady state, which can be obtained properly by setting the constants p x and p y . The authors found that it was not possible to choose between the models based on fits to the data of 69 patients who had experienced relapse or remission of the disease. They suggested experiments directly probing the haematopoietic stem-cell niche that could help in choosing the best model. Finally, [69] described another model of cancer initiation for CML. The authors assumed that the clonal expansion of mutant cells is given by a logistic equation where a is the time since mutation happened, x(a) the frequency of mutant clones with r relative fitness, and N the total cell population with generation time τ . q was the rate of detection and u the probability per cell division of producing a mutant cell. Letting c = r −1 τ , and b = N u c r , the probability of detecting cancer before time t was given by for with m a small probability that the first mutant arose early. Interestingly, this simple model, based only on the Philadelphia translocation, gave rise to cancer incidence curves with exponents of up to 3 as a function of age. This behaviour had been previously thought to be associated necessarily with three mutations, two of which were unknown. Thus, the model proved that CML incidence data were consistent with the hypothesis that the Philadelphia translocation alone could cause CML.

Other studies of myeloid leukaemias
Cancer initiation and maintenance are typically assumed to be related to cancer stem cells (CSCs) [15,119]. Two models for cancer initiation have been derived using this assumption. The first is a genetic mutation model, where mutations determine the phenotype of the tumour. In this conceptual framework different mutations may result in different tumour morphologies, even when starting from the very same stem cell. Cells inherit the molecular alteration and regain the ability for self-renewal, which leads to a population of cancer cells. The second model assumes that different cells serve as cells of origin for the different cancer subtypes, the so-called CSCs. This model proposes that oncogenic events occur in different cells, and these produce different kinds of cancer. In this model, self-renewal potential is limited for the CSCs. Both conceptual models are shown schematically in Fig. 7.
In [53] a stochastic model was constructed that considered drug resistance for CML, where the probability of treatment failure was approximated by for M 0 the initial non-mutant cells, n the quantity of drugs used, u the probability of mutation after cell division, and finally, the measurable parameters L, D and H as, respectively, the rate of growth, death and the drug-induced death rate. From the analysis of the mathematical model, the authors claimed that although drug resistance prevented successful treatment, resistance could be overcome with a combination of three targeted drugs. For the study of 38 patients diagnosed with Acute Promyelocytic Leukemia, the dynamics of leukocytes in peripheral blood and the effect of International Consortium on Acute Leukemia protocol treatment was studied in [126] by means of the following model  [121] where N = N (t) represents the number of normal leukocytes and A = A(t) is the number of leukemic cells. Function D = D(t) stands for chemotherapy concentration at time t. The reproduction of normal and leukemic cells was given respectively by r N and r A , while μ N and μ A represent the natural mortality rates. Parameter K A is the carrying capacity of leukemic cells. Constants β 1 and β 2 denoted the inhibitory effect of leukemic over normal cells and vice versa, while β 3 is the rate in which normal cells can mature into normal cells by the effect of all-trans retinoic acid (ATRA), a vitamin A derivative whose use have resulted in long-term overall survival rates. This treatment is described by a unit step function, for at specific final treatment time t F . Regarding chemotherapy, the considered daunorubicin doses of ρ = 60 mg with decay rate τ = 26.7 h. The model parameter values were selected from the literature and, for two selected patients, the model produces a good fit to the clinical data. The authors performed several simulations for each of these parameters, with the resulting model behavior ranging from very aggressive to non-aggressive leukemias. Induction treatment was modelled following the corresponding protocols and observing some similarities among the populations considered. The model was able to capture the dynamics of peripheral blood and was proven to be useful for the prediction of relapse. Finally, several authors have built models of leukaemias using graph-theoretical methods. Graphs can be used to describe the hierarchical organisation observed in haematopoiesis, as seen in Fig. 1. In [17], a PDE model of haematopoiesis was parametrised on a graph using publicly available RNA-Seq data in a high-dimensional space. The high-dimensional data were later reduced to R 2 or R 3 using reduction techniques, such as principal component analysis, diffusion maps and t-distributed stochastic neighbour embedding, and a PDE model on a graph G was constructed. u(x, t) denoted the cell distribution at the differentiation continuum space location x ∈ G and time t. Then, for every cell distribution u k (x, t) on an edge e k , the cell density was modelled with advection-diffusion-reaction equations for x ∈ e k = a k b k , where each edge e k was parametrised from a k to b k , and the following functions were considered: R k = R k (x) as the cell proliferation, V k = V k (x) as the advection coefficient, and apoptosis and diffusion terms D k = D k (x) and w k = w k (x), which respectively describe cell fluctuation and width of a narrow domain around an edge. Using this model, the authors performed simulations consistent with the evolution of AML populations. A similar approach was used in [25], where the graphs constructed presented the essential properties of functioning bone marrow.

Imatinib and its basic mathematical modelling
CML has been intensively studied in terms of therapy based on Imatinib. This drug is a 2phenyl amino pyrimidine derivative that inhibits a number of tyrosine kinase (TK) enzymes. Imatinib is specific for the TK domain in ABL (the Abelson proto-oncogene), c-kit and PDGF-R (platelet-derived growth factor receptor). In chronic myelogenous leukaemia, the Philadelphia chromosome leads to a fusion protein of ABL with the breakpoint cluster region, termed BCR-ABL. Imatinib decreases the BCR-ABL activity. CML treatments have been strongly influenced by the appearance of imatinib [28], that is now the standard first-line treatment against the disease. It is a very effective drug with up to about 70% of people having a complete cytogenetic response (CCyR) within 1 year of starting imatinib. After a year, even more patients will have had a CCyR. Many of these patients also have a complete molecular response (CMR). The capacity of the drug to impair the proliferation of leukaemic stem cells was the basic assumption behind the mathematical model of Michor and co-workers [68]. The model also included the development of resistance to therapy and was based on the following system of differential equations: Here, x i denotes the different populations of normal cells, y i the imatinib-sensitive leukaemic populations and z i the tumour clones resistant to imatinib. The indexes i = 0, 1, 2, 3, denote the subpopulations of stem cells, progenitors, differentiated and terminally differentiated cells in each compartment. The rate constants for each cell type (x, y, z) are given by a, b and c, and d i are the death rates for i = 0, 1, 2, 3. Cell division rates are given by r y and r z . The parameter u is the fraction of resistant cells produced per cell division. Finally, λ = λ(x 0 ) is a decreasing function of x 0 describing homoeostasis of normal stem cells. It models the feedback signals controlling haematopoiesis. Data from 169 CML patients were used to fit the mathematical model in [68]. The authors obtained numerical estimates for the turnover rates of leukaemic progenitors and differentiated cells and showed that imatinib dramatically reduced the rate at which these cells are being produced from leukaemic stem cells. They showed that the probability of harbouring resistance mutations increases with disease progression as a consequence of an increased leukaemic stem cell abundance, and proposed that the time to treatment failure caused by acquired resistance is given by the growth rate of the leukaemic stem cells. Their bottom line was that multiple drug therapy is especially important for patients who are diagnosed with advanced and rapidly growing disease.
A simplified version of the model (16) was studied in [29] by considering only the stem cell (0) and differentiated cell (1) compartments of healthy (x) and leukaemic (y) cells, i.e.
where φ = 1/ [1 + c x (x 0 + y 0 )] and ψ = 1/ 1 + c y (x 0 + y 0 ) are homeostasis functions for normal and tumour stem cells respectively, and c x , c y are Michaelis-Menten parameters. By a combination of analysis and simulation, the authors discussed how any successful therapy would require the eradication of the pool of leukaemic stem cells; otherwise, progressive disease is very likely. Thus, successful therapeutic agents must enhance the death rate of this rare population of cells. Therapies designed to target mitosis of malignant stem cells could not eradicate the disease quickly. Nevertheless, there has been some controversy surrounding the potential effectiveness of imatinib to achieve remission [67].
In [47], the immune response targeting leukaemic cells was added to Eq. (16). Using experimental data from the literature, a mathematical model was fitted in which immune response was described by delay differential equations. The authors considered that T cells targetting leukaemic cells could prevent relapse, and combine with imatinib therapy. The more simplified model in Eq. (17) was later used by [83] to study and numerically simulate treatment interruptions as a potential therapeutic strategy for CML patients. In many cases, strategic treatment interruptions led to the elimination of leukaemic cells in silico. The biological rationale behind these interruptions is to reduce the drug side effects, drug resistance, or to stimulate the immune response. The conclusion was that strategic treatment interruptions could be a feasible clinical approach to enhancing the effects of imatinib treatment for CML.
A number of extensions of the (16) model have been developed for CML. For example, in [79], four levels of cell differentiation were included and studied for the BCR-ABL1 gene, necessary for the pathogenesis of CML. In that study, data from 290 patients were used, 92 of them treated with dasatinib, 75 with nilotinib and 123 with imatinib. All treatments elicited similar responses. Another extension of the model was described in [41], with a focus on more theoretical aspects, including a stability analysis, and an existence proof for positive solutions.
The global dynamics of normal and CML haematopoietic stem cells and differentiated cells were also studied in [2]. The dynamic was assumed to be governed by the following system of Lotka-Volterra equations where x 0 (t) represents haematopoietic normal stem cells (HSC), x 1 (t) normal differentiated cells and y 0 (t) and y 1 (t) describe the same subpopulations of cancer cells. In Eq. (18) n, m, r , q are division rates, d 0 , d, g 0 , g death rates, K the carrying capacity and α ∈]0, 1[ is a constant. The production rates for differentiated cells are given by d 2 and g 2 . Several optimal control problems were solved for imatinib, whose effect on the division and mortality rates of cancer cells produces a suboptimal response. The effect of cyclic combination of two drugs in CML was studied in [52], and the modelling led to the conclusion that treatments should start with the stronger drug, and the weaker one should have cycles of longer duration. An interaction model between naïve T cells (mature T cells from thymus), effector T cells (cells which actively respond to stimuli) and CML cancer cells was described in [72], where Latin hypercube sampling was used to estimate parameter values due to the lack of data. This is a statistical technique for generating parameters from a multidimensional distribution. In their conclusion, the authors explained that the growth rate of CML and the natural death rate were the most important parameters, suggesting that treatment for CML patients should focus on these rates. Any drug with a high cost that is included in the model could be studied in order to obtain optimal treatment, and reduce not only radiation but also financial benefits. This model was later used in [8], focusing on cancer x = x(t) and effector y = y(t) cell population dynamics, by considering a combined treatment with imatinib and the interferonalpha (IFN-α) therapy. This last is a protein whose activation produces a cytogenetic response in CML patients. The model considered the following ODEs where β 1 , β 2 were the respective reproduction rates, K the maximal tumour population, η 1 , η 2 Michaelis-Menten terms and γ 1 , γ 2 the cell loss rates due to interaction. The death rate for effector cells is μ y , while tumour death is modelled by the constants ωγ 3 and a function h(t). This function is modelled as h(t) = t − θ e −λt , so that the influence of drugs tends to zero over time. The dose of IFN-α is modelled as in α γ 4 , which increases the effector cell population with a delay τ of about 7 days. The stability analysis proposed, as well as the results obtained, were able to describe the influence of two types of the treatment. The authors claimed that the dose of IFN-α has an inhibitory effect on the number of cancer cells, but its replacement with another type of treatment should be considered in order to avoid resistance. Finally, [76] studied optimal control problems for CML, in a model with a molecular targeted therapy such as imatinib. Naïve T cells, which are already differentiated T cells, but are precursors for more mature cells called effector cells, were also included in the model. The cancer cell population was then activated by the presence of the CML antigen. Aiming to minimise the cancer cell population and the detrimental effects of the drug, they found that a high dose level from the beginning was optimal. Also, combination therapy was better than single dosing.

Modelling the effect of quiescence on imatinib treatments
Quiescence, which corresponds to the G 0 phase of the cell cycle, and its relationship to drug therapy (in this case, imatinib) is an important factor in leukaemia because quiescent cells might not be affected by therapy, as drugs target proliferative cells, and a possible relapse may occur.
Imatinib treatment was studied using Roeder model [93][94][95] accounting for quiescent and proliferative cell compartments. Firstly, in [95] a stochastic model of haematopoiesis was developed. On the basis of that model, another was built to describe imatinib-treated patients [94].
A more advanced model based on partial differential equations (PDEs) was studied in [93]. This model considered quiescent and cycling stem cells, denoted by A and , respectively. The authors included a cell-intrinsic function a(t), which determined the affinity of a cell for residing in A or . With a(t) ∈ [a min , a max ], a quiescent stem cell would enter the cell cycle with probability ω and a cycling cell would become quiescent with probability α. These terms were modelled as where the sigmoidal functions f ω and f α were defined by for specific values of the parameters ν j , μ j , for j = 1, 2 and the scaling factors N ω and N α . The dynamics of the HSCs, quiescent (A) and proliferating ( ), were governed by these equations: The functions n A = n A (a, t) and n = n (a, t) represent the cell densities at affinity a and time t within A and , respectively. Also, v A = v A (a) and v = v (a) were the corresponding velocities that make v A · n A and v · n the corresponding cell fluxes for each compartment. Finally, τ was a parameter which simulates average cell division depending on cell cycle duration. Equations (20e) and (20f) were the basis for studying leukaemia and how the imatinib treatment affect its dynamics, in a highly efficient way when it comes to huge cell populations. They considered the dynamics for every cell subpopulation in the following system: ∂n (1) ∂t +v (1) ∂n (1) ∂n (2) ∂t +v (2) ∂n (2) where the super indexes i represent the different cell populations as normal cells (i = 1), imatinib-affected leukaemic cells (i = 2) and non-affected leukaemic cells (i = 3). Induced cell death is denoted by a constant r deg , while the constant r inh denotes the proliferation inhibition on the proliferating cells n (2) . The model in Eq. (20) was proved to qualitatively and quantitatively reproduced the results of the agent-based approach for imatinib-treated patients in [94]. This was fitted to 894 peripheral blood samples, where the authors claimed that the therapeutic benefits of imatinib can, under certain circumstances, be accelerated by being combined with proliferation-stimulating treatment strategies. [48] described an extension of the (20) model. This was done by considering the cycling cells to be dependant, among other variables, on a counter c(t), that indicates the position in the cell cycle, with a 49-h cell cycle. An imatinib treatment was then incorporated into the model. The authors conclude that PDE formulation provided a more efficient way of simulating the dynamics of the disease. In fact, in simulations of imatinib treatment, the PDE and the discrete-time models diverged more, as in this case a continuous-time description of the disease dynamics may be more realistic than discrete-time models. This latter model was later extended [19] by including feedback from cells and asymmetric division for stem cells and precursors. The general idea for this work was also to combine imatinib with a drug that induced cancer stem cells to cycle. Furthermore, the fact that many patients do relapse after being taken off imatinib motivates the study methods by which this therapy can be improved. Doumic-Jauffret et al. [30] performed a stability analysis of the model in [48], where the authors could set differences between AML and CML in terms of transition from stable equilibrium to unstable periodic behaviour.

Some examples on ibrutinib treatment for chronic lymphocytic leukaemia
Ibrutinib is a kinase inhibitor mainly used in CLL, used to block the protein signalling that allow cancer cells to multiply. Regarding CLL, in [125], 10 patients serial complete blood counts and estimated numbers of tumor cells in tissue were used to fit a mathematical model. The authors aimed to quantitatively understand the redistribution during ibrutinib treatment of lymphocytes in CLL. They denoted the number of CLL lymphocytes in the tissues and blood by x and y, respectively. The model was given by which described the time evolution of these populations during treatment. In the tissue compartment x, CLL cells were assumed to die with rate d 1 , and redistribute to the blood with rate m. In the blood, CLL cells died with, rate d 2 . Finally, parameter c described phenomenologically the observation that the rate of decline of lymphocytes slows down over time and stabilizes around a steady state that can be higher than healthy levels. To fit the model, several cell counts were used to fit while blood counts y, while CLL in tissues were fitted by computed tomography scans at two time points during treatment. The authors concluded that ibrutinib causes a significant amount of cell death in tissue (a larger death rate than in blood), and that a relatively small fraction of the tissue cell burden redistributes to the blood, which is 3 and 5 times higher than without treatment. Nevertheless, they did not include the bone marrow, which could increase tissue cell burden and lower the estimated fraction. System from Eq. (21) was later used in [12] to check the effects of ibrutinib on leukemia cell proliferation and death in a clinical study of 30 patients with CLL, now considering the use of stable isotopic labeling with deuterated water. This study confirmed the inhibition role of ibrutinib in CLL cell proliferation and its promotion of high rates of CLL cell death.

Whole body mathematical description of leukaemia and its treatment
Leukaemia treatment may affect blood flux in several tissues on the body. In order to understand the behaviour of these body parts during therapy, we set out a highly descriptive model of leukaemia, chemotherapy and blood flux throughout the entire body [85]. The inflow rate of drug j is where u j is the drug dose over duration j . This equation was then incorporated into the following equation, which models drug concentration in the blood C B, j : In this equation, V B is total patient blood volume, and Q i the blood flow in the organs i, such as heart (H ), liver (Li), bone marrow (M), lean muscle (Le) and kidneys K , and so C i, j was Fig. 8 Schematic view of the use of mathematical models to help in patient treatment. Personalized characteristics from the patient's disease, as well as specific ratios for each of the biological processes involved, could be implemented in the mathematical model as parameters and equations. After the simulation, several parameters might arise that could be advantageous for specific treatment protocol, optimised for the specific patient, who could benefit from the personalised drug. This cycle could be useful for following the disease in the patient.
for every organ i and drug j, where k k, j is the urine excretion rate, k L, j the elimination rate in the liver, and V i,T the volume of organ tissue where drug metabolism occurs. This model, along with many others, are useful in clinical terms, as it could provide guidance for optimising treatment for each patient in terms of their characteristics, as explained in Fig. 8. This pharmacokinetic model was reinforced by a pharmacodynamic model, which took into account the effect of the drug. Drug concentration at the location of the tumour, which for leukaemia would be the concentration of drug in the bone marrow (C M, j ), was considered for the j effect of the drug as the function effect j . It was included in the cell cycle as where P y was the cell population in phase y (G 1 , S, G 2 , M) and k y the transition term from phase y to y + 1. Although these equations are described in a general sense, for the specific case of chemotherapy cycles of intravenous (I V ) daunorubicin (DN R) and cytarabine (Ara − C), typical drugs in leukaemia treatment, the reactions occurred at a subcutaneous level. That is, the drug is injected under the skin and not below muscle tissue. This drug and its subcutaneous effect have also been addressed in other studies, such as [44], fitting data from 44 AML patients during consolidation therapy to a pharmacokinetic mathematical model, obtaining optimised treatment schedules. However, the authors of [85] considered, when simulating the subcutaneous effect of the therapy, that Eq. (22c) could then be replaced by the following two: where S is the subcutaneous tissue drug delivery, k a the absorption delay and k b the drug bioavailability. However, the simulations performed were adapted for two acute myeloid leukaemia patients. Sensitivity analysis method was applied on the model to identify the most crucial parameters that control treatment outcome. This model included data from two patient case studies and two chemotherapy protocols. The model is an example of the application of mathematics based on previous models, including interpretable parameters, performing sensitivity analysis and optimizing therapy protocols.

Mathematical models of acute lymphoblastic leukaemia treatments with cytotoxic drugs
The current standard treatment of acute lymphoblastic leukaemia involves different treatment stages: induction, consolidation, re-induction whenever needed, and maintenance [24]. The aggressiveness of treatments depends on the classification of patients into risk groups: standard, average or high (Fig. 9). The goal of the induction stage is to achieve a rapid reduction in tumour cell numbers. Next, the consolidation phase should ideally remove any trace of leukaemic cells in flow-cytometry or blood cell count studies. Re-induction is considered whenever leukaemic clones reappear early. The maintenance phase is administered when the first two phases are completed, and is intended to kill any possible remaining non-measurable quantities of cancer cells. Every phase includes specific treatments, the doses and timings of drugs depending on the patient's risk group. Using one mathematical model or another to describe therapy may lead to a different understanding of how treatment affects cells in terms of relapse [55]. For example, if relapse occurs and we consider a Cancer Stem Cell (CSC) model, a drug might not affect CSCs, or might only affect cells with specific mutations (in the genetic mutation model). This can be better seen in Fig. 10. In ALL, two drugs are used as part of these treatment phases: 6-Mercaptopurine and Methotrexate. Some mathematical models of their actions are now summarised.

Describing the effect of mercaptopurine
Mercaptopurine (6MP) is an antimetabolite antineoplastic agent with immunosuppressant properties. It interferes with nucleic acid synthesis by inhibiting purine metabolism and is used, usually in combination with other drugs, in the treatment of or in remission maintenance programmes for leukaemia.
A mathematical model of the effect of 6MP in leukaemia cells was described in [80]. In this model, the number of cells in the G 0 /G 1 -phases was denoted by G; S in the S-phase, and M in the G 2 /M-phase. The suffixed variables G I , S I and M I were the equivalent variables for the thioguanine (TGN) nucleotides, which were considered as the main active metabolites. That is, the most active molecules involved in the metabolic process. Apoptotic cells A, and non-viable cells N (cells that are unable to live), were also included in the model. The equations for the viable phases of cells are Finally, the apoptotic and non-viable phases are modelled as The dynamics of the model are summarised in Fig. 11. The model parameters describe the transition between phases, except for f ∈ [0, 1], which measures the fraction of cells continuing the cell cycle after TGNs were incorporated into the cell DNA. To estimate these parameters, the model was fitted to data for different cell lines treated with MP. The mathematical model provided a quantitative assessment to compare the cell cycle effects of MP in cell lines with varying degrees of MP resistance.
In a different study [43], semi-mechanistic mathematical models were also designed and validated for MP metabolism, by studying red blood cell mean corpuscular volume (MCV) dynamics, a biomarker of treatment effectiveness and leukopenia, a major side effect related to very low percentages of leukocytes. The model was validated with real patient data obtained from literature and a local institution. Models were individualised for each patient using nonlinear model-predictive control. The authors claimed that their approach could be implemented with routinely measured complete blood counts (CBC) and a few additional metabolite measurements. This would allow model-based individualised treat-ment, as opposed to a standard dose for all, and to prescribe an optimal dose for a desired outcome with minimum side-effects.

Mathematics of methotrexate treatments
Methotrexate (MTX) is an antimetabolite of the antifolate type. It is thought to affect cancer by inhibiting dihydrofolate reductase, an enzyme that participates in the tetrahydrofolate synthesis. This leads to an inhibitory effect on the synthesis of DNA, RNA, thymidylates, and proteins. A first mathematical model of MTX effect in ALL was constructed in [82]. The authors based their approach on the fact that within cells, MTX is metabolised to more active methotrexate polyglutamates (MTXPG), and these polyglutamates are subsequently cleaved in lysosomes by glutamyl hydrolase (GGH). GGH acts as either an endopeptidase or an exopeptidase. To better define the in-vivo functions of GGH in human leukaemia cells, GGH activity was characterised with different MTXPG substrates in human T-and B-lineage leukaemia cell lines and primary cultures. Parameters estimated from fitting a series of hypothetical mathematical models to the data revealed that the experimental data were best fitted by a model where GGH simultaneously cleaved multiple glutamyl residues, with the highest activity on cleaving the outermost or two outermost residues from a polyglutamate chain. The model also revealed that GGH has a higher affinity for longer chain polyglutamates.
Further research led to the development of an improved model in [81]: Finally, Le et al. [57] constructed a model involving a combination of several drugs, for chemotherapy-induced leukopenia in paediatric ALL patients. The model accounted for the action of both 6-MP and MTX and their cytotoxic metabolites 6-TGNc and MTXPGs during maintenance therapy. The equations were built on the basis of the previously discussed models [43,81]. The model predicted WBC counts for the available patient data surprisingly well, given the large variation of individual response patterns in the clinical data. The mathematical model and algorithmic procedure proposed could be used to guide personalised clinical decision support in childhood ALL maintenance therapy. Another model based on Refs. [43,81] gave rise to a compartmental model in [45], including pharmacokinetics and a myelosuppression model for ALL, considering both 6-MP and MTX. The model was crossvalidated with data from 116 patients, and simulations of different treatment protocols were performed to exploit the optimal effect of maintenance therapy on survival.

Immune response mathematical models
Immunotherapy is a type of therapy that stimulates cells within the immune system in order to help the body fight against cancer or infections. Interactions between cells are key in understanding processes such as, for example, proliferation or resource competition between cells. The immune system is one way in which the body may influence external agents and a greater understanding of it could be useful in fighting leukaemia.
An extension of the model already described, (16), was introduced in [21], where the CML populations were distributed as stem cells (y 0 ), progenitor (y 1 ) and mature leukaemic cells (y 2 ). In this study, the concentration of immune cells was also included and denoted as z. The authors designed a mathematical model integrating CML and an autologous immune response to the patients' data by considering the following system where a 0 , b 1 represents transition terms; d z and d i , for each cell type i = 1, 2, 3, denotes cell death; and a logistic growth for progenitor cells y 1 was included, with a reproduction rate r . The immune system action rate μ was included in the mass action term "μ y i z" in the last term of the leukaemic population equations from Eqs. (25a)-(25d). The proliferation of the immune cell pool included a constant factor s z and was activated by mature leukaemia cells with the term "α y 3 z" in Eq. (25e). These latter terms included an inhibition of the immune cells expansion, as they were divided by "1 + y 2 individual autologous immune response. The use of immunotherapy was then considered to be a useful complement to the usual treatment, playing a significant role in eliminating the residual leukaemic burden. A general mathematical model for tumour immune resistance and drug therapy was proposed in [27]. By including tumour cells, immune cells, host cells and drug interaction, an optimal control problem was constructed. This would provide a basis for the study of leukaemia immune cell interaction, shedding some light on the modelling for B leukaemia. For B-cells, fundamental in both acute and chronic lymphocytic leukaemia diagnosis, a more extensive model was presented in [75], including four different cell populations in the peripheral blood of humans: B cells, able to bind to antigens which will initiate antibody responses; NK cells, critical to the immune system; cytotoxic T cells, able to kill cancer cells; and helper T cells, which may help other immune cells by releasing T cell cytokines. This model was considered a tool that may shed light on factors affecting the course of disease progression in patients, and focused on sensitivity analysis for parameters and bifurcation analysis. Based on [27], an immunotherapy approach was considered in [92] by developing a model focused on B and T lymphocytes and their relation with a chemotherapeutic agent. The ODE system for this model was the following where N = N (t) represented the neoplastic B lymphocytes, I = I (t) the healthy T lymphocytes (this is, the immune cells), and Q = Q(t) the amount of a chemotherapeutic agent in the bloodstream. In Eq. (26a) N follows a logistic growth with a proliferation rate r , and dies due to both interaction with immune cells at a rate c 1 and with the chemotherapeutic agent at a rate μ. Immune cells in Eq. (26b) have a constant source s 0 and die naturally at a constant rate d and also due to interaction with cancer cells at a rate c 2 , and with drugs at a rate δ. However, there is a production rate ρ of immune cells stimulated by cancer cells. Both N and I have Michaelis-Menten terms with rates a, γ and b. For the case of the chemotherapeutic agent Q in Eq. (26c), λ is considered as the washout rate of a given cycle-nonspecific chemotherapeutic drug with λ = ln(2)/t 1 2 , where t 1 2 is the drug elimination half-life. Finally, the functions s(t) and q(t) are source terms, which can be considered to be constants. These parameters were all taken from the literature and claimed to simulate CLL behaviour. This model reinforces the option of combining treatments such as chemo-and immunotherapy, where the first may decrease cells to a point where immune cells may act.
A model for AML was considered in [77] by including the role of leukaemic blast cells (L), mature regulatory T cells (T reg ) and mature effector T cells (T eff ), this last also including cytotoxic T lymphocytes and Natural Killers. As hypothetical immunotherapy, instantaneous increases and decreases in T eff and T reg simulated infusion of effectors and depletion of regulatory T cells. The interaction between these cells was modelled by the following system: where a L , a T eff , a T reg represented influx rates, and d L , d T eff , d T reg the decay rates. Intracellular interactions were modelled as Hill functions with threshold constants (k 1 , k 2 , k 3 ) with strength p. Two existing steady states were found for this model in [77], corresponding either to leukaemia diagnosis or relapse, or to complete remission. The authors considered that the model explained the influence of the duration of complete remission on the survival of patients with AML after allogeneic stem cell transplantation. In [78], Monte Carlo simulations of trajectories in the phase plane were performed for the prior model. The authors concluded that the duration of complete remission influences the survival of patients with AML after allogeneic stem cell transplantation. Besides, they generated relapse-free survival curves, which were then compared with clinical data, yet to be further assessed by future immunotherapy clinical studies.

Including interleukins in mathematical models
Interleukins (ILs) are a group of cytokines first seen to be expressed by white blood cells (leukocytes). The immune system depends on interleukins as these signals between cells are useful for acting against several pathogens. The interaction between the actively responding effector cells E = E(t), tumour cells (T = T (t)) and the concentration of the cytokine IL-2 (I L = I L (t)) was the basis for the latter study, influenced by [50]. The reason behind the modelling of this cytokine is due to the fact that IL-2 might boost the immune system to fight tumours. This was described via the following system: In this model, c was antigenicity or ability to provoke an immune response, 1 μ 2 was the average natural lifespan, a the loss of tumour cells by interaction, μ 3 the degraded rate of IL-2, and s 1 , s 2 were treatment terms. The fraction terms were of the Michaelis-Menten form, to indicate saturation effects. The function r 2 (T ) could be described as a constant for linear growth, or with limiting-growth as logistic or Gompertz terms. With this model, the authors concluded that with only IL-2 treatment, the immune system might not be enough to clear tumours. These and other models were reviewed in [113] in terms of equilibrium points, considering T lymphocytes and their interaction with other cells, and it was found that there are two stable equilibrium points, one where there is no tumour, and the other where there is a large one.
Interaction between cells via interleukins was also studied in [14], as IL-21 is being developed as an immunotherapeutic cancer drug. Its effect has been studied in relation to Natural Killer (NK) cells, and CD8 + T-cells, which have the ability to make cytokines, with the model the function of the memory factor m and g(n) the dynamics of tumour cell number, which is constructed separately for each tumour type according to the observed growth curves. Parameters were estimated in terms of certain values from the literature, so that simulations were run to show IL-21 as a promising antitumour therapeutic. For more immunotherapeutic approaches towards cancer modelling, we highlight the work in [88], where some general aspects of cancer were also reviewed, including diffusion, angiogenesis and invasion. Finally, for the case of immune response to leukaemia, other studies have been undertaken, though not specially in the form of an ODE or PDE system. Some numerical simulations were run in [51] by proposing an integro-differential equation model. This study proposed a new possibility for defining the activation states for cancer, cytotoxic T and T helper cells. Using these definitions, the authors suggested that it would be easier to organise experiments suitable for measuring cell states. They also claimed that cell-mediated immunity is one of the most crucial components of antitumour immunity. Immune T-cells were studied in [18] in terms of a stochastic model from which was derived a Fokker-Planck equation. Stability analysis and behaviour of the solutions of the model led to the conclusion that more accurate simulations of cancer genesis and treatment were needed. Lastly, in [99], cytotoxic T cells were dynamically and structurally analysed in terms of a Boolean network model for T cell large granular lymphocyte leukaemia. Nineteen potential therapeutic targets were found, and these were versatile enough to be applicable to a wide variety of signals and regulatory networks related to diseases.

Novel therapies for leukaemia models: CAR-T cells
Immunotherapy based on chimeric antigen receptor T (CAR-T) cells has been especially successful in patients who did not respond to the usual types of chemotherapy. This technique is based on the patient's own T-cells, which are extracted from them, genetically modified and reinfused. This modification allows T-cells to kill tumour cells in a more effective way than the usual chemotherapies.
We have designed a general model for CAR-T cells in [87] considering several cell compartments. Firstly, for T-cell leukaemia, the number of CAR-T cells was denoted by C, leukaemic T cells by L, and normal T-cells by T . The dynamics of the model were as follows The parameter ρ L represents leukaemic proliferation rate, while ρ C represents stimulation of CAR-T cell mitosis after encounters with target cells; τ C is the finite lifespan of CAR-T cells; the parameter α represents death due to encounters with CAR-T cells; parameter ρ I is the external cytokine signal strength used for division of CAR-T cells; finally, the function g(T , L, C) denotes the rate of production of normal T cells, assumed to contribute only at a minimal residual level. The stability analysis of the cell dynamics leads to several conclusions: firstly, CAR-T cells allow for control of T-cell leukaemia in the presence of fratricide; secondly, the initial number of CAR-T cells injected, as well as re-injections, does not affect the outcome of therapy, while higher mitotic stimulation rates do; lastly, tumour proliferation rates have an impact on relapse time. A second, similar model was constructed for B cells, in [58], where CAR-T and now leukaemic B cells where again denoted as C and L, but the inclusion of mature healthy B cells B, CD19 -B cells P, and CD19 + cells I was considered. The initial autonomous system of differential equations was where parameters ρ C , τ C , ρ L and α were the same as the considered in the previous model. Parameter ρ B = βρ C , where 0 < β < 1, accounts for the fact that represented B cells are located mostly in the bone marrow and encounters with CAR-T cells will be less frequent. Parameters ρ P and ρ I represent growth rates for P and I cells, while τ I , τ B and τ P represent the finite lifespan of I , B and P cells respectively. A signalling function s(t) = 1/[1+k s (P + I )], with k s > 0 was constructed as in [65], also including the asymmetric division rates a P  Fig. 12 Illustration of the dynamics in model (31). B cells (in blue) develop in the bone marrow, arising from progenitor CD19cells (P), then turning, with rate τ P , into CD19 + cells I and reaching, with rate τ I a mature stage of healthy B cells B, finally dying after a time τ B . During this process, a signalling effect s(t) affects the proliferation rates of the early stages ρ P and ρ I . Leukaemic cells L develop in the bone marrow with rate ρ L , invading this tissue as well as the blood compartment. CAR-T cells C attack mature B cells and leukaemic cells with rate α, also inducing growth, with rate ρ C . In the bone marrow, they also attack CD19 + cells I , with a lower rate αβ. This interaction induces growth with rate ρ β . Solid lines represent cell growth and change between compartments. The dotted line represents the natural death of the healthy B cells. Dashed lines represent cell death due to CAR-T cell interaction and a I for P and I . This general model is reduced, in order to understand the dynamics of the expansion of CAR-T cells and their effect on the healthy B and leukaemic cells, neglecting the contribution of the haematopoietic compartments. Parameters are estimated from the literature and the main conclusion obtained is that not only does CAR-T cell persistence depend on T-cell mean lifetime, but also that reinjection may allow the severity of relapse to be controlled. The dynamics of the model from Eq. (31) are summarised in Fig. 12. A further iteration of this line of models is presented in [66] which splits CAR-T's population in effector and naive memory cells, focusing on bone marrow dynamics. Considering parameters from the literature, we concluded the importance of the product attributes rather than the initial amount of CAR-T cells, tumor burden or second dosages. A general model taken from the literature and applied to CAR-T cells is set out in [46]. The authors include a population of leukaemic cells (abnormal cells c = c(t)), and w = w(t) as the population of white blood cells or immune cells. Besides, they consider a compartment of healthy/susceptible blood cells s = s(t), which after contact with cancer cells become dysfunctional or "infected" i = i(t). The dynamics are modelled as where a 0 , β 0 , k 0 , and b 0 are the natural death rate of susceptible blood cells, infected cells, cancer cells, and immune cells, respectively; for susceptible cells, A is the recruitment rate and β is the loss rate of susceptible blood cells due to infection; β 1 is the decay rate parameter of infected cells; k is the constant recruitment rate of cancer cells, while k 1 and b 1 are the loss rates of cancer and immune cells due to interaction; finally, parameter B is considered as the external re-infusion rate of immune cells (CAR-T). This model was studied in terms of stability, and it was observed that the external re-infusion of immune cells by adoptive T-cell therapy reduces the concentration of cancer cells and infected cells in the blood.
With the success of T-cell-engaging immunotherapeutic agents, there has been growing interest in the so-called cytokine release syndrome (CRS), as it represents one of the most frequent serious adverse effects of these therapies. CRS is a systemic inflammatory response that can be caused by a variety of factors, such as infections and certain drugs. A more specific model that included the action of cytokines was studied by considering Tisagenlecleucel, a personalised cellular therapy of CAR-T cells for B-cell ALLs, associated with a high remission rate. It was modelled in [73] by considering the interaction of a CAR-T cell population c T = c T (t) with B-cell leukaemic population l = l(t), as well as with healthy B cells h = h(t), both marked with CD19, a characteristic of B lymphocytes. Other circulating lymphocytes were denoted as c = c(t), while the number of cytokines, key to understanding inflammatory processes, was generally considered as s = s(t). The dynamics of the model were as follows: where Eq. (33a) represented the dynamics of CAR-T cells with growth rate d 1 and natural death rate d 2 , while α 1 and β 1 were cell death given by interaction with leukaemic and healthy cells, respectively. Equation (33b) includes a growth rate of leukaemic cells k and a cell death α 2 by interaction with c T . Equation (33c) described a logistic growth of healthy cells with rates a and b, as well as a natural death rate d 3 and death β 2 due to interaction with c T . Circulating lymphocyte dynamics were considered in Eq. (33d) to have a constant input λ, death rate σ and growth dependant on c T , attenuated via a Hill function with constants α 3 and β 3 . Finally, for Eq. (33e), cytokines were secreted at a maximum rate α 4 and altered by a negative feedback mechanism corresponding to the term −β 4 s. Furthermore, the stimulation of CAR-T cells increased the levels of cytokines with rate d 4 and a constant m from the correspondent Hill function. Optimal control theory was applied for this model, controlling the injection of CAR-T cells and cytokines, to finally minimise the level of cancer cells and to keep healthy cells above a desired level.
Effector T cells are a group of cells including several T-cell types that actively respond to a stimulus. Following an infection, memory T cells are antigen-specific T cells that remain in the long term. This distinction is considered to help understand the dynamics of CAR-T cells in several models. For instance, a general description of Tisagenlecleucel was performed in [104], where data from 91 paediatric and young adult B-ALL patients were used for the analysis. The model describes the expansion of CAR-T cells up to a time T max , and then two phases: a first contraction phase, with rapid decline; and a second persistence phase, declining more gradually. This was represented by a dynamic system considering effector E and memory CAR-T cells M, as and M = 0, for T ≤ T max . After T max , effector cells rapidly decline at a rate α and convert to memory cells at a rate k, which decline at a rate β. However, before T max , only effector cells grow at a rate ρ and proportionally to a function F(t) which simulates the inclusion with step-wise functions of the co-medication of corticosteroids and tocilizumab (anti-IL-6 receptor antibody). This simple model was able to show the long-term persistence used in CAR-T therapies. The authors in [6] also considered a division between tumour T , effector CAR-T cells C T and memory CAR-T cells C M in the following model where f (T ) is the density dependence growth of tumour cells, and respectively for effector and memory CAR-T cells, we have the following: p C T (C T ) and p C M (C M ) as cell production functions, d C M (T , C M ) and d C T (T , C T ) as cell inhibition functions, and a C M (C M ) and a C T (C T ) as natural death functions. For this model, most functions were considered to be linear, except for the tumour growth function, considered to be logistic growth. Simulations were run for mice data found in the literature, showing different outcomes depending on tumour burden or initial therapy dose. The authors considered that a high CAR-T cell inhibition from tumour leads to tumour escape and absence of CAR-T cell memory. The same CAR-T cell division was considered in the model from [39], not only showing a distinction between effector and memory, but also between the cytotoxic (CD8 + ) and helper (CD4 + ) cells. Again, parameter values were not obtained from actual data, but from simulated clinical data. Their results suggest the hypothesis that initial tumour burden is a stronger predictor of toxicity than the initial dose of CAR-T cells. Also, the authors considered an inflammatory immune response regulated via a Hill function to maintain a realistic bound on the activation rate of T cells. This function gave rise to tumour-burden-correlated toxicity, while the correlation of CAR-T cell dose alone and toxicity was poor.

The authors in [49] developed a mathematical model with tumour B cells (B = B(t)), the normal memory T cells (N = N (t)) and CAR-T effector (E = E(t)) and memory CAR-T cells ((M = M(t))
where r N , r M are production rates, while K N and K M are carrying capacities with K M < K N . Effector and tumour cells have a predator-prey interaction term but with a predator rate coming from asymmetric differentiation from memory to effector CAR T cells. This rate r E is considered to be antigen-dependent but for simplicity assumed as linear. The rates γ E γ B are the effector-tumour interaction rates, where an effector or tumour cell dies upon interaction, respectively. The CAR effector death rate is denoted as d E and r B is the tumour reproduction rate. This model becomes stochastic when cell number are below a certain threshold and the authors conclude that therapy could be improved by optimizing the tumour killing rate and the CAR T cells' carrying capacity. The pharmacological model in [40] considered both the influence of CAR-T cells in inflammatory responses with cytokines (such as interleukins I L 6 , I L 10 or interferon I F N γ ), as well as the distinction between CAR-T cells into effector and memory cells. This was also done in order to understand toxicity related to cytokine release syndrome. In the model, the variable B represents CLL tumour B cells in peripheral blood (PB). CAR-T cells in PB are divided into effector E P B and memory M P B cells. This division is also performed for the CAR-T cells in the tissue compartments (E T and M T ). The complete mathematical model is shown in Fig. 13 Fig. 13 Illustration of the dynamics in model (37). The grey irregular shape represents tumour B-cells in CLL. Yellow triangles represent the inflammatory cytokines I L 6 = I L 6 (t), I L 10 = I L 10 (t) and interferon I F N γ = I F N γ (t). The green circle and pentagon represent respectively effector CAR-T cells from the PB and from the tissue. Squared, blue shapes represent memory CAR-T cells, also from PB and tissue. Solid lines represent promotion of cell production, while dotted lines represent cell loss due to natural death or due to encounters between cells. Short-dashed lines represent exchange between cell compartments, and finally long-dashed lines represent constant production in the cell compartments In this model, parameters r B , r E and r M represent growth rates, while d B , d E and d M are death rate constants, respectively for B, E P B and M P B cells. Parameter K BC is the is the effector CAR-T-mediated B-cell CLL degradation rate constant in peripheral blood. For the inflammatory immune responses we have, respectively for I L 6 , I L 10 and I F N γ the following constants: ρ endo IL 6 , ρ endo IL 10 and ρ endo IFN γ as endogenous synthesis rates; parameters ρ max IL 6 , ρ max IL 10 and ρ max IFN γ as production rates; and finally, d IL 6 , d IL 10 and d IFN γ are the natural death rates by the activated CAR-T cells. Constants a G and b G are the inhibitory parameters of I L 10 on I F N γ production. PB and tissue compartments are distributed via rate constants k in and k out after intravenous infusion. Peripheral blood effector memory CAR-T cells are activated via activation rates a E and a M . Finally, function f (B P B ) is chosen as a Hill function such that f (B P B ) = B P B B P B +h , with h the half-saturation constant of the tumour. This model was adjusted to data from 3 patients obtained from the literature. Its main conclusion is that toxic inflammatory response is correlated to disease burden, i.e. the number of tumour cells in bone marrow, and not with CAR-T cells doses, contrary to what is observed with most cancer chemotherapies. Other models have also considered these hypotheses, such as the discretised model in [115] for CAR-T cells. In this study, a logistic equation of growth was considered to explain the interaction between CAR-T cells and malignant tumour cells. The binding affinity of the CAR-T cell construct (the so-called single-chain variable fragment) and the antigenic epitope (the molecule binding to the antibody) on the malignant target was considered a critical parameter for all T-cell subtypes modelled. Both studies show the need for CAR-T cell doses to account for tumour burden, which would require a relatively low number of infused CAR-T cells to achieve the desired target.

Theoretical studies of leukaemia treatment models
In previous sections we have described leukaemia growth and response to therapy models that are careful to account for experimental facts or available data. There have been also many studies of models that focus their attention more on methodological mathematical aspects, and provide insight of a more fundamental type. For instance, some of them do not specify which type of leukaemia or treatment they describe.
As an example, some optimal control problems for general leukaemia treatment models have been discussed in the literature. In [10] the authors describe the dynamics of a healthy cell population N (t), a leukaemic population L(t) and a drug h(t) governed by the equations for L(0) = L 0 , N (0) = N 0 , h(0) = 0 and γ h the drug dissipation rate. The effect of the drug was described differently for diseased and healthy cells by the therapy functions f l (h) and f n (h), respectively. Here, L a and N a were the maximum number of diseased and healthy cells respectively, and γ l and γ n were respectively the death rates for the two kinds of cells. Interaction between these subsets was expressed by the parameter c. Finally, the control function u(t) is the quantity of drug given to the patient. The authors solved the optimal control problem using the Pontryagin maximum principle. Later research provided additional results along these lines in [114], by using a non-Gompertz interaction term and several phase constraints. Analysis of the switching points was performed, as well as several simulations. Some optimal therapy protocols are shown by introducing a 'shifting-variable', which avoids the violation of the normal cell constraint.
Other studies have considered the combined effect of Haematopoietic Inducing Agents (HIA) and Chemotherapeutic Agents (CTA) on stem cells, with the goal of minimising leukopenia [74]. Proliferating (P) and non-proliferating cells (N ) were included in the model: for t < τ, where τ was the time for a cell to complete one cycle of proliferation, γ the apoptosis rate, and δ the random cell loss. The expression N τ stood for N (t − τ ), introducing a time delay into the equation, and β(N ) = β 0 θ n θ n + N n ; (39c) β H I A (P) = β 0,H I A θ m 1 θ m 1 + P m g H I A (t); (39d) β C (P) = β 0,C P w θ w 2 + P w g C (t); (39e) were Hill functions measuring the rate of cell re-entry into proliferation, the effect of HIA, and the effect of CTA on stem cells, respectively. Also, simulated the time decay of HIA. Finally, CTA time decay was modelled by Using this set of equations the authors found that HIA administration increases the nadir observed in the proliferative cell line compared with when CTA treatment alone is administered. This is significant in preventing patients undergoing chemotherapy treatment from experiencing secondary effects. Furthermore, the steady state value of the proliferating cells was found to be significantly lower in silico after CTA treatment. The model and accompanying analysis give rise to an interesting question: Is concurrent administration of an HIA during chemotherapy a prudent approach for reducing toxicity during chemotherapy? There is substantial clinical evidence to suggest that HIAs could be useful in cases of anemia. They argued that prophylactic benefits of HIAs use together with chemotherapeutic agents at the onset of treatment, although rational, should be balanced with the treatment cost and the risk that HIAs will cause adverse side effects such as venous thromboembolism and tumour progression.
In the context of therapy, we should also highlight the existence of hybrid models such as the one in [54], where cells are simulated as discrete elements whose dynamics depend on continuous intracellular and extracellular processes. Individual cells have the ability to move, grow, divide, differentiate, and die by apoptosis. Some intracellular processes can be explained by two proteins, Erk (E i ) and Fas, (F i ) in each cell i at a time t, which are involved in cell fate decision and are described by the following ordinary differential equations: The activation of both proteins depend respectively on the functions α(E po) and γ (Fas L), which are dependant on the extracellular concentration of Fas-ligand, denoted Fas L, and erythropoietin, denoted E po. Both Erk and Fas proteins are deactivated with rate constants a and d, with inhibition terms b and c. The term β E k i denotes a a positive feedback loop to activate itself with a cooperativity coefficient k and a rate β. Moreover, the concentration of Ara-C w i in each cell i was also described by the following equation dw i dt = k 1 (w − w i ) − k 2 w i (t) − r p − r dp + r da (41) where k 1 is related to the cell membrane penetration, k 2 accounts for degradation and w is the time-dependent extracellular concentration of Ara-C, assumed to be uniformly distributed in the bone marrow. Functions r p , r dp and r da depended on w i as a Michaelis-Menten term, and represent respectively the phosphorylation rate of Ara-C into Ara-CTP, Ara-CTP dephosphorylation rate and a deamination rate. Then, the authors ran simulations of the hybrid model for the action of Ara-C on either normal (with circadian rhythm) or leukemic erythroid progenitors. Using parameters and data from the literature, they were able to compare the leukemic cell number from AML patients with Ara-C intravenous infusion with the simulations from leukemic hematopoiesis from the hybrid model. Finally, they highlight the importance of using chronotherapeutic treatments to treat leukemia by analysing the delivery time on the outcome as well as the influence of the period of treatment.

Conclusion
Mathematical models have proved to be an essential asset in biomedicine. Haematological diseases are well suited to mathematical modelling, not only with differential equations, but also with stochastic models or other techniques. Therefore, there is a huge amount of data to combine with the mathematical models already in the current literature. Notably, singlecell information, molecular and omics data are revealing a new picture of the structure of biological systems, at many different levels. These data allow for a more detailed description of the heterogeneity present in these systems such as different cellular subtypes and functions.
Besides, single-cell and molecular data analysis allows for a more precise quantification of the variables involved. A few works have already integrated this data in the modelling process (see e.g. [17] or [79]), exposing at the same time the challenge that this recently unveiled heterogeneity poses. The use of multi-omics data implies several challenges from the clinical point of view, such as the need for standardized protocols for data collection, panel creation and data storage. This would require an ongoing communication with the modelling community in order to precise the format that optimizes both clinical care and research. Other issues related to clinical data are the inaccessibility of bone marrow samples in comparison with peripheral blood ones. This impairs our capacity to follow bone marrow specific processes and at the same time demands a deeper understanding of the relationship between these two compartments. This would require the incorporation of data from healthy patients, which also poses a challenge from the clinical perspective. More directly related to the modelling process, coupling absolute count data (e.g., concentrations) with heterogeneity measures (e.g., flow cytometry proportions, molecular biology data) would facilitate the quantitative analysis.
From the mathematical point of view, the incorporation of such data calls for more elaborated and probably more complex modelling approaches. The classical ODE framework may give way to elaborate PDE models in which the additional variable represents molecular or phenotypical information. Discrete and agent-based models may also be a powerful alternative since they allow for the specification of the different species or agents that participate in the process. The other advantage of these models is that, rather than equations, they make use of rules of behavior that are easier to translate for the biomedical community, facilitating the communication between both fields and thereby fostering the advancement of mathematical oncology.
With respect to future directions, the most recent literature is mainly concerned with immunotherapy and with the modelling of continuum states of differentiation, going beyond rigid, hierarchical schemes. Works related to the first have the opportunity to describe and simulate data coming directly from clinical trials, enhancing the parallel development of both research approaches (clinical and mathematical) and thereby promoting the employment of mathematical models during the establishment of new therapeutic options. With respect to the possibility of modelling continuum states, the opportunities lie in the above-mentioned necessity of integrating complex data. This invites the use of artificial intelligence methods and machine learning algorithms and coupling them with more classical approaches. In general, a majority of works focus on myeloid malignancies, probably due to their higher incidence, so another direction of improvement would be to translate those findings to lymphoblastic leukemias. This may require multicenter studies with standardized protocols in order to include a higher a number of patients. This highlights again the importance of close collaboration not only for improving personalized therapies based on higher quality datasets but also for possible mathematical modelling approaches. Ultimately, mathematical models could be refined by their inclusion in hospital protocols as a diagnostic or prognostic tool, and this can only be achieved by cooperation between the mathematical and medical world.