A Mathematical Model of the Transition from Normal Hematopoiesis to the Chronic and Accelerated-Acute Stages in Myeloid Leukemia

Abstract

A mathematical model given by a two-dimensional differential system is introduced in order to understand the transition process from the normal hematopoiesis to the chronic and accelerated-acute stages in chronic myeloid leukemia. A previous model of Dingli and Michor is refined by introducing a new parameter in order to differentiate the bone marrow microenvironment sensitivities of normal and mutant stem cells. In the light of the new parameter, the system now has three distinct equilibria corresponding to the normal hematopoietic state, to the chronic state, and to the accelerated-acute phase of the disease. A characterization of the three hematopoietic states is obtained based on the stability analysis. Numerical simulations are included to illustrate the theoretical results.

1. Introduction

Mathematics can offer qualitative and quantitative tools for a better understanding, prediction, and control of biological processes. Particularly, such tools and models have been given for blood cell production process and hematological disorders. Among the first papers in this direction one can mention the works of Rubinow and Lebowitz [1,2], Mackey and Glass [3], Mackey [4], and Djulbegovic and Svetina [5]. For more recent contributions we refer to Fokas et al. [6], Neiman [7], Andersen and Mackey [8], Colijn and Mackey [9], Adimy et al. [10], Dingli and Michor [11], Kim et al. [12], Cucuianu and Precup [13], Doumic-Jauffret et al. [14], Komarova [15], Stiehl and Marciniak-Czochra [16], MacLean et al. [17,18], Radulescu et al. [19], Bianca et al. [20,21], Ragusa and Russo [22], and to the references therein. For some models concerning stem cell transplantation, we mention the papers of Vincent et al. [23], De Conde et al. [24], Kim et al. [25], Marciniak-Czochra and Stiehl [26], Precup et al. [27,28], Precup [29], and Stiehl et al. [30]. Some reviews that work on mathematical models for cancer, particularly for chronic myeloid leukemia, are the papers of Afenya [31], Michor [32], Foley and Mackey [33], and Clapp and Levy [34].

Hematological stem cells (HSC) from bone marrow are situated at the origin of the process of cell formation. They have the self-renewal capacity and the ability to differentiate and produce various types of blood cells. Perturbations of this complex biological process are at the root of several hematological diseases, such as chronic myeloid leukemia (CML), a malignant disease arising from mutant stem cells affecting the line of myeloid cells and progressing in three phases: chronic (also called indolent) phase, accelerated or transitory phase, and acute or blast phase. The separation of the last two phases being difficult, we shall refer to them as a whole, naming it the accelerated-acute phase.

An analogy can be drawn between the evolution of normal and abnormal stem cell populations and the spread of two species into ecological habitats. Mathematically this evolution is expressed by a differential system involving a number of specific parameters. Such a model for CML, proposed by Dingli and Michor [11], yields to a mathematical characterization of two hematological states: normal state and leukemic state, without being able to make distinction between the chronic and accelerated-acute phases of the disease. Further analysis of the model and its extensions for bone marrow transplantation have been undertaken by Precup and co-workers (Cucuianu and Precup [13], Precup et al. [27,28], and Precup [29]).

In this paper we propose a refined version of the Dingli–Michor differential system that allows us to make the distinction between the following three hematopoietic states related to chronic myeloid leukemia: normal hematopoietic state, chronic leukemic state, and accelerated-acute state. We assume normal hematopoietic state as a biological state in which the population of one or more mutant hematopoietic stem cells (mHSC) tends to zero due to random events that lead to the extinction of mutant cells, a fact demonstrated experimentally in stem cell lineages (see Jilkine and Gutenkunst [35], Driessens et al. [36], Klein et al. [37], Lopez-Garcia et al. [38], and Snippert et al. [39]). We accomplish our purpose by operating with distinct sensitivity parameters for normal and abnormal cells, instead of a single common parameter, like the one that was used in the original model. Our modeling choice is biologically justified by the asymmetry of mutual influences and interactions between normal and leukemic cells (Gou et al. [40]). The new model, also expressed by a two-dimensional differential system, has three nontrivial steady states, and their stability analysis shows that only one of them is asymptotically stable depending on the value of a cumulative parameter D that incorporates the growth rate, the cell death rate, and the sensitivity rate of the abnormal cells, and represents the equilibrium amount of abnormal cells. Values of D under some threshold correspond to the normal hematological state; values of D lying in some interval characterize the chronic phase, while larger values of D stay for the accelerated-acute phase (see Figure 1 below). In this way, from a mathematical point of view, the transition from one hematological state to another occurs as the result of a change of the basic leukemic cell parameters accumulated by $D$. The transition is progressive as the disease gets worse, and regressive as a result of treatment. For the last case, our analysis could be a guide to the improvement of therapeutic agents and strategies, as suggested in Precup et al. [41].

The paper is structured as follows: In Section 2, we consider the new model and we carry out its mathematical analysis. Based on this model, the normal hematopoietic state, and the chronic and accelerated-acute phases of CML are mathematically characterized in terms of parameter $D$. In Section 3, we provide some numerical simulations of the model. Next, in Section 4, the model is upgraded by six additional equations so that, as in Michor et al. [42] and Dingli and Michor [11], it becomes able to describe the cell evolution on four levels: primitive stem cells, progenitors, differentiated cells and terminally differentiated cells. The numerical simulations show the parallelism between the dynamics of primitive stem cells and that of the succeeding lines, which allows the analysis to be performed at any level, particularly for terminally differentiated cells for which laboratory data can be obtained easier. Finally, in Section 5, some medical remarks and conclusions are included. We conclude this introductory section with medical background addressed to those readers that would like to have more biological and medical information about hematopoiesis, malignant disorders, and related literature.

Medical Background

Hematopoiesis is the process of blood cell formation. The process starts in the intrauterine life in the mesoderm of the yolk sac and continues in the liver and the spleen during the second and seventh month. It then takes place at the level of the bone marrow, where it carries on after birth. During childhood, hematopoiesis takes place in almost all bones, gradually being replaced with growth by fat tissue. In adults, hematopoiesis occurs only in the pelvis, vertebrae, sternum (Howard et al. [43] and Young [44]), ribs, skull, proximal humerus, and femur epiphysis (Kaushansky et al. [45]).

Hematopoiesis and the differentiation process into various blood cell types can be viewed as an evolutionary tree that grows from one single hematopoietic stem cell (HSC). A HSC can have one of the following functions: it can renew itself; it can generate two other HSCs; or it can lose the ability of self-renewal and begin the maturation and differentiation pathway resulting two progenitor cells; or finally, it can give birth to a HSC and a progenitor cell. Progenitor cells are capable of initiating the differentiation towards one of the pathways that lead to the formations of various types of blood cells: common lymphoid progenitor (CLP) that will differentiate and maturate into B or T-lymphocytes, and the common myeloid progenitor (CMP) that will differentiate and maturate into leukocytes (white blood cells), erythrocytes (red blood cells), and platelets (Young [44] and Kaushansky et al. [45]). According to recent studies, there is growing evidence that hematopoietic stem cells produce a common myeloid progenitor and a common myelo-lymphoid progenitor (CMLP) that in their turn produce a bipotential myeloid T progenitor and a myeloid B progenitor (Kawamoto et al. [46]).

Although HSCs have the unique properties of leading to the formation of blood cells (self-renewal capacity and the ability to differentiate and produce different blood cells), they are dependent on several other factors: the environment that enables cell-to-cell and cell-to-matrix interaction (micromedia, niche), cytokines (growth, proliferation, differentiation and maturation factors) and humoral feedback from peripheral target tissue (Howard et al. [43], Young [44], Kaushansky et al. [45], Abkowitz [47], Cucuianu and Precup [13]). These unique properties are dependent upon two major groups of control factors: intrinsic cell factors (DNA alterations) and extrinsic cell factors (microenvironment factors, humoral feedback, cytokines) (Ramalingam et al. [48] and Zon [49]).

Even though mammals, humans included, have a stock of only $2\times {10}^4$ HSCs (Abkowitz [47]), they can give birth and release into the blood stream approximately 2.5 billion erythrocytes/kg/day, 2.5 billion platelets/kg/day and 1 billion granulocytes/kg/day.

Leukemias are an heterogeneous group of malignant disorders, also known as cancer, arising from one mutant hematopoietic stem cell (mHSC) (Howard et al. [43], Jilkine and Gutenkunst [35] Driessens et al. [36], Klein et al. [37], Lopez-Garcia et al. [38], and Snippert et al. [39]). In this study, we analyze the dynamics of HSCs and mHSCs by assuming that at least one mHSC can be found in the human body. Therefore, the complex biological processes on which the hematopoiesis is based are not completely involved.

There are four main types of leukemia (based on their progression-chronic or acute, and on the type of the affected cell-myeloid or lymphoid) (Neiman [7]). Although current guidelines include a more comprehensive and detailed classification of leukemias (with types and subtypes of cells, mutations acquired, etc.), they do not make the subject of the current study, nor do they bring useful information to our mathematical model (Arber et al. [50]).

The mHSCs have an abnormal process of differentiation and particular characteristics compared to HSCs, due to their acquired genetic and epigenetic abnormalities: increased proliferation/growth advantage; lower sensitivity to apoptosis and to the environment; poor differentiation; squeezing out of normal HSCs from bone marrow (Roeder and d’Inverno [51]). All these can lead to a wide range of clinical results with major impacts on patients’ health.

Chronic Myeloid Leukemia (CML) is an acquired myeloproliferative disorder (Howard et al. [43] and Neiman [7]). The CML is probably the first recognized leukemia, dating back to the 1840s (Young [44]).

CML represents 15% of all types of leukemia, annually occurring in 2 out of every 100,000 men and 1.1 out of every 100,000 women (Hemminki and Jiang [52]). Most frequently, the diagnosis is established during routine blood tests, occurring more frequently after the fifth decade of life (Howard et al. [43]). The signs and symptoms may include anemia, splenomegaly, weight loss, and dyspnea on exertion. Among the risk factors, ionizing radiation was observed to have a role in developing CML (Young [44]).

The hallmark of CML is the Philadelphia chromosome (Ph). Ph is characterized by a mutation in the normal hematopoietic stem cell population (Howard et al. [43]), generated by one abnormal stem cell, with the t(9;22)(q34;q11) (Young [44]) mutation, a reciprocal translocation of the ABL gene from chromosome 9 to chromosome 22, next to the BCR gene. This newborn BCR-ABL gene codes a fusion protein with tyrosine kinase activity, which apparently influences whether the cell lives or dies, and proliferates or not.

This type of leukemia typically undergoes three phases: chronic (also called indolent) phase, accelerated or transitory phase, and acute/blast phase (Abkowitz [47] and Arber et al. [50]). Due to the difficulty of separating the accelerated phase from the blast phase, we will refer to them as a whole, naming it the accelerated-acute phase. Most cases are diagnosed during the chronic phase and rarely during one of the other two phases. Once the mutation has occurred in one HSC, it starts a series of divisions, followed by differentiation and maturation that no longer obey the feedback and control mechanisms that apply to healthy HSCs. Therefore, the mutant cells divide at a quicker rate, producing a large number of thrombocytes and leukocytes, resulting in a population of cells where mHSCs are dominant.

At a certain point of CML, the occurrence of other events (most probably the acquisition of other genetic mutations) leads to an instability of the mHSC population, which consequently leads to an accelerated-acute phase, resulting in an exponential increase of the number of immature stem cells. Cells multiply in a more accelerated manner and do not undergo differentiation, resulting in a blast phase, similar to various types of acute leukemia (myeloid—70%, lymphoid—20%, and mixt type—10%) in terms of symptoms and clinical findings (Young [44], Abkowitz [47], and Neiman [7]). After reaching the accelerated phase or the blast phase, patients have a median survival of 4.8 years or 6 months, respectively (Kantarjian et al. [53]).

In terms of treatment, even though tyrosinkinase inhibitors are effective in 70%–80% of cases of CML, stem cell transplantation (SCT) seems to be the only curative treatment, involving, nevertheless, a high mortality rate due to complications (Howard et al. [43] and Thomas [54]).

2. The Mathematical Model

In this section, we give the mathematical model, we obtain the associated steady states, and we investigate their stability. Based on this analysis, we characterize the normal hematopoietic state, and the chronic and accelerated-acute stages in CML.

2.1. The Normal-Leukemic Dynamic System

The mathematical modeling of the time evolution of a population p of any nature begins in a first approximation (assuming no constraints exist) with the Malthusian equation

$$\frac{dp}{dt}=ap-cp,$$

where $p=p\left(t\right)$ is the population size at time $t$, and a and c are the growth and death (per capita) rates, respectively. Assuming that the growth rate is bigger than the death rate, the population will increase exponentially according to the law

$$p\left(t\right)=p\left(0\right)exp\left(\left(a-c\right)t\right),$$

which is non-realistic in the long run, particularly for limited biological populations. Thus, a more realistic approach is to consider that the growth (and/or death) rate will change during evolution by a self-limiting mechanism or exterior influences. For instance, the growth rate of a self-limiting population can be $a/\left(1+bp\right)$ depending on the population size itself. Here, b is a proportionality factor that shows how sensitive that population is with respect to its own size. This shows that the growth rate decreases as the population size p increases. In addition, the influence over p of a competitive population $q=q\left(t\right)$ can be simulated in the model by a growth rate of the form $a/\left(1+b_{1}p+b_{2}q\right)$, where the ratio $b_2/b_1$ shows how strong the diminishing effect due to population q is, compared to that of self-limitation.

Applied to the normal and abnormal stem cell populations denoted by x and $z$, the above modeling choice leads to the following differential system

$$\begin{array}{ccc}\hfill \frac{dx}{dt}& =& \frac{a}{1+{\beta }_{1}x+{\beta }_{2}z}x-cx\hfill \\ \hfill \frac{dz}{dt}& =& \frac{A}{1+{\gamma }_{1}x+{\gamma }_{2}z}z-Cz.\hfill \end{array}$$

(1)

Here, since abnormal cells have a stronger diminishing effect on the growth rate of normal cells than on their own growth rate, it is natural to assume that

$${\beta }_2>{\gamma }_2.$$

(2)

Additionally, the almost negligible effect of normal cells over the growth rate of population z (i.e., ${\gamma }_1$ is much smaller than ${\beta }_1$) justifies the inequality

$$\frac{{\gamma }_2}{{\gamma }_1}>\frac{{\beta }_2}{{\beta }_1}.$$

(3)

For the mathematical analysis of most models, it is often convenient that the number of parameters is reduced as much as possible. Thus, in our case, we can reduce the number of parameters ${\beta }_1,{\beta }_2,{\gamma }_1,{\gamma }_2$ to three by making the change of variable (equivalently, by rescaling the abnormal cell population)

$$y=\frac{{\gamma }_2}{{\gamma }_1}z.$$

Substituting in (1) yields the system

$$\begin{array}{ccc}\hfill \frac{dx}{dt}& =& \frac{a}{1+b_{1}x+b_{2}y}x-cx\hfill \end{array}$$

(4)

$$\begin{array}{ccc}\hfill \frac{dy}{dt}& =& \frac{A}{1+B(x+y)}y-Cy.\hfill \end{array}$$

(5)

where

$$b_1={\beta }_1,b_2={\beta }_2\frac{{\gamma }_1}{{\gamma }_2},B={\gamma }_1.$$

Then (3) yields $b_1>b_2,$ which shows the different contributions of the normal and abnormal stem cells to the diminution of the nonrestrictive growth rate of the normal cell population, while (2) guarantees that $b_2>B,$ reflecting the eventual advantage of the abnormal cells of being less sensitive to the bone marrow microenvironment.

The system (4)–(5) is our basic model for normal-leukemic cell evolution. Here, the model parameters a and A are the nonrestrictive growth rates (due to self-renewal) of normal and abnormal stem cells, respectively; $b_1,b_2$ and B are the bone marrow microenvironment sensitivities; while c and C stand for their cell death rates (due to differentiation, apoptosis, and other elimination mechanisms) (see Alenzi et al. [55], Cisneros et al. [56], Domen [57], Riether et al. [58], Vivier et al. [59]). The terms

$$\frac{1}{1+b_{1}x+b_{2}y}\mathrm{and}\frac{1}{1+B(x+y)}$$

quantify the impact of the induced crowding in the bone marrow microenvironment, introduce competition between normal and abnormal cells, and guarantee homeostasis at the level of cell population. We assume that for both cell populations, the growth rate is greater than the death rate, that is,

$$a>c\mathrm{and}A>C.$$

Note that an alternative model for normal-abnormal dynamics can be found in Neiman [7], where the role of the parameter $b_1/b_2$ is given by a parameter denoted by g and assumed greater than or equal to one.

The case $b_1=b_2$ was considered by Dingli and Michor [11] and Cucuianu and Precup [13]. In this case there are only two non-zero steady states, $(d,0)$ and $(0,D)$, where d and D represent the homeostatic amounts of normal and abnormal stem cells, and they are given by

$$d=\frac{1}{b_1}\left(\frac{a}{c}-1\right)andD=\frac{1}{B}\left(\frac{A}{C}-1\right).$$

(6)

In this paper we assume that $b_1>b_2$. As we shall see, in this case, besides the non-zero steady states $(d,0)$ and $(0,D)$, a steady state $(x^{*},y^{*})$ could also exist with both positive components, i.e., $x^{*}>0$ and $y^{*}>0$. This makes the new model able to differentiate between chronic and accelerated-acute phases in CML.

We continue to analyze the system (4)–(5).

(a) Monotonicity of the solutions. The function $x\left(t\right)$ increases during the time intervals where $dx/dt>0$, i.e., $a/(1+b_{1}x\left(t\right)+b_{2}y\left(t\right))-c>0$, or equivalently $x\left(t\right)+(b_2/b_1)y\left(t\right)<d$. Hence,

$$\begin{array}{cc}x\left(t\right)\mathrm{increases}\mathrm{as}\mathrm{long}\mathrm{as}\hfill & x\left(t\right)+\frac{b_2}{b_1}y\left(t\right)<d,\hfill \\ x\left(t\right)\mathrm{decreases}\mathrm{as}\mathrm{long}\mathrm{as}\hfill & x\left(t\right)+\frac{b_2}{b_1}y\left(t\right)>d.\hfill \end{array}$$

From a biological point of view, the increasing or decreasing of the normal cell population occurs on those intervals of time where the weighted total of cells remains under or upper the normal homeostatic level, respectively.

Similarly,

$$\begin{array}{cc}y\left(t\right)\mathrm{increases}\mathrm{as}\mathrm{long}\mathrm{as}\hfill & x\left(t\right)+y\left(t\right)<D,\hfill \\ y\left(t\right)\mathrm{decreases}\mathrm{as}\mathrm{long}\mathrm{as}\hfill & x\left(t\right)+y\left(t\right)>D.\hfill \end{array}$$

(b) Steady states. A steady state (or an equilibrium) is a constant solution, i.e., a solution for which $dx/dt=dy/dt=0$. Hence, the steady states are obtained by solving the algebraic system

$$\begin{array}{ccc}\hfill \frac{a}{1+b_{1}x+b_{2}y}x-cx& =& 0\hfill \end{array}$$

(7a)

$$\begin{array}{ccc}\hfill \frac{A}{1+B(x+y)}y-Cy& =& 0.\hfill \end{array}$$

(7b)

The solutions of the system (7a)–(7b) are the pairs

$$(0,0),(d,0),(0,D)\mathrm{and}(x^{*},y^{*}),$$

where $d,$D are given by (6),

$$x^{*}=-\frac{b_{2}c(A-C)-BC(a-c)}{BCc(b_1-b_2)}$$

and

$$y^{*}=\frac{b_{1}c(A-C)-BC(a-c)}{BCc(b_1-b_2)}.$$

Direct calculation leads to

$$x^{*}=\frac{b_1}{b_1-b_2}d-\frac{b_2}{b_1-b_2}D,y^{*}=\frac{b_1}{b_1-b_2}\left(D-d\right).$$

It is easy to see that under the assumption that $b_1>b_2$, both numbers $x^{*}$ and $y^{*}$ are positive (acceptable values from a biological point of view) if and only if

$$d<D<\frac{b_1}{b_2}d.$$

Therefore, in addition to the non-zero steady states $\left(d,0\right)$ and $\left(0,D\right)$, a positive steady state $\left(x^{*},y^{*}\right)$ appears, contrary to the case where $b_1=b_2$.

(c) Stability. We study the stability of the steady states of the system (4)–(5) using the standard first approximation method (for details see Kaplan and Glass [60], Coddington and Levinson [61] and Jones et al. [62]). According to this method, an equilibrium $(\alpha ,\beta )$ is asymptotically stable if the Jacobian matrix $J(\alpha ,\beta )$ is a Hurwitz matrix, i.e., Re $\lambda <0$ for all its characteristic roots $\lambda$, and is unstable if Re $\lambda >0$ for at least one of its characteristic roots.

It should be emphasized that the local stability that we analyze below coincides with the global stability of the corresponding equilibria. This guarantees that the predictions are robust to uncertainties in the initial cell counts.

Note that the steady state $(0,0)$ is unstable as can be shown based on the assumptions $a>c$ and $A>C$.

For the steady state $(d,0)$, the eigenvalues of the Jacobian matrix $J(d,0)$ are

$${\lambda }_1=-\frac{c(a-c)}{a},{\lambda }_2=\frac{b_{1}c(A-C)-BC(a-c)}{b_{1}c+B(a-c)}.$$

Obviously ${\lambda }_1<0$. Thus, the steady state $(d,0)$ is asymptotically stable if and only if ${\lambda }_2<0$, or equivalently, if $D<d$. On the contrary, if $D>d$, then the equilibrium $(d,0)$ is unstable.

The eigenvalues of the Jacobian matrix $J(0,D)$ associated to the equilibrium $\left(0,D\right)$ are

$${\lambda }_1=-\frac{b_{2}c(A-C)-BC(a-c)}{b_2(A-C)+BC},{\lambda }_2=-\frac{C(A-C)}{A}.$$

Here ${\lambda }_2<0$, and so the steady state $(0,D)$ is asymptotically stable if and only if ${\lambda }_1<0$. This happens if $D>(b_1/b_2)d$. Contrarily, if $D<(b_1/b_2)d$, then the steady state $(0,D)$ is unstable.

In the chronic case $d<D<(b_1/b_2)d$, by direct computation, we easily find the following representation of the Jacobian matrix $J(x^{*},y^{*})$

$$J\left(x^{*},y^{*}\right)=\left(\begin{array}{cc}-\frac{b_1{c}^2}{a}{x}^{*}& -\frac{b_2{c}^2}{a}{x}^{*}\\ -\frac{B{C}^2}{A}{y}^{*}& -\frac{B{C}^2}{A}{y}^{*}\end{array}\right)$$

whose characteristic polynomial is

$${\lambda }^2-\lambda \mathrm{tr}J\left(x^{*},y^{*}\right)+\det J\left(x^{*},y^{*}\right),$$

where

$$\mathrm{tr}J\left(x^{*},y^{*}\right)=-\frac{b_1{c}^2}{a}{x}^{*}-\frac{B{C}^2}{A}{y}^{*},\det J\left(x^{*},y^{*}\right)=\frac{B{c}^2{C}^2}{aA}\left(b_1-b_2\right)x^{*}{y}^{*}.$$

Obviously $\mathrm{tr}J\left(x^{*},y^{*}\right)<0$ and $\det J\left(x^{*},y^{*}\right)>0$. These imply that $J\left(x^{*},y^{*}\right)$ is a Hurwitz matrix. Hence the equilibrium $\left(x^{*},y^{*}\right)$ is asymptotically stable.

The above discussion shows a qualitative change of the system’s behavior, i.e., a change of stable equilibrium where the parameter D varies. The values of this parameter at which the stable equilibrium changes (called bifurcation points) are $D=d$ and $D=(b_1/b_2)d$. We can summarize the bifurcation analysis of our system as follows (see Figure 1):

• If $D<d$, then the steady state $(d,0)$ is asymptotically stable, and the steady state $(0,D)$ is unstable.
• If $d<D<(b_1/b_2)d$, then the steady state $(x^{*},y^{*})$ is positive and asymptotically stable, and the steady states $(d,0)$ and $(0,D)$ are unstable.
• If $D>(b_1/b_2)d$, then the steady state $(0,D)$ is asymptotically stable, and the steady state $(d,0)$ is unstable.

From a biological point of view, as long as normal cells have an advantage over abnormal cells $\left(d>D\right)$, the healthy state may not present any perceptible changes. A balance between normal and abnormal cells will exist as long as abnormal cells have a not too large advantage over normal cells $\left(d<D<\left(b_1/b_2\right)d\right)$, and the balance is disturbed, leading to an accelerated-acute phase of the disease, once the advantage of leukemic cells becomes significant $\left(D>\left(b_1/b_2\right)d\right)$ (see Stine and Matunis [63]).

2.2. The Mathematical–Biological Interpretation

In view of the above discussion, we can claim that the relationship $D<d$ characterizes the normal hematopoiesis; the relationship $d<D<(b_1/b_2)d$ stands for the chronic phase of leukemia; meanwhile, the inequality $D>(b_1/b_2)d$ characterizes the accelerated-acute phase of the disease. Indeed (see Figure 2), in the case where $D<d$, the normal cell population $x\left(t\right)$ approaches the equilibrium abundance d (normal homeostatic state) while the abnormal cell population $y\left(t\right)$ tends in time to zero; in the case where $d<D<(b_1/b_2)d$, the cell populations $x\left(t\right)$ and $y\left(t\right)$ approach their equilibrium abundances $x^{*}$ and $y^{*}$, respectively; finally, if $D>(b_1/b_2)d$, then the leukemic cell population becomes dominant approaching its equilibrium abundance D (leukemic homeostatic state) and leads to the elimination of the normal cells, i.e., $x\left(t\right)$ tends to zero.

Finally, a condition like $D=d$ or $D=(b_1/b_2)d$ is physiologically very unstable, since small variations of the kinetic parameters can switch the normal state into the chronic leukemic state and vice-versa, if $D=d$, and can switch the chronic state into the accelerated-acute phase and vice-versa, if $D=(b_1/b_2)d$. Additionally, from a medical point of view, the situations $D=d$ and $D=(b_1/b_2)d$ are practically undetectable.

In terms of the system’s biological growth parameters, the hematological states are characterized as follows:

$$\begin{array}{cc}\frac{1}{B}\left(\frac{A}{C}-1\right)<\frac{1}{b_1}\left(\frac{a}{c}-1\right)\hfill & \left(\mathrm{normal}\mathrm{state}\right);\hfill \\ \frac{1}{b_1}\left(\frac{a}{c}-1\right)<\frac{1}{B}\left(\frac{A}{C}-1\right)<\frac{1}{b_2}\left(\frac{a}{c}-1\right)\hfill & \left(\mathrm{chronic}\mathrm{phase}\right);\hfill \\ \frac{1}{b_2}\left(\frac{a}{c}-1\right)<\frac{1}{B}\left(\frac{A}{C}-1\right)\hfill & \left(\mathrm{accelerated}-\mathrm{acute}\mathrm{phase}\right).\hfill \end{array}$$

A diagram of the transition from normal hematopoiesis to chronic and accelerated-acute stages in CML is presented in Figure 1. Notice that the length of the interval $\left[d,\left(b_1/b_2\right)d\right]$ that corresponds to the chronic phase is $\left(b_1/b_2-1\right)d$ and depends on the ratio $b_1/b_2$. The bigger the ratio $b_1/b_2$ is, the larger the interval in which a patient’s parameter D can lay in chronic phase. Notice that the patient-related parameter $b_1/b_2$ gives the ratio of contributions of the normal and abnormal stem cells to the diminution of the nonrestrictive growth rate of the normal cell population, aiming to restore homeostasis.

According to our model, values of parameter D close to d correspond to early stages of the disease, while values of D close to $\left(b_1/b_2\right)d$ indicate advanced stages of the disease moving towards the accelerated-acute phase.

3. Numerical Simulation of the Model

Further on, we illustrate the theoretical results on the system (4)–(5) by numerical simulations using the Maple package.

3.1. Parameter Estimations

The parameters employed by our model depend on a large number of biophysical and biochemical mechanisms. The latter make the exact estimation of these parameters almost impossible. Instead, one may expect that from any estimation procedure with confidence intervals, parameters are obtained. For a qualitative analysis such as ours, parameter estimation is not essential, and as already seen, relationships between parameters are enough. Parameter estimation becomes essential when the model is applied for real-time predictions and individual patients.

According to the paper of Dingli and Michor [11], the number of the stem cells in a healthy adult body is approximately $d=2\times {10}^4$ (normal homeostatic state), and the growth and death rates of normal stem cells could be taken $a=0.005$ and $c=0.002$. However, recent studies have shown different growth rates of HSC. One study concluded that HSCs divide on average every 40 weeks, with a range from 25 to 50 weeks (see Catlin et al. [64]). Others even suggested the existence of two types of HSCs with different replication rates (dormant HSCs divide about every 145 days, active HSCs divide about every 36 days) (see Wilson et al. [65]). Regarding the life span of HSC, recent studies have suggested a death rate ranging from 10 to 60 months (see Sieburg et al. [66] and Sieburg et al. [67]). The factor $b_1/b_2$ allows for the possibility that the abnormal stem cells are less sensitive to the bone marrow microenvironment than the normal cells. The parameter $b_1$ that stands for the bone marrow microenvironment sensitivity of the normal stem cells can be estimated from the expression of d, namely, $b_1=\left(a/c-1\right)/d=0.75\times {10}^{-4}$.

Parameters $A,B,C$ and $b_2$ vary from patient to patient, and so do parameter D and the length of the chronicity interval $\left[d,\left(b_1/b_2\right)d\right]$. For our numerical simulations we choose value 2 for $b_1/b_2$ and we then have $b_2\simeq 0.38\times {10}^{-4}$. Additionally, as in Dingli and Michor [11], we assume that the value of parameter B is approximately half of $b_2$; hence, $B\simeq 0.19\times {10}^{-4}$. As regards the parameters A and $C$, several values are considered in our simulations such that all the previous relationships between the model parameters hold.

3.2. Numerical Simulations

We simulate numerically system (4)–(5) in order to investigate the behavior of normal and abnormal stem cell populations in each of the cases: $D<d$ (normal state); $d<D<(b_1/b_2)d$ (chronic state); and $(b_1/b_2)d<D$ (accelerated-acute state). The graphs of $x\left(t\right)$ (blue solid line) and $y\left(t\right)$ (red broken line) for a time interval $0\le t\le T$, are represented for different sets of values of the model parameters: $a,b_1,b_2,c,A,B,C$; initial values $x\left(0\right)$, $y\left(0\right)$; and length T of the time interval.

We shall restrict our simulations to the situations presented in

Table 1. The numerical simulation cases. $a,A$ = growth rates; $b_1,b_2,B$ = bone marrow microenvironment sensitivity; $c,C$ = death rates; $a,b_1,b_2,c$ = normal stem cell parameters; and $A,B,C$ = abnormal (leukemic) stem cell parameters.

Case I	Case II	Case III	Case IV
$a<A$	$a<A$	$a>A$	$a>A$
$c<C$	$c>C$	$c<C$	$c>C$
$b_1>b_2>B$	$b_1>b_2>B$	$b_1>b_2>B$	$b_1>b_2>B$

. In all of the cases we assume that the abnormal stem cells are less sensitive to environmental crowding than the normal stem cells, that is, $b_1>b_2>B$.

Case I: In this case the growth and death rates of normal stem cells are smaller than the growth and death rates of abnormal stem cells. Figure 3a shows the behavior in time ($T=3000$ days) of the two cell populations for the parameter values provided in the first line in

Table 2. Parameter values for simulations. S-S = steady state; $d=2\times {10}^4$ (normal); D = variable parameter (leukemic).

Figure	a	${\mathit{b}}_1\times {10}^{-4}$	${\mathit{b}}_2\times {10}^{-4}$	c	A	$\mathit{B}\times {10}^{-4}$	C	$\mathit{S}-\mathit{S}$
Figure 3a	$0.005$	$0.75$	$0.38$	$0.002$	$0.01$	$0.19$	$0.009$	$(d,0)$
Figure 3b	$0.005$	$0.75$	$0.38$	$0.002$	$0.01$	$0.19$	$0.007$	$(x^{*},y^{*})$
Figure 3c	$0.005$	$0.75$	$0.38$	$0.002$	$0.01$	$0.19$	$0.004$	$(0,D)$
Figure 4	$0.005$	$0.75$	$0.38$	$0.002$	$0.007$	$0.19$	$0.001$	$(0,D)$
Figure 5a	$0.005$	$0.75$	$0.38$	$0.002$	$0.004$	$0.19$	$0.003$	$(d,0)$
Figure 5b	$0.005$	$0.75$	$0.38$	$0.002$	$0.0045$	$0.19$	$0.003$	$(x^{*},y^{*})$
Figure 5c	$0.005$	$0.75$	$0.38$	$0.002$	$0.0045$	$0.19$	$0.0025$	$(0,D)$
Figure 6a	$0.005$	$0.75$	$0.38$	$0.002$	$0.0012$	$0.19$	$0.001$	$(d,0)$
Figure 6b	$0.005$	$0.75$	$0.38$	$0.002$	$0.0015$	$0.19$	$0.001$	$(x^{*},y^{*})$
Figure 6c	$0.005$	$0.75$	$0.38$	$0.002$	$0.0025$	$0.19$	$0.001$	$(0,D)$

, values that correspond to the normal hematopoietic state ($D<d$). The normal stem cell population $x\left(t\right)$ (blue solid line) tends to the value d while the abnormal stem cell population $y\left(t\right)$ (red broken line) tends towards 0. Biologically, this mutant/abnormal cell extinction, due to random events, has been explained and demonstrated in several studies using stem cell lineages (see Jilkine and Gutenkunst [35], Driessens et al. [36], Klein et al. [37], Lopez-Garcia et al. [38], and Snippert et al. [39]). Figure 3b shows the behavior in time ($T=$ 25,000 days) of the two cell populations for the parameter values from the second line of Table 2, values that correspond to the chronic state $d<D<(b_1/b_2)d$. The normal and leukemic stem cell populations $x\left(t\right)$, $y\left(t\right)$ tend toward $x^{*}$ and $y^{*}$, respectively. Figure 3c shows the behavior in time ($T=8000$ days) of the two cell populations for the corresponding parameter values from Table 2, values that lead to the accelerated-acute state $D>(b_1/b_2)d$. The normal stem cell population $x\left(t\right)$ tends towards $0$, while the leukemic stem cell population $y\left(t\right)$ tends to the value $D$.

Case II: Here the growth rate of normal stem cells is smaller than the growth rate of abnormal stem cells, and the death rate of normal stem cells is greater than the death rate of abnormal stem cells. Then $A/C>a/c$, and since $1/B>1/b_2$, we immediately see that $\left(1/B\right)\left(A/C-1\right)>\left(1/b_2\right)\left(a/c-1\right)$, or equivalently $D>\left(b_1/b_2\right)d$. Hence, in this case only the accelerated acute state is possible. Figure 4 shows the behavior in time ($T=6000$ days) of the two cell populations for the corresponding parameter values from Table 2. The normal stem cell population $x\left(t\right)$ tends towards $0$, while the leukemic stem cell population $y\left(t\right)$ tends to the value $D$.

Case III: In this case, the growth rate of normal stem cells is greater than the growth rate of abnormal stem cells, and the death rate of normal stem cells is smaller than the death rate of abnormal stem cells. Figure 5a shows the behavior in time ($T=$ 25,000 days) of the two cell populations for the corresponding parameter values from Table 2, values that correspond to the normal hematopoietic state $D<d$. The normal stem cell population $x\left(t\right)$ becomes arbitrarily close to the value $d$, while the abnormal stem cell population $y\left(t\right)$ tends towards $0$. Figure 5b shows the behavior in time ($T=$ 25,000 days) of the two cell populations for the corresponding parameter values from Table 2, values that correspond to the chronic state $d<D<(b_1/b_2)d$. The normal and leukemic stem cell populations $x\left(t\right)$, $y\left(t\right)$ tend toward $x^{*}$ and $y^{*}$, respectively. Figure 5c shows the behavior in time ($T=$ 40,000 days) of the two cell populations for the corresponding parameter values from Table 2, values that lead to the accelerated acute state $D>(b_1/b_2)d$. The normal stem cell population $x\left(t\right)$ tends towards $0$, while the leukemic stem cell population $y\left(t\right)$ tends to the value $D$.

Case IV: In this case, the growth rate of normal stem cells is greater than the growth rate of abnormal stem cells, and the death rate of normal stem cells is greater than the death rate of abnormal stem cells. Figure 6a shows the behavior in time ($T=$ 25,000 days) of the two cell populations for the corresponding parameter values from Table 2, values that correspond to the normal hematopoietic state $D<d$. Normal stem cell population $x\left(t\right)$ tends to the value $d$, and abnormal stem cell population $y\left(t\right)$ approaches $0$. Figure 6b shows the behavior in time ($T=$ 25,000 days) of the two cell populations for the corresponding parameter values from Table 2, values that correspond to the chronic state $d<D<(b_1/b_2)d$. The normal and leukemic stem cell populations $x\left(t\right)$, $y\left(t\right)$ tend toward $x^{*}$ and $y^{*}$, respectively. Figure 6c shows the behavior in time ($T=$ 25,000 days) of the two cell populations for the corresponding parameter values from Table 2, values that lead to the accelerated-acute state $D>(b_1/b_2)d$. The normal stem cell population $x\left(t\right)$ approaches $0$, while the leukemic stem cell population $y\left(t\right)$ tends to the value $D$.

4. The Model Extended to Terminally Differentiated Cells

Working at the level of primitive stem cells, there is not an ordinary method to determine the sizes of cell populations. Therefore, it would be useful to have a reflection of the primitive stem cell evolution at the level of terminally differentiated cells, since the latter can be easily estimated by current blood tests.

The idea appears in Michor et al. [42], and when applied to our model yields the extended system

$$\begin{array}{cccc}\frac{d{x}_1}{dt}=\frac{a_1}{1+b_1{x}_1+b_2{y}_1}{x}_1-c_1{x}_1\hfill & \left(\mathrm{NSC}\right)\hfill & \frac{d{y}_1}{dt}=\frac{A_1}{1+B(x_1+y_1)}{y}_1-C_1{y}_1\hfill & \left(\mathrm{ASC}\right)\hfill \\ \frac{d{x}_2}{dt}=a_2{x}_1-c_2{x}_2\hfill & \left(\mathrm{NPC}\right)\hfill & \frac{d{y}_2}{dt}=A_2{y}_1-C_2{y}_2\hfill & \left(\mathrm{APC}\right)\hfill \\ \frac{d{x}_3}{dt}=a_3{x}_2-c_3{x}_3\hfill & \left(\mathrm{NDC}\right)\hfill & \frac{d{y}_3}{dt}=A_3{y}_2-C_3{y}_3\hfill & \left(\mathrm{ADC}\right)\hfill \\ \frac{d{x}_4}{dt}=a_4{x}_3-c_4{x}_4\hfill & \left(\mathrm{NTC}\right)\hfill & \frac{d{y}_4}{dt}=A_4{y}_3-C_4{y}_4\hfill & \left(\mathrm{ATC}\right).\hfill \end{array}$$

Here, $x_2\left(t\right)$, $y_2\left(t\right)$ stand for the normal (N) and abnormal (A) progenitor cell (PC) populations; $x_3\left(t\right)$, $y_3\left(t\right)$ stand for the normal and abnormal differentiated cell (DC) populations; and $x_4\left(t\right)$, $y_4\left(t\right)$ stand for the normal and abnormal terminally differentiated cell (TC) populations, respectively.

Notice the different form of the additional six equations associated to the succeeding cell compartments. It is the balance expression of the transfer from one compartment to the next; namely, it states that the rate of change of cell population equals the amount of new cells produced by the anterior compartment minus the amount of cells that leave (by differentiation or apoptosis) that compartment. Thus, the new parameters $a_2$, $A_2$ are the rates at which normal and abnormal progenitor cells are produced from normal and abnormal stem cells; $a_3$, $A_3$ are the rates at which normal and abnormal differentiated cells are produced from normal and abnormal progenitor cells; $a_4$, $A_4$ are the rates at which normal and abnormal terminally differentiated cells are produced from normal and abnormal differentiated cells; and $c_2$, $c_3$, $c_4$, $C_2$, $C_3$, $C_4$ are the death rates of normal and abnormal progenitors, and differentiated and terminally differentiated cells.

In the equilibrium state, we assume that in a healthy adult body the number of stem cells is $d=x_1^{*}=2\times {10}^5$, the number of progenitor cells is $x_2^{*}=1\times {10}^8$, the number of differentiated cells is $x_3^{*}=1\times {10}^{10}$, and the number of terminally differentiated cells is $x_4^{*}=1\times {10}^{12}$ (see Michor et al. [42]). Consequently, in the equilibrium state, we have for progenitor cells $a_2{x}_1^{*}-c_2{x}_2^{*}=0$, whence $a_2/c_2=x_2^{*}/x_1^{*}=5\times {10}^2$, for differentiated cells $a_3{x}_2^{*}-c_3{x}_3^{*}=0$, whence $a_3/c_3=x_3^{*}/x_2^{*}=1\times {10}^2$, and for terminally differentiated cells $a_4{x}_3^{*}-c_4{x}_4^{*}=0$; hence $a_4/c_4=x_4^{*}/x_3^{*}=1\times {10}^2$. Therefore, if $c_2=0.008$, $c_3=0.05$ and $c_4=1$ (see Michor et al. [42]), then $a_2=4$, $a_3=5$, and $a_4=100$.

Note that if

$$(x_{1E},y_{1E})$$

is any equilibrium (E) of the initial system (4)–(5), then

$$\begin{array}{c}(x_{1E},y_{1E},a_2{x}_{1E}/c_2,A_2{y}_{1E}/C_2,a_2{a}_3{x}_{1E}/c_2{c}_3,A_2{A}_3{y}_{1E}/C_2{C}_3,\hfill \\ a_2{a}_3{a}_4{x}_{1E}/c_2{c}_3{c}_4,A_2{A}_3{A}_4{y}_{1E}/C_2{C}_3{C}_4)\hfill \end{array}$$

is an equilibrium of the extended system, and the two equilibria have the same stability property. Therefore, working at the level of primitive stem cells is equivalent to working at the level of any one of the succeeding classes of cells. Thus, if by blood tests one estimates the steady state ration $x_{4E}/y_{4E}$ between healthy and unhealthy terminally differentiated cells as being equal to $\lambda$, then we can immediately calculate the analogue steady state ratios of differentiated, progenitor, and stem cells, as follows:

$$\frac{x_{3E}}{y_{3E}}=\lambda \frac{A_4{c}_4}{a_4{C}_4},\frac{x_{2E}}{y_{2E}}=\lambda \frac{A_3{A}_4{c}_3{c}_4}{a_3{a}_4{C}_3{C}_4},\frac{x_{1E}}{y_{1E}}=\lambda \frac{A_2{A}_3{A}_4{c}_2{c}_3{c}_4}{a_2{a}_3{a}_4{C}_2{C}_3{C}_4}.$$

Figure 7a–d show that there is indeed a parallelism between the behaviors of normal and abnormal cell populations in all four cell compartments.

5. Discussion and Conclusions

By resuming and adding a new parameter to the system created by Dingli and Michor, the present study tried to bring the model and the numerical simulations one step closer to the complex reality of leukemic stem cell evolution, while also maintaining the model simple enough to draw conclusions.

We were able to analyze elements and characteristics of both normal and mutant stem cells, unproven or unnoticed on the original model, but hopefully of some use for the upcoming research and mathematical modeling of leukemic pathology.

By analyzing the numerical simulation for the first case, in which leukemic stem cells had a rate of multiplication higher than that of normal stem cells (in our case two times higher), we observed that the evolution of leukemic pathology was dependent on the ratio given by death rates of the two types of stem cells lines. Therefore, at a mHSC death rate four times higher than that of normal HSC, the disease progressed and led to the disappearance of mHSCs, while in case of a lower death rate of the leukemic stem cells ranging between $4\times c$ and $3\times c$, the disease evolved towards a stabilization of the two stem cell lines; that is, to the chronic phase of CML. Finally, the disease evolution shifted towards the accelerated-acute phase when the mHSC death rate decreased below $3\times c$, a situation clinically similar to acute myeloid leukemia (AML). The analysis of the second case reveals to us a very logical conclusion. The leukemic stem cells, having all the advantages (proliferation rate higher than normal stem cells, death rate lower than normal stem cells, and lower sensitivity towards the microenvironment) apparently leads us to only one possible case, that of the accelerated-acute phase of CML. This allowed us to confirm our previously mentioned conclusion, which stated that the disease evolved towards the specific course of the accelerated-acute phase of chronic myeloid leukemia in the absence of a death rate of mHSC several times higher than the death rate of normal stem cells. In terms of medical treatment, this case appears to underline the importance of microenvironment in the evolution of the leukemic disease. Furthermore, assuming that medical treatments fail to decrease the proliferation rate or increase the death rate of mHSC, the microenvironment seems to be the only therapeutic target left available. Additionally, we noticed that the accelerated-acute phase of CML developed approximately 6 years after the occurrence of the first leukemic stem cell (if the initial parameters were maintained constant).

A special situation is represented by the analysis of Cases III and IV, when the leukemic stem cells multiplication rate is lower than that of normal stem cells. Although the advantage gained by leukemic stem cells proliferation and growth rate is mentioned by many medical sources (Howard et al. [43], Young [44], Kaushansky et al. [45], Abkowitz [47], Cucuianu and Precup [13], Hemminki and Jiang [52], and Thomas [54]), the analysis of numerical simulations for our model showed the possibility of chronic myeloid leukemia occurring and shifting between its three stages. Therefore, in Cases III and IV, we observed the occurrence of all three stages, the disease following a similar course to that seen in Case I, but characterized by a longer evolutionary time frame, probably due to the slowness of the mHSCs proliferation rate. In clinical practice, though rare, such cases can be seen in elderly people that present a very slow evolution of the disease.

The analysis of the mathematical model and numerical simulations led us to clues and conclusions that are otherwise difficult to notice, explain, or even measure, relying solely on the gross figures produced by laboratory studies or clinical observation. According to the analysis of our mathematical model and numerical simulations, we can state that:

• The mHSCs proliferation rate is a predictive factor for the development of the accelerated-acute state: an increased rate of proliferation of these cells in comparison to normal stem cells determines the accelerated-acute phase to occur earlier;
• The death rate of leukemic stem cells is predictive for the global evolution of the disease, influencing the shifts between the different phases of the chronic myeloid leukemia.

Therefore, the clinical judgment, treatment plan, and research to improve therapy for leukemic disease could be influenced or based on the importance of these two factors. In terms of treatment, according to our mathematical model, we should probably focus more on controlling values of these parameters of mHSC in order to reach a coexisting phase of the two populations, taking into account patients’ symptoms and quality of life, rather than being aggressive and trying to eradicate all leukemic cells. In light of recent research and opinions (Enriquez-Navas et al. [68] and Gerlinger et al. [69]), trying to live and collaborate with cancer could be a more intelligent strategy than trying to eradicate it.

Certainly, real clinical experience regarding the shift of leukemic disease from one phase to another is more complex and may comprise other parameters that have not been taken into account by this study, either because they were impossible to include in our mathematical model or because they are still unknown by medical researchers. However, mathematical models and numerical simulations that are applied in the biomedical field can reveal aspects and ideas that deserve to be closely investigated in correlation with the medical practice and research. Furthermore, these interdisciplinary collaborations may be considered groundwork for the necessary reasoning process that will consequently identify research directions for improved treatments of extensively investigated pathologies, such as leukemia.