MATHEMATICAL ANALYSIS OF A MODEL FOR CHRONIC MYELOID LEUKEMIA

In this paper, a mathematical analysis of a model describing the evolution of chronic myeloid leukemic with effect of growth factors is considered. The corresponding dynamics are represented by a system of ordinary differential equations of dimension 5. This system described the interactions between hematopoietic stem cells (H.S.C), hematopoietic mature cells (M.C), cancer hematopoietic stem cells, cancer hematopoietic mature cells and the associated growth factor concentration. Our research is, henceforth, carried out on the existence and the uniqueness of the solution of this system. The next substantive concern will be a discussion on the local and global stability of the corresponding steady states. Three scenarios, however, corresponding to different actions of hematopoiesis on stem cells (differentiate cells or both cells) are considered.


Introduction
Chronic myeloid leukemia (CML) is a cancer of the bone marrow and blood. It accounts about 15 percent of newly diagnosed cases in leukemia in the world. Most of these are adults with an average of diagnosis in 64 years, but rare are those cases in children. CML grows very slowly over years in the sense that a patient may have it for long time before symptoms are noticed. Common signs of CML are anemia, splenomegaly, tiredness, weight loss, discomfort and goes through 3 phases: chronic phase, accelerated phase and blast phase. It affects a specific type of white blood cells which are known as the myelocytes. In fact, it is myelo-proliferative disorder originating from the myeloid hematopoietic stem cell. In turn, this results from the clonal expansion of pluripotent hematopoietic stem cells containing the active BCR-ABL fusion gene produced by a reciprocal translocation of the ABL gene to the BCR gene in chromosomes 22 and 9 (see [13,17]). This new chromosome is named Philadelphia chromosome [4]. Through mitosis or division, thanks of growth factor, the multiplication of healthy cells is abundant and cancerous cells do not respect the cellular mechanisms in their proliferation. After being stimulated by a physiological signal, stem cells at rest start their renewing and differentiating. The balance between them is named Homeostasis. Thus CML has become one of the most extensively studied human malignancy. Even if many interesting studies on CML are given (see [10,19,22]), many information are unknown about the dynamic of how cancer cells are propagated. During the last decade, many several advances are mainly based, particulary, on advances in scientific solutions to capture the dynamics of CML, one of those promising dynamics approaches includes mathematical modeling by identifying interactions between all types of cells playing rule in propagation of leukemia and using estimate parameters based on experimental advances. CML mathematical models do not suggest an exact solution but provide useful results. Therefore, it is still a great need to continue studying and developing existing models in order to adapt biological discovers and recent clinical results and associate therapies [15,25]. Those models are often represented by ordinary differential equations, partial differential equations or delay differential equations that represent different states of stem cells. For more details see for example [5-7, 11, 20, 23, 27, 31, 32].
The models given in [1,2,14,26] describe particularly the dynamics of H.S.C population. Dingli and Michor [15] assumed that the H.S.C cancerous cells compete with normal H.S.C cells. In their model, the regeneration of H.S.C is governed by homeostasis, which controls the process of dividing according to the total number of H.S.C. They have developed model that involve H.S.C cells and M.C mature cells, where proliferation and death rates are introduced and interactions are also considered between x 0 , x 1 , y 0 and y 1 .
The regeneration of (H.S.C) is generated by Homeostasis of x 0 represented by a function ϕ whereas that of y 0 is represented by function ψ (see [27]).
They have considered the following model: The parameters used in (1.1) -(1.2) are given in Table 1 given bellow and come from [9]. The value of K (the bone morrow receiving capacity) changes in each scenario as we shall see after.  In this paper, we propose an extension of this model taking into account the influence of the growth factor E on leukemia disease. This factor is studied particularly in [30] and plays an important rule in evolution of CML.
Our paper is organized as follows: In Section 2, the mathematical model of leukemia is proposed in details and three possible scenarios are set: • scenario 1 corresponds to ε 1 = 1 and ε 2 = 0 where homeostasis acts only on stem cells x 0 and y 0 . • scenario 2 corresponds to ε 1 = 0 and ε 2 = 1 where homeostasis acts only on differentiate cells x 1 and y 1 . • scenario 3 corresponds to ε 1 = 1 and ε 2 = 1 where homeostasis acts on stem cells x 0 and y 0 and on differentiate cells x 1 and y 1 .
In Section 3, the existence and uniqueness of a positive local solution of our model is proved for all scenarios, existence and uniqueness of a positive global solution is proved for scenario 1 and 3. The existence of steady states for the three scenarios is studied in Section 4 (trivial, blast, no pathologic and chronic states). In Section 5, local stability analysis for those steady states is developed followed in Section 6 by a global stability analysis. Biological interpretation of the obtained results is proposed and some perspectives are set in Section 7.

Mathematical Model
In this paper, we consider the following five-dimensional model which is an extension of the one proposed in [3,9]: where safe (M.C) cells x 1 proliferate at a rate d 2 and are eliminated at a rate In fact, the growth factor concentration E follows the evolution equation [32] .
The function f acts as negative feedback of the non proliferating (H.S.C) population on the production of growth factor.
We assume that f is positive and decreasing and according to [2]: The associated growth factor concentration is given by The studied population is divided into two compartments: hematopoietic stem cells (H.S.C) and differentiated cells (M.C).
In model (2.3), the death rate of normal differentiated cells d(E) undergoes a Hill function evolution and according to [2], d(E) is given by In addition to the Table 1, the Table 2 lists the parameters used in our model. All those parameters are assumed to be positive.

Parameter Explanation d(E)
Death rate of norma(l M.C) x 1 K 0 Clearing rate of growth factors K 1 Rate of maximum saturation of growth factors K 2 Rate of half saturation of growth factors r 0 Oscillation rate a Absorption rate of E by cells 3. Basic properties of model (2.3)

Existence of a positively invariant attracting set.
Proposition 3.1. The system (2.3) is positively invariant in the following cone: Proof. One has at (0, x 1 , y 0 , y 1 , E) : at (x 0 , 0, y 0 , y 1 , E) : Then the vector field is pointing in the direction of D and does not leave it, so according to [28], this system is positively invariant.

Existence and uniqueness of solution.
The Cauchy problem associated to the model (2.3) is given by

. A local solution of the Cauchy problem associated to (2.3) exists and is unique in D.
Proof. Since the function F is of class C 1 , then a local solution of the Cauchy problem associated to (2.3) exists and is unique in D. This is due to Picart-Lindeloff theorem [24]. Now, let us prove the global existence of the solution of the Cauchy problem associated to (2.3) for scenarios 1 and 3. In fact, it is sufficient to prove that the corresponding solution is bounded in a convenient set Γ of D. Here we shall consider

Proposition 3.3. The Cauchy problem associated to model (2.3) admits a global and unique solution defined on Γ for scenarios 1 and 3.
Proof. The first and third equations of model (2.3) are given in scenario 1 by and in scenario 3 by For those two scenarios, one has those corresponding majorations The solution of (3.2) and (3.4) can be compared to the solution of the following Bernoulli equation with the same initial condition According to the comparison Theorem [24], the solutions of (3.2) and (3.4) satisfy for all t ≥ 0, The solutions of (3.3) and (3.5) can be compared to the solution of the Bernoulli equation with the same initial condition and y 0 (0) is assumed to be different from 0.
According to the comparison Theorem, the solutions of (3.3) and (3.5) satisfy for all t ≥ 0, .
The second equation of (2.3) for scenarios 1 and 3 is given by In this case, the solution of (3.6) can be compared to the solution of the following differential equation with the same initial condition which is According to the comparison theorem, the solution of (3.6) satisfies for all t ≥ 0: The forth equation of (2.3) for scenario 1 and 3 is given by In this case, the solution of (3.7) can be compared to that of the following differential equation with the same initial condition which is According to the comparison theorem, the solution of (3.7) satisfies for all t ≥ 0 .
The fifth equation of (2.3), for scenarios 1 and 3, is given by In this case, the solution of (3.8) can be compared to the solution of the following differential equation with the same initial condition which is According also to the comparison Theorem, the solution of (3.8) satisfies for all t ≥ 0 Hence, Finally, for scenarios 1 and 3, any solution of (2.3) that starts in R 5 + is confined in Γ and since Γ is compact and positively invariant for model (2.3), according to [28], there exists an unique global solution of the Cauchy problem associated to (2.3) in Γ (for scenarios 1 and 3).

Existence of steady states
The steady states of (2.3) are the following: • The trivial steady state S 0 = (0, 0, 0, 0, E 0 ) corresponds to the extinction of cell population, where E 0 = a K 0 > 0.
• The no pathologic steady state S np = (x 0,np , x 1,np , 0, 0, E np ) corresponds to presence of normal cells without leukemic cells, where • The blast steady state S b = (0, 0, y 0,b , y 1,b , E b ) corresponds to presence of leukemic cells without normal cells, where • The chronic steady state S c = (x 0,c , x 1,c , y 0,c , y 1,c , E c ) corresponds to simultaneous presence of normal cells and leukemic cells, where Theorem 4.1. For all the three scenarios: • The trivial steady state S 0 always exists.
• If n > d 0 then the no pathologic steady state S np exists.
• If m > g 0 then the blast steady state S b exists.
Then the non pathologic steady state S np exists if n > d 0 .
Then the chronic steady state S c exists if T 1 < 1 < T 2 .
Then the chronic steady state S c exists if T 1 < 1 < T 2 .
Then the no pathologic steady state S np exists if n > d 0 .
Then the chronic steady state S c exists if T 1 < 1 < T 2 .

Local stability analysis
The Jacobian matrix J of system (2.3) is given by Proof. In this case, J is rewritten (and denoted J 1 ) for this scenario as

Denote by
The corresponding eigenvalues of J 1 (X) are: Hence J 1 (X) has three negatives eigenvalues λ 1 , λ 2 and λ 3 for all steady states.
• At the trivial steady state. One has Then the trivial steady state is LAS if n < d 0 and m < g 0 .
• At the no pathologic steady state. One has Then, the no pathologic steady state is LAS if T 1 < 1.
• At the blast steady state. One has Then the blast steady state is LAS if T 2 > 1.
• At the chronic steady state. One has Then λ 4 and λ 5 satisfy Then the chronic steady state is unstable.

The trivial steady state is LAS if n < d 0 and m < g 0 . 2. The no pathologic steady state and the blast steady state are LAS if
The blast steady state is the unique LAS state if T 2 > 1 and T 1 > 1.

The no pathologic steady state is the unique LAS state if T 1 < 1 and
T 2 < 1.

The chronic steady state is unstable.
Proof. The Jacobian matrix is rewritten (and denoted J 2 ) for this scenario as
• At chronic steady state. One has The corresponding characteristic polynomial is given by P (λ) = a 5 λ 5 + a 4 λ 4 + a 3 λ 3 + a 2 λ 2 + a 1 λ + a 0 where a 5 = 1, The associated Hurwitz matrix is given by a 1 a 4 ), So from Hurwitz criterion [21], as at least one element of the first column of M (here e 0 ) is negative then the chronic equilibrium state is unstable.

Scenario 3.
Proposition 5.3. Proof. The Jacobian matrix is rewritten (and denoted J 3 ) for this scenario as

Denote by
Then • At the trivial steady state. One has J 3 (X) has three negatives eigenvalues: So the trivial steady state is LAS if n − d 0 < 0 and m − g 0 < 0. • At the no pathologic steady state.
The corresponding characteristic polynomial is given by The associated Hurwitz matrix M ′ is given by . Then from the Hurwitz criterion since one of the first column of M ′ is negative (here c ′ 0 ), the no pathologic steady state is unstable. • At the blast steady state. One has 3 . J 3 (X) has four negatives eigenvalues: Then the blast steady state is L.A.S if T 2 > 1.
• At the chronic steady state. One has The corresponding characteristic polynomial is given by The associated Hurwitz matrix M ′′ is given by So at least one element of the first column of M ′′ is negative(here e ′′ 0 ) then the chronic steady state is unstable(from Hurwitz criterion ).

Global stability analysis
In this section, global stability analysis for steady states of model (2.3) is proposed for scenarios 1 and 3 .

Study of global analysis of no pathologic and blast steady states for Scenario 1.
Recall that according to the theorem 4.1 given in section 4, the no pathologic steady state exists for n > d 0 and the blast steady state exists for m > g 0 .So, this is assumed in this section to deal with those two states. To analyze the global stability of those steady states of system (2.3), we shall use the following theorem given in [28].
In order to apply this theorem, let us split our system (2.3) into two subsystems: the first one is a subsystem in (x 0 , y 0 ) and the second one in (x 1 , y 1 , E) So, first consider the following subsystem of (2.3) in (x 0 , y 0 ) under initial conditions x 0 (0), y 0 (0) From Proposition 3.3, one has that the solution of (6.2) satisfies x 0 (t) ≤ m 1 and y 0 (t) ≤ m 2 for all t ≥ 0. According to this result, let us define the positively invariant compact set as Lemma 6.2. The system (6.2) has no limit cycle in intB, where intB is the interior of B.
Proof. As the subsystem (6.2) has no limit cycle in the bounded set intB ⊂ R 2 + , the GAS results are obtained from a direct application of the Poincaré-Bendixon theorem [8] to this subsystem. Now let us consider the second subsystem Using, the obtained results on L.A.S of no pathologic and blast steady states in section 5 for system (2.3), we underline that those results are available also for its subsystems (6.3).
Moreover, one can directly see that by direct integration of subsystem (6.3) and majorations that (x 0,np , y 0,np , x 1,np , y 1,np Hence according to Theorem 6.1, the non pathologic steady state is GAS if T 1 < 1 < T 2 and the blast steady state is GAS if T 2 > 1 and T 2 > 1 for model (2.3 The components of the blast steady state are given by . Denote by y 0,b + y 1,b = K α 1 − g0 m and g 1 = q y 0,b y 1,b . Then the system (6.4) is rewritten as To prove the global stability of blast steady state, let us construct an appropriate Lyapunov function and consider the following function V defined by V (x 0 , x 1 , y 0 , y 1 where α i , i = 1, ..., 5 are positive constants (that we will choose later).
Knowing that thus for all and also Moreover, and so after replacing, one has The coefficients α i where i = 1, 2, 3, 4, 5 will be chosen such that α 2 r = α1n K and α4q In this case, T and consider the following matrix: So finally, let us choose α 1 = α 3 = α 5 = 2, nα 1 = mα 3 , α 2 = 2n rK and α 4 = . Thus the matrix Π becomes One obtains −Π = R(−W ) + (−W T )R. Since R is a diagonal matrix and (−W ) has all its eigenvalues negatives, thus (−Π) is positive definite matrix and then Π is a negative definite matrix, this implies that dV (x 0 , x 1 , y 0 , y 1 , E) dt < 0.
Then, according to [29], the function V is a Lyapunov function for the system (6.4) and the blast steady state is G.A.S for scenario 3.

Interpretation results and some perspectives
The main results of this study are resumed in this two tables: Steady states Existence conditions Stability conditions trivial steady state always exists n < d 0 and m < g 0 no pathologic steady state n > d 0 T 2 > 1 blast steady state m > g 0 T 1 < 1 chronic steady state T 1 < 1 < T 2 unstable Table 3. Local and global stability in scenario 1 and local stability in scenario 2 In this study, we have shown that for any positive initial conditions the system (2.3) admits a unique global positive solution for scenarios 1 and 3. The trivial steady state is LAS if n and m (the proliferation rates) are lower than the mortality rates d 0 and g 0 respectively, however these conditions are clinically impracticable then the trivial balance is unstable. Furthermore, the chronic steady state being a saddle point, the coexistence of normal and

Steady states
Existence conditions Stability conditions trivial steady state always exists n < d 0 and m < g 0 no pathologic steady state n > d 0 unstable blast steady state m > g 0 T 1 < 1 chronic steady state T 1 < 1 < T 2 unstable Table 4. Local and global stability in scenario 3 cancer cells does not maintain for long time and this is realist with biological observations. The blast steady state in the third scenario is also a saddle point. Otherwise the local and global stability of the blast and no pathologic steady states is linked to coefficients T 1 and T 2 : If T 1 > 1 and T 2 > 1, the no pathologic steady state is the only stable state in this case there is no leukemia where all types of cells are in no pathological steady state therefore there is no cancer cells or their proliferation is practically zero.
If T 1 < 1 < T 2 , the blast and no pathologic steady states are both LAS. In this case the state converges towards no pathologic steady where the number of healthy cells is higher than that corresponding to cancer cells. There may even be a conversion of cancerous cells into healthy ones or else the dynamics of our model converge towards blast steady state and cancer cells are in majority (if not all cells are diseased).
If T 1 > 1 and T 2 < 1, all steady states are unstable, it is corresponds to the final phase of leukemia.
Note that the overall stability of the steady states has been demonstrated in scenario 1 and scenario 3. The study of global stability in the scenario 2 is not possible by considering mathematical classic tools. We let it in a future work using simulations.
Biological interpretation of the results must be also developed. In parallel and in the case where the blast and no pathological steady states are LAS at the same time, a control could possibly be introduced on the growth factors so that the system converges towards the no pathological steady state, this will be also investigated in a future work.