Analysis of mathematical model of leukemia

In this paper, a model describing the dynamic of chronic myeloid leukemia is studied. By analyzing the corresponding characteristic equations, the local stability of trivial and nontrivial equilibria are discussed. By establishing appropriate Lyapunov functions, we prove the global stability of the positive constant equilibrium solutions.


Introduction
Several recent mathematical models have been developed to study the dynamics of chronic myeloid leukemia (CML) under IMATINIB treatment, see ( [2,3,5,6]). In all of these studies, the authors conclude that IMATINIB does not completely eliminate leukemic cells, and propose that IMATINIB therapy should be combined with an additional form of treatment.
In [5], the authors analyze a four-compartment differential equation model based on the known biology of the hematopoietic system to investigate the dynamics of CML with IMATINIB treatment. In their model, Michor et al. [5] study the interaction between leukemic cells and IMATINIB, they assume that these cells differentiate through four stages of their life cycle, beginning with leukemia stem cells. IMATINIB reduces the rate at which leukemic cells pass from one stage to the next, causing a rapid drop in the leukemia population. Based on their assumptions and analysis, they propose that leukemia inevitably persists, because IMATINIB hinders the differentiation of differentiated leukemic cells, but does not affect the leukemia stem cells. In particular, Michor et al. [5] hypothesizes that there is always a steadily growing population of leukemia stem cells despite IMATINIB treatment. As a result, based on their model, the leukemia population under IMATINIB eventually relapses, regardless of whether the model considers IMATINIB resistance mutations.
In [2], the authors incorporate the anti-leukemia immune response in CML patients on IMATINIB therapy to the model proposed in [5] by adding interactions with anti-leukemia T-cells. They formulate their mathematical model as a system of delay differential equations, and prove that the immune ITM Web of Conferences response may play a critical role in determining the length of time that CML patients under IMATINIB treatment remain in remission.
In [6], Roeder et al. develop a similar model of CML and IMATINIB. However, they subdivide the leukemia stem cells into two compartments: proliferating and quiescent cells. Proliferating leukemia stem cells are affected by IMATINIB, while quiescent leukemia stem cells are not affected. Due to this additional assumption, the leukemia population under IMATINIB does not relapse without the effects of IMATINIB resistance mutations. Instead, under IMATINIB treatment, the leukemia stem cell population restabilizes at lower equilibrium level and does not continue growing as in the Michor model.
Both [5] and [6] propose that IMATINIB does not eliminate the leukemia stem cell population. Consequently, the papers conclude that IMATINIB therapy should be combined with an additional treatment that either directly impacts leukemia stem cells or causes leukemia stem cells to become vulnerable to IMATINIB.
As an alternative approach, Komorova and Wodarz develop a model that focuses on the drug resistance of leukemia cells [3]. In their model, they implicitly assume that IMATINIB affects all leukemia cells including stem cells and that inevitable relapse is a result of acquired IMATINIB resistance mutations. Komorova and Wodarz consider the possibility of treating patients with multiple drugs to reduce the probability of any leukemic cell eventually acquiring resistance-mutations to all drugs. They determine that a treatment strategy consisting of three leukemia-targeted drugs of different specificity might have a strong chance of eliminating the disease.
The four approaches discussed above present a variety of hypotheses for the dynamics of IMATINIB treatment on leukemic cells. These papers also propose potential treatment strategies to enhance the effectiveness of IMATINIB. However, the difficulty with these treatments is that it is unclear what kind of drug could be used to target leukemia stem cells or what alternative drugs could be used in addition to IMATINIB for a multiple-drug strategy.
In this work, we consider the following more general mathematical model which is an extension of a model proposed in [5]. In our model, we assume that normal (resp. leukemic) cells differentiate through two stages of their life cycle, beginning with leukemic stem cells which produce produce progenitors.
The mathematical form of the system we shall investigate satisfieṡ with the following initial conditions The population of interest is divided into three compartments coming from dictated by the epidemiological stages; normal cells, sensitive leukemic cells and resistant leukemic cells. We assume that normal (resp. leukemic) cells differentiate through two stages of their life cycle, beginning with leukemic stem cells which produce. Tables 1 and 2 list the definitions and symbols (populations and parameters) used in our model. This paper is organized as follows. In the next section, we state and prove general criterion for the existence positive solutions of system (1). In Sect. 3, we will study the global dynamics of (1) by constructing a suitable Lyapunov function and using LaSalle's invariance principle rather than by using the theory of competitive systems, as has been done in [6]. This will 01005-p.2 production rate of the normal stem cells a y produce rate of the leukemic stem cells d 1 death rates of the normal progenitors cells d 2 death rates of the leukemic progenitors cells d 3 (r) death rates of the normal leukemic progenitors cells r resistant parameter 0 < < 1 enable us to obtain the global asymptotic stability of the equilibrium point under some hypotheses and by a simpler method. This will enable us to obtain the global asymptotic stability of the equilibrium point under some hypotheses and by a simpler method. In Sect. 4, we use numerical simulations of our model to discuss biological significance of our results and indicate possible extensions to the study of more comprehensive models. The model (1) becomes as follow

WMLS 2014
with the following initial conditions Throughout this paper, we assume that conditions 01005-p.3 and hold in order that the state variables x 0 , x 1 , y 0 , y 1 , z 0 , z 1 have a biological meanings. The compartmental diagram of the system (3) appears in Fig. 1.

Global Wel-Posedness
First we show that solutions of system (3) are positive and bounded.

WMLS 2014
Observe that F is Lipschitz continuous on bounded sets of R 6 . Furthermore, for t ≥ 0 and Hence we obtain the boundedness and t ≥ 0, this implies boundedness of solutions.
Thus the existence of a unique global nonnegative bounded solution is proved.

Case r = 0
In this section, we consider the case where the resistant parameter r = 0 with z 0 = z 1 = 0. We obtain the following system which we call the endemic equilibrium. The expression of i and i are given in Table 3.

Local asymptotic stability
To analyze the stability of these steady states we compute their linearizations.
At the point , we obtain the linearization 2. At point E 1 we get 1 = − 0 1 and 2 = − a y − a x d 1 1 . The eigenvalues of A are 2 , −d 2 and the roots of the characteristic polynomial of the minor matrix of A given by This equation is equivalent to P (u) = u 2 + e 1 u + e 2 = 0, where e 1 = d 1 − 1 and e 2 = −d 1 1 + a x 1 . Since e 1 > 0 and e 2 > 0, the Ruth-Hurwitz criterion implies that all roots of P have negative real parts. Hence eigenvalues of A have negative real parts if 2 < 0, which shows that the equilibrium E 1 is locally asymptotically stable if q < 1 and d 1 < d * 1 . 3. At point E 2 we have 1 = − a x − a y d 2 2 and 2 = − 0 2 . The eigenvalues of A are 1 , −d 1 and the roots of the characteristic polynomial of the minor matrix of A given by This equation is equivalent to where e 1 = d 2 − 2 and e 2 = −d 2 2 + a y 2 . Since e 1 > 0 and e 2 > 0, the Ruth-Hurwitz criterion implies that all roots of P have negative real parts. Hence eigenvalues of A have negative real parts if 1 < 0, which shows that the equilibrium E 2 is locally asymptotically stable if q > , d 2 < d * 2 . 4. At point E 3 we have 1 = − 0 3 and 2 = − 0 3 . The eigenvalues of A are the roots of the characteristic polynomial of the minor matrix of A given by where In summary we have Table 4. Summary of the model with r = 0.
Case 1: q ≤ Case 1: < q < 1 Case 3: 1 ≤ q region stability region stability region stability Bifurcation diagram for the model when < q < 1. In region I, the disease free equilibrium E 1 is locally asymptotically stable, in region II, the healthy free equilibrium E 2 is locally asymptotically stable, in region III, the endemic equilibrium E 3 is locally asymptotically stable and region IV represents the local asymptotic stability of disease free equilibrium E 1 and healthy free equilibrium E 2 .

Concept of R 0
The disease free equilibrium (DFE) for this nondimensionalized general model of chronic myeloid leukemia may be used to find the basic reproduction number R 0 , which indicates the average number of new infections. The basic epidemiological reproductive number is given by However, this nondimensional number is not enough to characterize the dynamics of model (4), (10).

Analysis at
To analyze the local stability in critical cases, we use the center manifold theory, as described in Castillo-Chavez and Song (Theorem 4.1) [1]. Center manifold theory has been used to decide the local stability of a nonhyperbolic equilibrium (linearization matrix has at least one eigenvalue with zero real part). This theory can not only determine the local stability of the nonhyperbolic equilibrium but also settles the question of the existence of another equilibrium (bifurcated from the nonhyperbolic equilibrium).
To apply this method, the system (10) becomeṡ At R 0 = 1, we obtain = − a y = ( −a x ) a x 0 d 1 + a x . From the linearization matrix of system (13) around the disease free equilibrium E 1 when R 0 = 1, we can show that zero is a simple eigenvalue and all other eigenvalues have negative real parts. A right eigenvector w corresponding to the zero eigenvalue is , and the left eigenvector satisfying v.w = 1 is v = 0 0 , 1 0 .
For system (13) the second partial derivatives of f 3 are given by  Bifurcation diagram for the model when q < q 1 (Resp. q > q 2 ). In region I, the disease free equilibrium E 1 is globally asymptotically stable, in region III, the endemic equilibrium E 3 is globally asymptotically stable (Resp. In region II, the healthy free equilibrium E 2 is globally asymptotically stable, in region III, the endemic equilibrium E 3 is globally asymptotically stable).
It follows that 1. If d 2 < d 2 < d * 2 , the unique endemic equilibrium disappears whenever R 0 > 1 and is close to 1. 2. If d 2 > d * 2 > d 2 , the unique endemic equilibrium is locally asymptotically stable whenever R 0 > 1 and is close to 1. 3. If d 2 < d 2 and d 1 < d 1 , the unique endemic equilibrium disappears whenever R 0 > 1 and is close to 1.

Global stability
Next, we study the global stability of the steady states. The proof relies on the construction of a global Lyapunov functions.
Theorem 3.5: , then E 2 is globally asymptotically stable in R 4 Proof: 1. Let us consider the Lyapunov function It is easily seen that V 1 ≥ 0 and V 1 = 0 if and only if x 0 = 1 , x 1 = a x d 1 1 and y 0 = y 1 = 0. Calculating the time derivative of V 1 along the positive solutions of model (10), we obtaiṅ . It follows from q < 1 and d 1 ≤ d * 1 thatV 1 ≤ 0 for all x 0 , x 1 , y 0 , y 1 andV 1 = 0, when x 0 = 1 , x 1 = a x d 1 1 and y 0 = y 1 = 0. Hence La Salle's theorem [4] implies convergence of the solutions to this equilibrium, for all initial values not in the set {0} × R 3 + . This shows that the disease free equilibrium E 1 is globally asymptotically stable in R 4 + / {0} × R 3 + . If the initial data starts from {0} × R 3 + , then the solution obviously converges to the equilibria E 0 or E 2 . 2. Let us consider the Lyapunov function It is easily seen that V 2 ≥ 0 and V 2 = 0 if and only if x 0 = x 1 = 0, y 0 = 2 and y 1 = a y d 2 2 . Calculating the time derivative of V 2 along the positive solutions of model (10), we obtaiṅ . It follows from q > and d 2 ≤ d * 2 thatV 2 ≤ 0 for all x 0 , x 1 , y 0 , y 1 andV 2 = 0, when x 0 = x 1 = 0, y 0 = 2 and y 1 = a y d 2 2 , hence La Salle's theorem [4] implies convergence of the solutions to this equilibrium, for all initial values not in the set R 2 + × {0} × R + . This shows that the disease equilibrium E 2 is globally asymptotically stable in R 4 + /R 2 + × {0} × R + . If the initial data starts from R 2 + × {0} × R + , then the solution obviously converges to the equilibria E 0 or E 1 . and d 2 > max(d * 2 , d • 2 ) thatV 3 ≤ 0 for all x 0 , x 1 , y 0 , y 1 andV 3 = 0 for x 0 = 3 , x 1 = a x d 1 3 , y 0 = 3 and y 1 = a y d 2 3 . Hence La Salle's theorem [4] implies convergence of the solutions to this equilibrium, for all initial values not in the set {0} × R 3 + R 2 + × {0} × R + . This shows that the endemic equilibrium E 3 is globally asymptotically stable in R 4 + / {0} × R 3 + R 2 + × {0} × R + . If the initial data starts from {0} × R 3 + R 2 + × {0} × R + , then the solution obviously converges to the equilibria E 0 , E 1 or E 2 .