*Article* **Stability Analysis of Equilibria for a Model of Maintenance Therapy in Acute Lymphoblastic Leukemia**

**Irina Badralexi 1,\*,†, Andrei-Dan Halanay 2,† and Ragheb Mghames 3,†**


**Abstract:** In this paper, we study two mathematical models, involving delay differential equations, which describe the processes of erythropoiesis and leukopoiesis in the case of maintenance therapy for acute lymphoblastic leukemia. All types of possible equilibrium points were determined, and their stability was analyzed. For some of the equilibrium points, conditions for parameters that imply stability were obtained. When this was not feasible, due to the complexity of the characteristic equation, we discuss the stability through numerical simulations. An important part of the stability study for each model is the examination of the critical case of a zero root of the characteristic equation. The mathematical results are accompanied by biological interpretations.

**Keywords:** delay differential equations; critical case for stability; acute lymphoblastic leukemia; maintenance therapy

#### **1. Introduction**

#### *1.1. Mathematical Background*

The primary goal of this paper is to present the stability analysis of two different mathematical models for the processes of erythropoiesis and leukopoiesis in the case of maintenance therapy in acute lymphoblastic leukemia (for more details, see [1–6]).

The models represent the original work of some of the authors and were first introduced in [7]. They consist of systems of delay differential equations (often used to capture cell dynamics). In [7], the authors presented the models and demonstrated the positivity of solutions and the existence of equilibria. The novelty of this paper is the thorough stability study of all of the equilibrium points.

A challenging situation was the presence of critical cases for stability analysis. For these, we used an original theorem proven in [8] by some of the authors of this paper.

#### *1.2. Biological Background*

The models presented in this paper depict the erythropoiesis and leukopoiesis processes. These are part of a bigger process called hematopoiesis: the process through which all blood cells are created. Hematopoietic stem cells become red blood cells (which transport oxygen), white blood cells (which fight infections), or platelets (which stop bleeding). There is a complex network of cytokines and growth factors that regulate the production of blood cells.

Hematopoietic stem cells generate two major progenitors cell lineages: myeloid and lymphoid. The myeloid line contains cells such as granulocytes, monocytes, erythrocytes,

**Citation:** Badralexi, I.; Halanay, A.-D.; Mghames, R. Stability Analysis of Equilibria for a Model of Maintenance Therapy in Acute Lymphoblastic Leukemia. *Mathematics* **2022**, *10*, 313. https:// doi.org/10.3390/math10030313

Academic Editors: Ioannis G. Stratis and Dumitru Baleanu

Received: 4 November 2021 Accepted: 12 January 2022 Published: 20 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

or platelets, while the lymphoid line contains cells associated with the immune system (i.e., natural killer cells, T-cells, and B-cells).

Erythropoiesisis the process through which red blood cells (also called erythrocytes) are produced. Erythropoietic stem cells mature to red blood cells when a decrease in oxygen levels is detected by the kidneys. The kidneys secrete a hormone called erythropoietin, which stimulates the production of erythrocytes. It is the action of this hormone (also considered in [2,9]), together with the effects of therapy, that complicates the stability study of equilibria. Since mortality rates are no longer constant, a new variable is introduced. This variable is the source of a zero root of the characteristic equation.

White blood cells are generated through the process of leukopoiesis. Depending on which progenitor line they come from, there are two groups of white blood cells: myelocytes and lymphocytes.

For more information about hematopoiesis and its underlying processes, please see [2,3,10,11], pp. 1–18.

Acute lymphoblastic leukemia (ALL) is a type of cancer that affects white blood cells [4]. In the case of ALL, blood cell production is disrupted, due to the presence of a large number of immature lymphocytes (called lymphoblasts) in the circulatory system. These immature cells do not function properly and overcrowd the healthy cells. Lymphoblasts eventually invade the liver, spleen, and lymph nodes. For more information, statistics, diagnoses, and prognoses, see [12], pp. 1556–1576, 1616–1636, and [10], pp. 173–192.

Maintenance therapy is conducted after an initial treatment is administered [6]. In patients with ALL, we consider that this therapy consists of oral administration of mercaptopurine (6-MP). The biologically inactive substance 6-MP metabolizes into the active 6-thioguanine nucleotide (6-TGN).

Studies have shown that maintenance therapy contributes greatly to the overall evolution and life expectancy of the patient ([12], p. 1566, [13]). The need to mathematically study the effects of this treatment is very important, as it may help to fill the existing gaps in knowledge [13].

The models studied in this paper capture the effects of this medication on the creation of red and white blood cells.

It is worth mentioning that, although we concentrated on patients with ALL, the same maintenance therapy is administrated for some autoimmune diseases. Thus, the models, computations, and discussions could also be suitable for those cases.

#### **2. The Modeling of Erythropoiesis**

#### *2.1. The Mathematical Model*

The mathematical model for erythropoiesis consists of seven delay differential equations (DDEs) with two delays. We considered the dynamics of the stem-like short-term erythroid cells (*u*1), the mature erythrocytes (*u*2), the concentration of erythropoietin (*u*3), the amount of 6-MP in the gut (*u*5), the amount of 6-MP in plasma (*u*5), and the concentration of thioguanine nucleotide (6-TGN) in red blood cells (*u*7) (see [7,9,14,15]). The variable *u*<sup>4</sup> represents the loss suffered during the cell cycle.

Each stem-like cell is considered to go through either asymmetric division (a fraction *η*1*<sup>e</sup>* of the stem-like population), which means that one of the daughter cells re-enters the stem-like population, while the other goes on to differentiate, symmetric differentiation (a fraction *η*2*<sup>e</sup>* of the stem-like population)—both daughter cells proceed to differentiate—or self-renewal (a fraction 1 − *η*1*<sup>e</sup>* − *η*2*<sup>e</sup>* of the stem-like population)—both daughter cells re-enter the stem-like pool. We used two feedback functions for the rates of self-renewal and differentiation:

$$
\beta\_{\mathfrak{c}}(\mathfrak{u}\_1, \mathfrak{u}\_3) = \beta\_{0\mathfrak{c}} \frac{1}{1 + \mathfrak{u}\_1^{\mathfrak{m}}} \cdot \frac{\mathfrak{u}\_3}{1 + \mathfrak{u}\_3}
$$

$$
k\_{\mathfrak{c}}(\mathfrak{u}\_3) = k\_{0\mathfrak{c}} \frac{\mathfrak{u}\_3}{1 + \mathfrak{u}\_3}
$$

The time necessary for a stem-like cell to go through self-renewal or differentiation and re-enter the stem-like population is considered to be the same and is represented by the delay *τ*1.

These aspects are included in the first equation of the model, which describes the evolution of stem-like short-term erythroid cells. As can be seen in this equation, both the erythropoietin and the medication have an impact on cell death. The drug concentration in the blood stream has a toxicity that increases the apoptosis rate of the stem-like cells (the second term). The erythropoietin regulates the process of erythropoiesis, either repressing (the first term) or stimulating (by influencing the rate of self-renewal and differentiation) the production of stem-like cells.

The constant *R*˜ *<sup>m</sup>* represents the maximum drug effect on erythrocytes, and *R*˜ <sup>50</sup> represents the drug saturation.

The healthy erythrocytes, illustrated by the second equation, have a constant mortality rate (the first term) and a supply controlled by the erythropoietin (the second term). We introduce the notation *A*˜*<sup>e</sup>* = *Ae*(2*η*1*<sup>e</sup>* + *η*2*e*), where *Ae* is an amplification factor. The time necessary for the cells to mature is given by the delay *τ*2.

The third equation captures the changes in the concentration of erythropoietin. This protein has a constant rate of clearance, and its production is influenced by the number of existing erythrocytes in the blood stream.

The fourth equation represents the cell loss sustained during the cell cycle. For more information regarding this equation, please consult [7].

The last three equations describe the dynamics of the drug after administration.

The model, which depicts the process of erythropoiesis under maintenance therapy, is:

$$\dot{u} = f\_i(u, u\_{r\_j}),\\ \dot{\iota} = \overline{1, 7},\\ \dot{\jmath} = \overline{1, 2} \tag{1}$$

*<sup>u</sup>*˙1 <sup>=</sup> <sup>−</sup> *<sup>γ</sup>*<sup>0</sup> 1 + *u<sup>α</sup>* 3 *<sup>u</sup>*<sup>1</sup> <sup>−</sup> *<sup>R</sup>*˜ *mu*<sup>7</sup> *R*˜ <sup>50</sup> + *u*<sup>7</sup> *u*<sup>1</sup> − (*η*1*<sup>e</sup>* + *η*2*e*)*ke*(*u*3)*u*<sup>1</sup> − (1 − *η*1*<sup>e</sup>* − *η*2*e*)*βe*(*u*1, *u*3)*u*<sup>1</sup> +2*u*4(1 − *η*1*<sup>e</sup>* − *η*2*e*)*βe*(*u*1*τ*<sup>1</sup> , *u*3*τ*<sup>1</sup> )*u*1*τ*<sup>1</sup> + *η*1*eu*4*ke*(*u*3*τ*<sup>2</sup> )*u*1*τ*<sup>1</sup> *<sup>u</sup>*˙2 <sup>=</sup> <sup>−</sup>*γ*2*u*<sup>2</sup> <sup>+</sup> *<sup>A</sup>*˜*eke*(*u*3*τ*<sup>2</sup> )*u*1*τ*<sup>2</sup> *<sup>u</sup>*˙3 <sup>=</sup> <sup>−</sup>*ku*<sup>3</sup> <sup>+</sup> *<sup>a</sup>*<sup>1</sup> 1 + *u<sup>n</sup>* 2 *u*˙4 = *u*<sup>4</sup> <sup>−</sup> *<sup>γ</sup>*<sup>0</sup> 1 + *u<sup>α</sup>* 3 <sup>−</sup> *<sup>R</sup>*˜ *mu*<sup>7</sup> *R*˜ <sup>50</sup> + *u*<sup>7</sup> <sup>+</sup> *<sup>γ</sup>*<sup>0</sup> 1 + *u<sup>α</sup>* 3*τ*<sup>1</sup> + *<sup>R</sup>*˜ *mu*7*τ*<sup>1</sup> *<sup>R</sup>*˜ <sup>50</sup> + *<sup>u</sup>*7*τ*<sup>1</sup> *u*˙5 = −*b*1*u*<sup>5</sup> + *a*<sup>2</sup> *<sup>u</sup>*˙6 <sup>=</sup> *<sup>b</sup>*1*u*<sup>5</sup> <sup>−</sup> *<sup>e</sup>*1*u*<sup>6</sup> <sup>−</sup> *<sup>c</sup>*1(<sup>1</sup> <sup>−</sup> *<sup>e</sup>*2) *c*<sup>2</sup> + *u*<sup>6</sup> *<sup>u</sup>*<sup>6</sup> <sup>−</sup> *<sup>m</sup>*2*e*<sup>2</sup> *m*<sup>1</sup> + *u*<sup>6</sup> *u*6 *<sup>u</sup>*˙7 <sup>=</sup> *<sup>b</sup>*2*c*1(<sup>1</sup> <sup>−</sup> *<sup>e</sup>*2) *c*<sup>2</sup> + *u*<sup>6</sup> *u*<sup>6</sup> − *e*3*u*7.

#### *2.2. The Equilibrium Points and Linearization*

Using some elementary, but tedious calculations, the following types of equilibrium points corresponding to the model of erythropoiesis in ALL under treatment were obtained (see also [7]): *E*<sup>1</sup> = (0, 0, *u*ˆ3, *u*ˆ4, *u*ˆ5, *u*ˆ6, *u*ˆ7), which corresponds to the "death of the patient", and *E*<sup>2</sup> = (*u*ˆ1, *u*ˆ2, *u*ˆ3, *u*ˆ4, *u*ˆ5, *u*ˆ6, *u*ˆ7), which corresponds to a "chronic phase of the disease".

In order to study the stability of these equilibrium points, we performed a linearization of the nonlinear system (1). We denote the matrix of partial derivatives for the undelayed variables by *<sup>A</sup>* <sup>=</sup> *<sup>∂</sup> <sup>f</sup> <sup>∂</sup><sup>u</sup>* = [*aij*] and the matrices of the partial derivatives with respect to the delayed variables by: *<sup>B</sup>* <sup>=</sup> *<sup>∂</sup> <sup>f</sup> ∂uτ*<sup>1</sup> = [*bij*] and *<sup>C</sup>* <sup>=</sup> *<sup>∂</sup> <sup>f</sup> ∂uτ*<sup>2</sup> = [*cij*]. For a complete list of the matrix elements, please consult [7].

#### *2.3. Stability Analysis of the Equilibrium Point E*<sup>1</sup>

Using the matrices defined above, the characteristic equation corresponding to *E*<sup>1</sup> is:

$$
\lambda \left(\lambda - a\_{22}\right) \left(\lambda - a\_{33}\right) \left(\lambda - a\_{55}\right) \left(\lambda - a\_{66}\right) \left(\lambda - a\_{77}\right) \left(\lambda - a\_{11} - b\_{11}e^{-\lambda \tau\_1}\right) = 0.
$$

We notice that *λ* = 0 is a root, and so, we are in a critical case for the stability of the nonlinear system (1).

#### 2.3.1. The Real Solutions of the Characteristic Equation

The real solutions of the characteristic equation are given by:

$$
\lambda\_1 = 0$$

$$
\lambda\_2 = a\_{22} = -\gamma\_2 < 0$$

$$
\lambda\_3 = a\_{33} = -k < 0$$

$$
\lambda\_4 = a\_{55} = -b\_1 < 0$$

$$
\lambda\_5 = a\_{66} = -\varepsilon\_1 - \frac{c\_1(1-\varepsilon\_2)c\_2}{(c\_2 + \hat{\imath}\_7)^2} - \frac{m\_2c\_2m\_1}{(m\_1 + \hat{\imath}\_7)^2} < 0$$

$$
\lambda\_7 = a\_{77} = -c\_3 < 0.$$

The existence of a zero eigenvalue implies a critical case for stability by the first approximation.

#### 2.3.2. Analysis of the Critical Case

The following theorem, developed and proven by some of the authors in [8], gives stability criteria in the critical case of a zero eigenvalue.

**Theorem 1** ([8], Theorem 2.1)**.** *Consider the following nonlinear system with time delays:*

$$\begin{aligned} \dot{\mathbf{x}}(t) &= A\_0 \mathbf{x}(t) + \sum\_{j=1}^{m} A\_j \mathbf{x}(t - \tau\_j) + \mathbf{F}[\mathbf{x}(t), \mathbf{x}(t - \tau\_1), \dots, \mathbf{x}(t - \tau\_m), y(t)] \\ \dot{y}(t) &= G[\mathbf{x}(t), \mathbf{x}(t - \tau\_1), \dots, \mathbf{x}(t - \tau\_m), y(t)], \end{aligned} \tag{2}$$

*where Aj* ∈ *Mn*(R)*, τ<sup>j</sup>* > 0 *for all* 1 ≤ *j* ≤ *m*, *G*(0, 0, ..., 0, *y*) = *F*(0, 0, ..., 0, *y*) = 0*,* ∀*y* ∈ R*, F takes values in* R*n, and G is scalar. F and G contain only powers of the variables with the sum greater than or equal to two. Then, for every δ* > 0*, there exist M*1(*δ*) *and M*2(*δ*) *with lim δ*→0 *M*1(*δ*) = *lim δ*→0 *M*2(*δ*) = 0 *so that, whenever* ||*x*(*t*)|| ≤ *δ*, ||*x*(*t* − *τj*)|| ≤ *δ*, 1 ≤ *j* ≤ *m*, |*y*| ≤ *δ,*

$$\begin{split} \left|| F(\mathbf{x}(t), \mathbf{x}(t-\tau\_1), \dots, \mathbf{x}(t-\tau\_m), y(t)) \right|| &\leq \\ \leq M\_1(\delta) (||\mathbf{x}(t)|| + ||\mathbf{x}(t-\tau\_1)|| + \dots + ||\mathbf{x}(t-\tau\_m)||) \\ &| G((\mathbf{x}(t), \mathbf{x}(t-\tau\_1), \dots, \mathbf{x}(t-\tau\_m), y(t))) \leq \\ \leq M\_2(\delta) (||\mathbf{x}(t)|| + ||\mathbf{x}(t-\tau\_1)|| + \dots + ||\mathbf{x}(t-\tau\_m)||). \end{split} \tag{3}$$

*Now, consider the following system with time delays. Suppose that the linear system:*

$$\dot{\mathbf{x}}(t) = A\_0 \mathbf{x}(t) + \sum\_{j=1}^{m} A\_j \mathbf{x}(t - \tau\_j) \tag{4}$$

*is asymptotically stable, that is, if λ is a root of the characteristic equation, then Re*(*λ*) < 0*. Then, the zero solution of (2) is simple stable, and if ϕ is the initial data of (2) in C* - [−*τ*, 0]; <sup>R</sup>*n*+<sup>1</sup> *with τ* = max *τ<sup>j</sup>* 1≤*j*≤*m , there exist δ* > 0 *so that, if* sup{||*ϕ*(*t*)||2/*t* ∈ [−*τ*, 0]} < *δ, then:*

$$\lim\_{t \to \infty} \mathfrak{x}\_i(t) = 0, i = 1, \dots, n \text{ and } \exists \lim\_{t \to \infty} \mathfrak{y}(t) = \mathfrak{y}.$$

In what follows, we transform our system so that Theorem 1 can be applied. We first perform a translation to zero: *yi* = *ui* − *u*ˆ*i*, for *i* = 3, 7.

The new system becomes:

*<sup>y</sup>*˙ = *fi*(*y*, *<sup>y</sup><sup>τ</sup><sup>j</sup>* ), *i* = 1, 7, *j* = 1, 2 (5) *<sup>y</sup>*˙1 <sup>=</sup> <sup>−</sup> *<sup>γ</sup>*<sup>0</sup> <sup>1</sup> + (*y*<sup>3</sup> <sup>+</sup> *<sup>u</sup>*ˆ3)*<sup>α</sup> <sup>y</sup>*<sup>1</sup> <sup>−</sup> *<sup>R</sup>*˜ *mu*<sup>7</sup> *R*˜ <sup>50</sup> + *y*<sup>7</sup> + *u*ˆ7 *y*<sup>1</sup> − (*η*1*<sup>e</sup>* + *η*2*e*)*ke*(*y*<sup>3</sup> + *u*ˆ3)*y*<sup>1</sup> −(1 − *η*1*<sup>e</sup>* − *η*2*e*)*βe*(*y*1, *y*<sup>3</sup> + *u*ˆ3)*y*<sup>1</sup> +2(*y*<sup>4</sup> + *u*ˆ4)(1 − *η*1*<sup>e</sup>* − *η*2*e*)*βe*(*y*1*τ<sup>l</sup>* , *y*3*τ<sup>l</sup>* + *u*ˆ3*<sup>τ</sup><sup>l</sup>* )*y*1*τ<sup>l</sup>* +*η*1*e*(*y*<sup>4</sup> + *u*ˆ4)*ke*(*y*3*τ*<sup>1</sup> + *u*ˆ3)*y*1*τ*<sup>1</sup> *<sup>y</sup>*˙2 <sup>=</sup> <sup>−</sup>*γ*2*y*<sup>2</sup> <sup>+</sup> *<sup>A</sup>*˜*eke*(*y*3*τ*<sup>2</sup> <sup>+</sup> *<sup>u</sup>*ˆ3*<sup>τ</sup>*<sup>2</sup> )*y*1*τ*<sup>2</sup> *<sup>y</sup>*˙3 <sup>=</sup> <sup>−</sup>*k*(*y*<sup>3</sup> <sup>+</sup> *<sup>u</sup>*ˆ3) + *<sup>a</sup>*<sup>1</sup> 1 + *y<sup>r</sup>* 2 *y*˙4 = (*y*<sup>4</sup> + *u*ˆ4) <sup>−</sup> *<sup>γ</sup>*<sup>0</sup> <sup>1</sup> + (*y*<sup>3</sup> <sup>+</sup> *<sup>u</sup>*ˆ3)*<sup>α</sup>* <sup>−</sup> *<sup>R</sup>*˜ *<sup>m</sup>*(*y*<sup>7</sup> <sup>+</sup> *<sup>u</sup>*ˆ7) *R*˜ <sup>50</sup> + *y*<sup>7</sup> + *u*ˆ7 <sup>+</sup> *<sup>γ</sup>*<sup>0</sup> 1 + (*y*3*τ*<sup>1</sup> + *u*ˆ3)*<sup>α</sup>* + *<sup>R</sup>*˜ *<sup>m</sup>*(*y*7*τ*<sup>1</sup> + *<sup>u</sup>*ˆ7) *<sup>R</sup>*˜ <sup>50</sup> + *<sup>y</sup>*7*τ*<sup>1</sup> + *<sup>u</sup>*ˆ7 *y*˙5 = −*b*1(*y*<sup>5</sup> + *u*ˆ5) + *a*<sup>2</sup> *<sup>y</sup>*˙6 <sup>=</sup> *<sup>b</sup>*1(*y*<sup>5</sup> <sup>+</sup> *<sup>u</sup>*ˆ5) <sup>−</sup> *<sup>e</sup>*1(*y*<sup>6</sup> <sup>+</sup> *<sup>u</sup>*ˆ6) <sup>−</sup> *<sup>c</sup>*1(<sup>1</sup> <sup>−</sup> *<sup>e</sup>*2) *c*<sup>2</sup> + *y*<sup>6</sup> + *u*ˆ6 (*y*<sup>6</sup> + *u*ˆ6)− *m*2*e*<sup>2</sup> *m*<sup>1</sup> + *y*<sup>6</sup> + *u*ˆ6 (*y*<sup>6</sup> + *u*ˆ6) *<sup>y</sup>*˙7 <sup>=</sup> *<sup>b</sup>*2*c*1(<sup>1</sup> <sup>−</sup> *<sup>e</sup>*2) *c*<sup>2</sup> + *y*<sup>6</sup> + *u*ˆ6 (*y*<sup>6</sup> + *u*ˆ6) − *e*3(*y*<sup>7</sup> + *u*ˆ7),

where we consider *y*˙4 = *g*(*y*3, *y*4, *y*7, *y*3*τ*<sup>1</sup> , *y*7*τ*<sup>1</sup> ), with *g*(0) = 0.

The matrices of partial derivatives into zero are, as before, *<sup>A</sup>*˜ <sup>=</sup> *<sup>∂</sup> <sup>f</sup> <sup>∂</sup><sup>y</sup>* = [*a*˜*ij*], *<sup>B</sup>*˜ <sup>=</sup> *<sup>∂</sup> <sup>f</sup> ∂yτ*<sup>1</sup> = [˜ *bij*], *<sup>C</sup>*˜ <sup>=</sup> *<sup>∂</sup> <sup>f</sup> ∂yτ*<sup>2</sup> = [*c*˜*ij*].

The characteristic equation for the zero solution of the new system (5) has exactly the same form as the one for *E*1. The terms that are different are:

$$\begin{array}{ll} \frac{\partial \mathcal{g}}{\partial y\_3}(0) = & \frac{\hat{\mathfrak{a}}\_4 \gamma\_0 \mathfrak{a} \hat{\mathfrak{a}}\_3^{\mathfrak{a}-1}}{(1+\hat{\mathfrak{a}}\_3^{\mathfrak{a}})^2} \\\\ \frac{\partial \mathcal{g}}{\partial y\_7}(0) = & -\frac{\hat{\mathcal{R}}\_m \mathcal{R}\_{50}}{(\hat{\mathcal{R}}\_{50}+\hat{\mathfrak{a}}\_7)^2} \\\\ \frac{\partial \mathcal{g}}{\partial y\_{3\tau\_1}}(0) = & -\frac{\hat{\mathfrak{a}}\_4 \gamma\_0 \mathfrak{a} \hat{\mathfrak{a}}\_3^{\mathfrak{a}-1}}{(1+\hat{\mathfrak{a}}\_3^{\mathfrak{a}})^2} \\\\ \frac{\partial \mathcal{g}}{\partial y\_{7\tau\_1}}(0) = & \frac{\hat{\mathcal{R}}\_m \mathcal{R}\_{50}}{(\mathcal{R}\_{S0}+\hat{\mathfrak{a}}\_7)^2} \end{array}$$

Notice that Theorem 1 is not directly applicable since the linear part is not equal to zero. In order to apply Theorem 1, we need to rewrite the system (5). Take *<sup>η</sup>* <sup>=</sup> *<sup>α</sup>*1*y*<sup>1</sup> <sup>+</sup> ··· <sup>+</sup> *<sup>α</sup>*7*y*<sup>7</sup> where *<sup>y</sup>*˙ <sup>=</sup> *Ay*˜ . This means that:

$$
\begin{bmatrix}
\dot{y}\_1\\ \dot{y}\_2\\ \dot{y}\_3\\ \dot{y}\_4\\ \dot{y}\_5\\ \dot{y}\_6\\ \dot{y}\_7
\end{bmatrix} = \begin{bmatrix}
\
\ddot{a}\_{11} & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & \ddot{a}\_{22} & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & \ddot{a}\_{33} & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & \ddot{a}\_{43} & 0 & 0 & 0 & \ddot{a}\_{47} \\ 0 & 0 & 0 & 0 & 0 & \ddot{a}\_{55} & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & \ddot{a}\_{65} & \ddot{a}\_{66} & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & \ddot{a}\_{76} & \ddot{a}\_{77}
\end{bmatrix} \begin{bmatrix}
y\_1\\ y\_2\\ y\_3\\ y\_4\\ y\_5\\ y\_6
\end{bmatrix}
$$

so,


We have *η*˙ = *α*1*y*˙1 + ··· + *α*7*y*˙7. Then,

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣

$$\begin{array}{rcl} \eta &=& \alpha\_1 \vec{a}\_{11} y\_1 + \alpha\_2 \vec{a}\_{22} y\_2 + \alpha\_3 \vec{a}\_{33} y\_3 + \alpha\_4 (\vec{a}\_{43} y\_3 + \vec{a}\_{47} y\_7) + \\ & \alpha\_5 \vec{a}\_{55} y\_5 + \alpha\_6 (\vec{a}\_{65} y\_5 + \vec{a}\_{66} y\_6) + \alpha\_7 (\vec{a}\_{76} y\_6 + \vec{a}\_{77} y\_7). \end{array}$$

Next, by forcing *η*˙ = 0, it follows that:

$$\begin{aligned} 0 &= \begin{array}{rcl} \alpha\_1 \overline{a}\_{11} y\_1 + \alpha\_2 \overline{a}\_{22} y\_2 + (\alpha\_3 \overline{a}\_{33} + \alpha\_4 \overline{a}\_{43}) y\_3 + (\alpha\_5 \overline{a}\_{55} + \alpha\_6 \overline{a}\_{65}) y\_5 \\ + (\alpha\_6 \overline{a}\_{66} + \alpha\_7 \overline{a}\_{76}) y\_6 + (\alpha\_4 \overline{a}\_{47} + \alpha\_7 \overline{a}\_{77}) y\_7, \end{array} \end{aligned}$$

from which we obtain:

$$\begin{aligned} \mathfrak{a}\_1 &= 0, \; \mathfrak{a}\_2 = 0 \\\\ \mathfrak{a}\_3 \vec{a}\_{33} + \mathfrak{a}\_4 \vec{a}\_{43} &= 0 \\\\ \mathfrak{a}\_5 \vec{a}\_{55} + \mathfrak{a}\_6 \vec{a}\_{65} &= 0 \\\\ \mathfrak{a}\_6 \vec{a}\_{66} + \mathfrak{a}\_7 \vec{a}\_{76} &= 0 \\\\ \mathfrak{a}\_4 \vec{a}\_{47} + \mathfrak{a}\_7 \vec{a}\_{77} &= 0. \end{aligned}$$

Thus, we have *<sup>α</sup>*<sup>4</sup> <sup>=</sup> 1, *<sup>α</sup>*<sup>3</sup> <sup>=</sup> <sup>−</sup>*a*˜43 *a*˜33 , *<sup>α</sup>*<sup>7</sup> <sup>=</sup> <sup>−</sup>*a*˜47 *a*˜77 , *<sup>α</sup>*<sup>6</sup> <sup>=</sup> *<sup>a</sup>*˜76*a*˜47 *a*˜66*a*˜77 and *<sup>α</sup>*<sup>5</sup> <sup>=</sup> <sup>−</sup>*a*˜65*a*˜76*a*˜47 *a*˜55*a*˜66*a*˜77 We remark that:

.

$$
\dot{y}\_{\text{3\tau}\_1} = \dot{a}\_{\text{33}} y\_{\text{3\tau}\_1} + R\_{\text{3\tau}\_1}
$$

*y*˙7*<sup>τ</sup>*<sup>1</sup> = *a*˜77*y*7*τ*<sup>1</sup> + *R*7*τ*<sup>1</sup> ,

with *R*3*τ*<sup>1</sup> and *R*7*τ*<sup>1</sup> containing terms of order higher than or equal to two. We introduce:

$$\begin{split} \eta\_{1} &= \quad a\_{3}y\_{3} + y\_{4} + a\_{5}y\_{5} + a\_{6}y\_{6} + a\_{7}y\_{7} - \frac{\tilde{b}\_{43}}{\tilde{a}\_{33}}y\_{3\tau\_{1}} - \frac{\tilde{b}\_{47}}{\tilde{a}\_{77}}y\_{7\tau\_{1}} \\ \text{Then, we have:}\\ \eta\_{1} &= \huge{\bar{b}}\_{43}y\_{3\tau\_{1}} + \huge{\bar{b}}\_{47}y\_{7\tau\_{1}} - \frac{\tilde{b}\_{43}}{\tilde{a}\_{33}}y\_{3\tau\_{1}} - \frac{\tilde{b}\_{47}}{\tilde{a}\_{77}}y\_{7\tau\_{1}} + R\_{4}^{1} \\ &= \huge{\bar{b}}\_{43}y\_{3\tau\_{1}} + \huge{\bar{b}}\_{47}y\_{7\tau\_{1}} - \frac{\tilde{b}\_{43}}{\tilde{a}\_{33}}(\tilde{a}\_{33}y\_{3\tau\_{1}} + R\_{3\tau\_{1}}) - \frac{\tilde{b}\_{47}}{\tilde{a}\_{77}}(\tilde{a}\_{77}y\_{7\tau\_{1}} + R\_{7\tau\_{1}}) \\ &= R\_{4}^{(2)}(y, y\_{\tau\_{1}}). \end{split}$$

We replace the fourth equation in (5) with the equation of *η*˙1 from above, so that the equation has a zero linear part.

Next, we consider:

*<sup>y</sup>*<sup>4</sup> <sup>=</sup> *<sup>η</sup>*<sup>1</sup> <sup>−</sup> *<sup>α</sup>*3*y*<sup>3</sup> <sup>−</sup> *<sup>α</sup>*5*y*<sup>5</sup> <sup>−</sup> *<sup>α</sup>*6*y*<sup>6</sup> <sup>−</sup> *<sup>α</sup>*7*y*<sup>7</sup> <sup>+</sup> ˜ *b*<sup>43</sup> *a*˜33 *<sup>y</sup>*3*τ*<sup>1</sup> <sup>+</sup> ˜ *b*<sup>47</sup> *a*˜77 *y*7*τ*<sup>1</sup> and we substitute *y*<sup>4</sup> in the equations of (5). The new system becomes: *<sup>y</sup>*˙1 <sup>=</sup> <sup>−</sup> *<sup>γ</sup>*<sup>0</sup> <sup>1</sup> + (*y*<sup>3</sup> <sup>+</sup> *<sup>u</sup>*ˆ3)*<sup>α</sup> <sup>y</sup>*<sup>1</sup> <sup>−</sup> *<sup>R</sup>*˜ *<sup>m</sup>*(*y*<sup>7</sup> <sup>+</sup> *<sup>u</sup>*ˆ7) *R*˜ <sup>50</sup> + *y*<sup>7</sup> + *u*ˆ7 *y*<sup>1</sup> − (*η*1*<sup>l</sup>* + *η*2*l*)*ke*(*y*<sup>3</sup> + *u*ˆ3)*y*<sup>1</sup> −(1 − *η*1*<sup>e</sup>* − *η*2*e*)*βe*(*y*1, *y*<sup>3</sup> + *u*ˆ3)*y*1+ <sup>+</sup>2(*η*<sup>1</sup> <sup>−</sup> *<sup>α</sup>*3*y*<sup>3</sup> <sup>−</sup> *<sup>α</sup>*5*y*<sup>5</sup> <sup>−</sup> *<sup>α</sup>*6*y*<sup>6</sup> <sup>−</sup> *<sup>α</sup>*7*y*<sup>7</sup> <sup>+</sup> ˜ *b*<sup>43</sup> *a*˜33 *<sup>y</sup>*3*τ*<sup>1</sup> <sup>+</sup> ˜ *b*<sup>47</sup> *a*˜77 *y*7*τ*<sup>1</sup> ) · (1 − *η*1*<sup>e</sup>* − *η*2*e*)*βe*(*y*1*τ*<sup>1</sup> , *y*3*τ*<sup>1</sup> + *u*ˆ3)*y*1*τ*<sup>1</sup> <sup>+</sup>*η*1*e*(*η*<sup>1</sup> <sup>−</sup> *<sup>α</sup>*3*y*<sup>3</sup> <sup>−</sup> *<sup>α</sup>*5*y*<sup>5</sup> <sup>−</sup> *<sup>α</sup>*6*y*<sup>6</sup> <sup>−</sup> *<sup>α</sup>*7*y*<sup>7</sup> <sup>+</sup> ˜ *b*<sup>43</sup> *a*˜33 *<sup>y</sup>*3*τ*<sup>1</sup> <sup>+</sup> ˜ *b*<sup>47</sup> *a*˜77 *y*7*τ*<sup>1</sup> ) · *ke*(*y*3*τ*<sup>1</sup> + *u*ˆ3*y*1*τ*<sup>1</sup> ) *<sup>y</sup>*˙2 <sup>=</sup> <sup>−</sup>*γ*2*y*<sup>2</sup> <sup>+</sup> *<sup>A</sup>*˜*eke*(*y*3*τ*<sup>2</sup> <sup>+</sup> *<sup>u</sup>*ˆ3)*y*1*τ*<sup>2</sup> *<sup>y</sup>*˙3 <sup>=</sup> <sup>−</sup>*k*(*y*<sup>3</sup> <sup>+</sup> *<sup>u</sup>*ˆ3) + *<sup>a</sup>*˜1 1 + *y<sup>r</sup>* 2 *η*˙1 = ˜ *<sup>b</sup>*43*y*3*τ*<sup>1</sup> + ˜ *<sup>b</sup>*47*y*7*τ*<sup>1</sup> <sup>−</sup> ˜ *b*<sup>43</sup> *a*˜33 *<sup>y</sup>*˙3*<sup>τ</sup>*<sup>1</sup> <sup>−</sup> ˜ *b*<sup>47</sup> *a*˜77 *y*˙7*<sup>τ</sup>*<sup>1</sup> + *R*<sup>1</sup> 4 = ˜ *<sup>b</sup>*43*y*3*τ*<sup>1</sup> + ˜ *<sup>b</sup>*47*y*7*τ*<sup>1</sup> <sup>−</sup> ˜ *b*<sup>43</sup> *a*˜33 (*a*˜33*y*3*τ*<sup>1</sup> <sup>+</sup> *<sup>R</sup>*3*τ*<sup>1</sup> ) <sup>−</sup> ˜ *b*<sup>47</sup> *a*˜77 (*a*˜77*y*7*τ*<sup>1</sup> + *R*7*τ*<sup>1</sup> ) = *R*(2) <sup>4</sup> (*y*, *yτ*<sup>1</sup> ) *<sup>y</sup>*˙5 <sup>=</sup> <sup>−</sup>˜ *b*1(*y*<sup>5</sup> + *u*ˆ5) + *a*˜2 *y*˙6 = ˜ *<sup>b</sup>*1(*y*<sup>5</sup> <sup>+</sup> *<sup>u</sup>*ˆ5) <sup>−</sup> *<sup>e</sup>*1(*y*<sup>6</sup> <sup>+</sup> *<sup>u</sup>*ˆ6) <sup>−</sup> *<sup>c</sup>*1(<sup>1</sup> <sup>−</sup> *<sup>e</sup>*2) *c*<sup>2</sup> + *y*<sup>6</sup> + *u*ˆ6 (*y*<sup>6</sup> + *u*ˆ6)− <sup>−</sup> *<sup>m</sup>*2*e*<sup>2</sup> *m*<sup>1</sup> + *y*<sup>6</sup> + *u*ˆ6 (*y*<sup>6</sup> + *u*ˆ6) *<sup>y</sup>*˙7 <sup>=</sup> ˜ *b*2*c*1(1 − *e*2) *c*<sup>2</sup> + *y*<sup>6</sup> + *u*ˆ6 (*y*<sup>6</sup> + *u*ˆ6) − *e*3(*y*<sup>7</sup> + *u*ˆ7). Consider: *f*1(*y*1, *y*3, *y*4, *η*1, *y*6, *y*7, *y*1*τ*<sup>1</sup> , *y*3*τ*<sup>1</sup> , *y*4*τ*<sup>1</sup> , *y*7*τ*<sup>1</sup> ) = <sup>−</sup> *<sup>γ</sup>*<sup>0</sup> <sup>1</sup> + (*y*<sup>3</sup> <sup>+</sup> *<sup>u</sup>*ˆ3)*<sup>α</sup> <sup>y</sup>*<sup>1</sup> <sup>−</sup> *<sup>R</sup>*˜ *<sup>m</sup>*(*y*<sup>7</sup> <sup>+</sup> *<sup>u</sup>*ˆ7) *R*˜ <sup>50</sup> + *y*<sup>7</sup> + *u*ˆ7 *y*1− −(*η*1*<sup>e</sup>* + *η*2*e*)*ke*(*y*<sup>3</sup> + *u*ˆ3)*y*<sup>1</sup> − (1 − *η*1*<sup>e</sup>* − *η*2*e*)*βe*(*y*1, *y*<sup>3</sup> + *u*ˆ3)*y*1+

+ 2(1 − *η*1*<sup>e</sup>* − *η*2*e*)*βe*(*y*1*τ*<sup>1</sup> , *y*3*τ*<sup>1</sup> + *u*ˆ3)*y*1*τ*<sup>1</sup> + *η*1*eke*(*y*3*τ*<sup>1</sup> + *u*ˆ3)*y*1*τ*<sup>1</sup> *η*<sup>1</sup> + *B*1(*yτ*1) and remark that the linear part of *f*<sup>1</sup> does not contain *η*1.

Since the other equations, except that of *η*1, do not contain *η*<sup>1</sup> at all, we conclude that Theorem 1 can be applied. Therefore, the stability of *E*<sup>1</sup> depends on the study of the transcendental terms in its characteristic equation.

2.3.3. The Transcendental Part of the Characteristic Equation

Consider the equation:

$$
\lambda - a\_{11} - b\_{11}e^{-\lambda \tau \mathbf{q}} = 0. \tag{6}
$$

The stability analysis of Equation (6) is classical (see, for example, [16]).

**Proposition 1.** *Let p*<sup>1</sup> <sup>=</sup> *<sup>γ</sup>*<sup>0</sup> 1 + *u*ˆ*<sup>α</sup>* 3 + *R*˜ *mu*ˆ7 *R*˜ <sup>50</sup> + *u*ˆ7 *. Assume that the following condition is true:*

$$(\eta\_{1\varepsilon}\pounds\_4 - \eta\_{1\varepsilon} - \eta\_{2\varepsilon})k\_{\varepsilon}(\pounds\_3) < p\_1 + (1 - 2\imath\partial\_4)(1 - \eta\_{1\varepsilon} - \eta\_{2\varepsilon})\pounds\_{\varepsilon}(0, \pounds\_3). \tag{7}$$

*Then, Equation (6) is stable for τ*<sup>1</sup> = 0 *and remains stable for τ*<sup>1</sup> > 0*.*

**Proof.** For the equilibrium point *E*1, we have:

$$a\_{11} = -p\_1 - (\eta\_{1\varepsilon} + \eta\_{2\varepsilon})k\_{\varepsilon}(\mathfrak{A}\_3) - (1 - \eta\_{1\varepsilon} - \eta\_{2\varepsilon})\beta\_{\varepsilon}(0, \mathfrak{A}\_3),$$

$$b\_{11} = 2\mathfrak{A}\_4(1 - \eta\_{1\varepsilon} - \eta\_{2\varepsilon})\beta\_{\varepsilon}(0, \mathfrak{A}\_3) + \eta\_{1\varepsilon}\mathfrak{A}\_3k\_{\varepsilon}(z\_4).$$

If *τ*<sup>1</sup> = 0, Equation (6) becomes:

$$(\lambda + p\_1 + (1 - 2\p\_4)(1 - \eta\_{1\varepsilon} - \eta\_{2\varepsilon})\beta\_\varepsilon(0, \mathbb{A}\_3) - (\eta\_{1\varepsilon}\p\_4 - \eta\_{1\varepsilon} - \eta\_{2\varepsilon})k\_\varepsilon(\mathbb{A}\_3) = 0$$

Equation (6) is stable for *τ*<sup>1</sup> = 0 if:

$$(\eta\_{1\varepsilon}\p\_4 - \eta\_{1\varepsilon} - \eta\_{2\varepsilon})k\_{\varepsilon}(\p\_3) < p\_1 + (1 - 2\p\_4)(1 - \eta\_{1\varepsilon} - \eta\_{2\varepsilon})\beta\_{\varepsilon}(0, \p\_3).$$

Since *b*<sup>11</sup> > 0, the following conditions from [16,17] must hold for stability when *τ*<sup>1</sup> > 0 :

$$1. \ a\_{11} < \frac{1}{\pi\_1};$$

2. *a*<sup>11</sup> + *b*<sup>11</sup> < 0.

Since *a*<sup>11</sup> = −*p*<sup>1</sup> − (*η*1*<sup>e</sup>* + *η*2*e*)*ke*(*u*ˆ3) − (1 − *η*1*<sup>e</sup>* − *η*2*e*)*βe*(0, *u*ˆ3) < 0 < 1 *τ*1 , the first condition holds true.

For the second condition we must have:

$$-p\_1 + (2\pounds\_4 - 1)(1 - \eta\_{1e} - \eta\_{2e})\pounds\_e(0, \pounds\_3) + (\eta\_{1e}\pounds\_4 - \eta\_{1e} - \eta\_{2e})k\_e(\pounds\_3) < 0\tag{8}$$

We remark that, if Condition (7) holds**,** Condition (8) will also hold.

**Remark 1.** *The equilibrium point E*<sup>1</sup> *is stable if Condition (7) holds.*

*2.4. Stability Analysis of the Equilibrium Point E*<sup>2</sup>

The characteristic equation corresponding to *E*<sup>2</sup> becomes:

$$\det(\lambda I - A - B e^{-\lambda \tau\_1} - \mathbb{C} e^{-\lambda \tau\_2}) = d\_1(\lambda) d\_2(\lambda),$$

where:

$$d\_1(\lambda) = (\lambda - a\_{77})(\lambda - a\_{66})(\lambda - a\_{55})$$

and:

$$\begin{split} d\_{2}(\lambda) &= -\lambda^{4} - \lambda^{3}(a\_{11} + a\_{22} + a\_{33} + b\_{11}e^{-\lambda\tau\_{1}}) - \\ &- \lambda^{2}[-a\_{11}a\_{22} - a\_{11}a\_{33} - a\_{22}a\_{33} - (a\_{22}b\_{11} + b\_{11}a\_{33})e^{-\lambda\tau\_{1}} + a\_{32}c\_{23}e^{-\lambda\tau\_{2}}] - \\ &- \lambda \left[ a\_{11}a\_{22}a\_{33} + a\_{22}b\_{11}a\_{33}e^{-\lambda\tau\_{1}} - (a\_{11}a\_{32}c\_{23} - a\_{13}a\_{32}c\_{21})e^{-\lambda\tau\_{2}} - \right. \\ &- (b\_{11}a\_{32}c\_{23} - a\_{32}b\_{13}c\_{21})e^{-\lambda\tau\_{1}}e^{-\lambda\tau\_{2}} \\ &- a\_{14}a\_{32}a\_{43}c\_{21}e^{-\lambda\tau\_{2}} - a\_{14}a\_{32}c\_{21}b\_{43}e^{-\lambda\tau\_{1}}e^{-\lambda\tau\_{2}}. \end{split}$$

**Remark 2.** *Since the characteristic equation corresponding to E*<sup>2</sup> *is complicated, the stability of E*<sup>2</sup> *is investigated using numerical procedures.*

#### *2.5. Numerical Simulations*

For the numerical simulations, we considered the time scale as days and consulted the available literature in order to set the values of the parameters. Table 1 shows the numerical values of the parameters presented in the erythropoiesis model (1).

**Table 1.** Parameter values for the erythropoiesis model.


Figure 1 depicts the equilibrium point *E*<sup>1</sup> in blue and a small perturbation from this equilibrium point in red. We obtained parameter conditions for the stability of this equilibrium point, but it would be preferable that it be unstable; remember that *E*<sup>1</sup> corresponds to the death of the patient.

We proved that the equilibrium point *E*<sup>1</sup> is stable if Condition (7) holds. For the parameter values in the table above, it is easy to verify that (7) does not hold. This can also be seen in Figure 1, where the evolution seems to be a favorable one for the patient.

Figure 2 depicts the equilibrium point *E*<sup>2</sup> in blue and a small perturbation from this equilibrium point in red. Although we could not obtain any parameter conditions for stability, it is easy to notice that the equilibrium point is unstable for the considered set of parameters. As *E*<sup>2</sup> corresponds to the chronic phase of the disease, the evolution obtained in Figure 2 is a desired one. The healthy cells grow in number, and the leukemic cells die out; this basically means that the patient recovers.

**Figure 2.** Stability of equilibrium point *E*2.

#### **3. The Leukopoiesis Model**

*3.1. The Mathematical Model*

In this model, *s*<sup>1</sup> represents the concentration of short-term stem-like white blood cells' precursors and *s*<sup>2</sup> the adult leukocytes. The treatment is presented through the function *l*<sup>1</sup> (see [15]):

$$l\_1(\mathbf{x\_6}) = \frac{\mathbf{x\_6}}{L\_{1S0} + \mathbf{x\_6}}$$

The equations that describe these variables mirror those in the erythropoiesis model, as the evolution of the leukocytes has many points in common. The difference between the two is the influence erythropoietin has on the production of the cells.

A non-constant rate of the elimination of stem cells due to treatment is encountered, and this leads to the consideration of a new, auxiliary variable *s*3.

Again, the last three equations depict the progress of the drug after administration. The following model describes a compartment of leukopoiesis coupled with the dynamics of 6-MP used in the maintenance therapy:

$$\dot{s} = \dot{f}\_1(s, s\_{\gamma\_j}), i = \overline{1,6}, j = \overline{3,4} \tag{9}$$

$$\dot{s}\_1 = -\gamma\_{11}s\_1 - T\_1l\_1(s\_6)s\_1 - \eta\_{11}k\_1(s\_2)s\_1 - \eta\_{21}k\_1(s\_2)s\_1 - (1 - \eta\_{11} - \eta\_{21})\beta\_1(s\_1)s\_1$$

$$\quad + 2e^{-\gamma\_{11}\tau\_{13}}s\_3(1 - \eta\_{11} - \eta\_{21})\beta\_1(s\_{1\tau\_1})s\_{1\tau\_1} + \eta\_{11}e^{-\gamma\_{11}\tau\_{13}}s\_3k\_1(s\_{2\tau\_1})s\_{1\tau\_1}$$

$$\begin{aligned} \dot{s}\_2 &= -\gamma\_{21}s\_2 + \bar{A}\_1k\_1(s\_{2\tau\_2})s\_{1\tau\_2} \\ \dot{s}\_3 &= -s\_3T\_1[l\_1(s\_{6\tau\_1}) - l\_1(s\_6)] \\ \dot{s}\_4 &= -b\_1s\_4 + a\_2 \\ \dot{s}\_5 &= -b\_1s\_4 - c\_1s\_5 - \frac{c\_1(1 - c\_2)}{c\_2 + s\_5}s\_5 - \frac{m\_2c\_2}{m\_1 + s\_5}s\_5 \\ \dot{s}\_6 &= -\frac{b\_2c\_1(1 - c\_2)}{c\_2 + s\_5}s\_5 - c\_3s\_6. \end{aligned}$$

#### *3.2. The Equilibrium Points and Linearization*

Following [7], we concluded that there are two types of equilibrium points corresponding to the model of leukopoiesis:

$$\begin{array}{l} \mathcal{E}\_1 = (0, 0, \mathfrak{s}\_{3\prime} \mathfrak{s}\_{4\prime} \mathfrak{s}\_{5\prime} \mathfrak{s}\_6) \\\\ \mathcal{E}\_2 = (\mathfrak{s}\_{1\prime} \mathfrak{s}\_{2\prime} \mathfrak{s}\_{3\prime} \mathfrak{s}\_{4\prime} \mathfrak{s}\_{5\prime} \mathfrak{s}\_6) .\end{array}$$

When linearizing System (9), the matrix of partial derivatives with respect to undelayed variables is denoted as *<sup>A</sup>*˜ <sup>=</sup> *<sup>∂</sup> <sup>f</sup> <sup>∂</sup><sup>s</sup>* = [*a*˜*ij*] and the matrices of the partial derivatives with respect to the delayed variables are denoted as *<sup>B</sup>*˜ <sup>=</sup> *<sup>∂</sup> <sup>f</sup> ∂sτ*<sup>3</sup> = [˜ *bij*] and *<sup>C</sup>*˜ <sup>=</sup> *<sup>∂</sup> <sup>f</sup> ∂sτ*<sup>4</sup> = [*c*˜*ij*]. For a complete list of the coefficients, please consult [7].

#### *3.3. Stability Analysis for Equilibrium Point E*˜ 1

The characteristic equation corresponding to *E*˜ <sup>1</sup> = (0, 0,*s*ˆ3,*s*ˆ4,*s*ˆ5,*s*ˆ6) is:

$$
\lambda (\lambda - a\_{22}) (\lambda - a\_{44}) (\lambda - a\_{55}) (\lambda - a\_{66}) \left(\lambda - a\_{11} - b\_{11} e^{-\lambda \tau\_1} \right) = 0.
$$

The real solutions of the characteristic equation corresponding to *E*˜ <sup>1</sup> are given by:

$$
\lambda\_1 = 0$$

$$
\lambda\_2 = a\_{22} = -\gamma\_{2I} < 0$$

$$
\lambda\_3 = a\_{44} = -b\_1 < 0$$

$$
\lambda\_4 = a\_{55} = -c\_1 - \frac{c\_1(1-c\_2)c\_2}{\left(c\_2 + \hat{s}\_5\right)^2} - \frac{m\_2 c\_2 m\_1}{\left(m\_1 + \hat{s}\_5\right)^2} < 0$$

$$
\lambda\_5 = a\_{66} = -c\_3 < 0.$$

The analysis of the critical case follows the same steps as in the previous case. As before, we can rewrite the system in a way that Theorem 1 can be applied. In order to avoid redundancy, we skip the computations. We concluded that the stability of *E*˜ <sup>1</sup> relies on the zeros of the transcendental equation:

$$
\lambda - a\_{11} - b\_{11}e^{-\lambda \tau \mathbf{q}} = 0.\tag{10}
$$

**Proposition 2.** *Assume that the following condition is true:*

$$(2\mathfrak{s}\_3 - 1)(1 - \eta\_{1l} - \eta\_{2l})\mathfrak{k}\_l(0) + (\mathfrak{s}\_3\eta\_{1l} - \eta\_{1l} - \eta\_{2l})\mathfrak{k}\_l(0) < \gamma\_{1l} + T\_1l\_1(\mathfrak{s}\_6). \tag{11}$$

*Then, Equation (10) is stable for τ*<sup>1</sup> = 0*, and it remains stable for all τ*<sup>1</sup> > 0*.*

The proof of Proposition 2 is similar to that of Proposition 1.

**Remark 3.** *The equilibrium point E*˜ <sup>1</sup> *is stable if Condition (11) holds.*

#### *3.4. Stability Analysis of the Equilibrium Point E*˜ 2

The characteristic equation corresponding to *E*˜ <sup>2</sup> is:

$$\begin{array}{l} \lambda \left(\lambda - a\_{44}\right) \left(\lambda - a\_{66}\right) \left(\lambda - a\_{55}\right) \\ \cdot \left[\left(\lambda - a\_{11} - b\_{11}\varepsilon^{-\lambda\tau\_{1}}\right) \left(\lambda - a\_{22} - c\_{22}\varepsilon^{-\lambda\tau\_{2}}\right) - c\_{21}\varepsilon^{-\lambda\tau\_{2}} \left(a\_{12} + b\_{12}\varepsilon^{-\lambda\tau\_{1}}\right)\right] = 0. \end{array}$$

#### 3.4.1. The Real Solutions of the Characteristic Equation

The nonzero real solutions of the characteristic equation corresponding to *E*˜ <sup>2</sup> are given by:

$$\begin{aligned} \lambda &= a\_{44} < 0 \\\\ \lambda &= a\_{66} < 0 \\\\ \lambda &= a\_{55} < 0 \end{aligned}$$

and

3.4.2. The Transcendental Part of the Characteristic Equation

The stability of the characteristic equation corresponding to *E*˜ <sup>2</sup> depends on the stability of the following equation:

$$(\lambda - a\_{11} - b\_{11}e^{-\lambda \tau\_1})(\lambda - a\_{22} - c\_{22}e^{-\lambda \tau\_2}) - c\_{21}e^{-\lambda \tau\_2}(a\_{12} + b\_{12}e^{-\lambda \tau\_1}) = 0\tag{12}$$

The stability analysis of Equation (12) follows the approach in [19], Theorem (1).

**Proposition 3.** *Assume the following conditions hold:*

$$\begin{aligned} a\_{11} + b\_{11} + a\_{22} + c\_{22} &< 0\\ a\_{11}a\_{22} + a\_{22}b\_{11} + a\_{11}c\_{22} + b\_{11}c\_{22} - a\_{12}c\_{21} - b\_{12}c\_{21} &> 0 \end{aligned} \tag{13}$$

*then Equation (12) is stable for τ*<sup>1</sup> = *τ*<sup>2</sup> = 0.

**Proof.** For *τ*<sup>1</sup> = *τ*<sup>2</sup> = 0, Equation (12) becomes:

$$\begin{aligned} \lambda^2 + \lambda(-a\_{11} - a\_{22} - b\_{11} - c\_{22}) + a\_{11}a\_{22} + a\_{22}b\_{11} \\ + a\_{11}c\_{22} + b\_{11}c\_{22} - a\_{12}c\_{21} - b\_{12}c\_{21} = 0 \end{aligned} \tag{14}$$

Suppose that *λ*<sup>1</sup> and *λ*<sup>2</sup> are two roots for Equation (14), then in order to have *Re*(*λ*1),*Re*(*λ*2) < 0, the following conditions must hold:

$$\begin{aligned} a\_{11} + b\_{11} + a\_{22} + c\_{22} &< 0\\ a\_{11}a\_{22} + a\_{22}b\_{11} + a\_{11}c\_{22} + b\_{11}c\_{22} - a\_{12}c\_{21} - b\_{12}c\_{21} &> 0 \end{aligned}$$

Now, let:

$$\begin{aligned} u\_1 &= a\_{11} + b\_{11} + a\_{22} \\ v\_1 &= a\_{11}a\_{22} + a\_{22}b\_{11} = a\_{22}(a\_{11} + b\_{11}) \\ u\_2 &= -c\_{22} \\ v\_2 &= a\_{11}c\_{22} + b\_{11}c\_{22} - a\_{12}c\_{21} - b\_{12}c\_{21} = c\_{22}(a\_{11} + b\_{11}) - c\_{21}(a\_{12} + b\_{12}) \end{aligned}$$

**Proposition 4.** *Consider the following conditions:*

$$\left(\mu\_1^2 - 2v\_1 - v\_2^2\right)^2 - 4\left(v\_1^2 - \mu\_2^2\right) > 0\tag{15}$$

$$
u\_1^2 - 2v\_1 - v\_2^2 < 0.\tag{16}$$

*If either Condition (15) or Condition (16) does not hold and if Equation (12) is stable for τ*<sup>1</sup> = *τ*<sup>2</sup> = 0, *it will remain stable for τ*<sup>1</sup> = 0 *and τ*<sup>2</sup> > 0.

**Proof.** Suppose *τ*<sup>1</sup> = 0 and *τ*<sup>2</sup> > 0. Equation (12) becomes:

$$
\lambda^2 - \mu\_1 \lambda + \upsilon\_1 + e^{-\lambda \tau\_2} (\mu\_2 \lambda + \upsilon\_2) = 0 \tag{17}
$$

The study of this equation follows the approach of Theorem 1 in [19]. Define:

$$\begin{aligned} P(\lambda) &= \lambda^2 - \mu\_1 \lambda + \upsilon\_1 \\ Q(\lambda) &= \mu\_2 \lambda + \upsilon\_2 \end{aligned}$$

Note that Conditions (i), (iv), and (v) of Theorem 1 in [19] are most likely to hold and Conditions (ii) and (iii) of Theorem 1 in [19] hold.

Then, the stability of Equation (17) depends on the zeros of the equation:

$$|P(iy)|^2 - |Q(iy)|^2 = 0\tag{18}$$

If Equation (18) has no positive roots (*y* > 0) and if (17) is stable with *τ*<sup>2</sup> = 0, it will be stable for all *τ*<sup>2</sup> > 0. If Equation (18) has a least one *y* > 0 as a root and all the roots are simple, as *τ*<sup>2</sup> increases, there might be stability switches. Therefore, if Equation (17) is stable at *τ*<sup>2</sup> = 0, it might become unstable when *τ*<sup>2</sup> = *τ*<sup>∗</sup> 2 .

Now, consider *P*(*iy*) = *PR*(*y*) + *iPI*(*y*) and *Q*(*iy*) = *QR*(*y*) + *iQI*(*y*) where *PR*, *PI*, *QR*, and *QI* are real-valued functions.

Equation (18) can be written as:

$$\left|P(iy)\right|^2 - \left|Q(iy)\right|^2 = 0$$

$$\left|P(iy)\right|^2 = \left|Q(iy)\right|^2$$

$$\left|P\_R^2(y) + P\_I^2(y) = Q\_R^2(y) + Q\_I^2(y).$$

We obtain the following polynomial:

$$y^4 + y^2(u\_1^2 - 2v\_1 - v\_2^2) + v\_1^2 - u\_2^2 = 0,\tag{19}$$

By setting *α* = *y*2, we have:

$$
\mu^2 + a(\mu\_1^2 - 2\upsilon\_1 - \upsilon\_2^2) + \upsilon\_1^2 - \mu\_2^2 = 0. \tag{20}
$$

In order that Equation (20) has at least one simple root *y* > 0, the following conditions must hold:

$$\left(u\_1^2 - 2v\_1 - v\_2^2\right)^2 - 4\left(v\_1^2 - u\_2^2\right) > 0\tag{21}$$

$$u\_1^2 - 2v\_1 - v\_2^2 < 0\tag{22}$$

Therefore, Equation (12) is stable if at least one of the conditions (21) or (22) is not met.

Now, we consider *τ*<sup>1</sup> = *τ*<sup>∗</sup> <sup>1</sup> (fixed) and *τ*<sup>2</sup> > 0. Equation (12) becomes:

$$(\lambda - a\_{11} - b\_{11}e^{-\lambda \tau\_1^\*})(\lambda - a\_{22} - c\_{22}e^{-\lambda \tau\_2}) - c\_{21}e^{-\lambda \tau\_2}(a\_{12} + b\_{12}e^{-\lambda \tau\_1^\*}) = 0\tag{23}$$

Equation (23) can be written as:

$$P(\lambda) + Q(\lambda)e^{-\lambda \overline{\nu}\_2} = 0,$$

with

$$\begin{aligned} P(\lambda) &= \lambda^2 - (a\_{11} + a\_{22})\lambda + a\_{11}a\_{22} - (b\_{11}\lambda + a\_{22}b\_{11})e^{-\lambda \tau\_1^\*} \\ Q(\lambda) &= -c\_{22}\lambda + a\_{11}c\_{22} - a\_{12}c\_{21} + (b\_{11}c\_{22} - b\_{12}c\_{21})e^{-\lambda \tau\_1^\*} .\end{aligned}$$

**Remark 4.** *Suppose that the equation:*

$$\left|P(iy)\right|^2 - \left|Q(iy)\right|^2 = 0$$

*has no positive real roots. Then, if Equation (23) is stable for τ*<sup>1</sup> = *τ*<sup>2</sup> = 0*, it will remain stable for τ*<sup>1</sup> = *τ*<sup>∗</sup> <sup>1</sup> *and all τ*<sup>2</sup> > 0.

Since *P*(*λ*) and *Q*(*λ*) are analytic functions, we can apply the result of Theorem 1 in [19]. Set *λ* = *iy*. We are interested in the roots of the equation:

$$F(y) = \left|P(iy)\right|^2 - \left|Q(iy)\right|^2.$$

We have:

$$F(y) = \begin{array}{c} y^4 + 2b\_{11}y^3\sin\left(y\tau\_1^\*\right) + r\_1y^2 + 2a\_{11}b\_{11}y^2\cos\left(y\tau\_1^\*\right) \\ + r\_2y\sin\left(y\tau\_1^\*\right) + r\_3\cos\left(y\tau\_1^\*\right) + k\_4 \end{array} \tag{24}$$

$$\begin{array}{l} r\_1 = & a\_{11}^2 + a\_{22}^2 + b\_{11}^2 - c\_{22}^2 \\ r\_2 = & 2a\_{22}^2 b\_{11} - 2b\_{11}c\_{22}^2 + 2b\_{12}c\_{21}c\_{22} \\ r\_3 = & 2a\_{11}b\_{12}c\_{21}c\_{22} + 2a\_{11}a\_{22}^2b\_{11} - 2a\_{11}b\_{11}c\_{22}^2 - 2a\_{12}b\_{12}c\_{21}^2 + 2a\_{12}b\_{11}c\_{21}c\_{22} \\ r\_4 = & a\_{11}^2a\_{22}^2 + a\_{22}^2b\_{11}^2 - a\_{11}^2c\_{22}^2 - a\_{12}^2c\_{21}^2 - b\_{11}^2c\_{22}^2 - b\_{12}^2c\_{21}^2 + 2a\_{11}a\_{12}c\_{21}c\_{22} + \\ & + 2b\_{11}b\_{12}c\_{21}c\_{22} \end{array}$$

If the function (24) has no *y* > 0 as a root and if Equation (12) is stable with *τ*<sup>1</sup> = *τ*<sup>2</sup> = 0, it will remain stable for all *τ*<sup>2</sup> > 0 and *τ*<sup>1</sup> = *τ*<sup>∗</sup> 1 .

#### **Remark 5.** *The equilibrium point E*˜ <sup>2</sup> *is stable if Equations (17) and (23) are stable.*

#### *3.5. Numerical Simulations*

Again, we considered the time scale as days and consulted the available literature in order to set the values of the parameters. Table 2 shows the numerical values of the parameters presented in the leukopoiesis model (9).


**Table 2.** Parameter values for the leukopoiesis model.

Due to the fact that the drug dynamics after administration is similar to the one obtained in the erythropoiesis model, we only focus on the representation of the white blood cells' precursors, *s*1, and the adult leukocytes, *s*2.

We would prefer that the equilibrium point *E*˜ <sup>1</sup> be unstable and perhaps attracted to a healthy state. Figure 3 depicts the equilibrium point *E*˜ <sup>1</sup> in blue and a small perturbation from this equilibrium point in red. The evolution of both the white blood cells' precursors and adult leukocytes is a favorable one for the patient.

In Figure 3, the equilibrium point *E*˜ <sup>1</sup> is clearly unstable. By checking Condition (11), we obtain the same result: the condition does not hold, and the equilibrium point is unstable.

**Figure 3.** Stability of equilibrium point *E*˜ 1: (**a**) Evolution of the white blood cells' precursors; (**b**) Evolution of the adult leukocytes

The equilibrium point *E*˜ <sup>2</sup> corresponds to a healthy state of the patient. Figure 4 shows the equilibrium point *E*˜ <sup>2</sup> in blue and a small perturbation from this equilibrium point in red. The equilibrium point is clearly unstable, but the evolution of the patient is still promising. The white blood cells' precursors and adult leukocytes both show an increase in number.

Oscillatory trajectories are common in blood cell evolution, as there is a natural inner oscillatory dynamics.

**Figure 4.** Stability of equilibrium point *E*˜ 2: (**a**) Evolution of the white blood cells' precursors; (**b**) Evolution of the adult leukocytes

#### **4. Conclusions**

The objective of this article was to present the stability study of two mathematical models that describe the processes of erythropoiesis and leukopoiesis (which are responsible for the production of red and white blood cells) in the case of maintenance therapy for acute lymphoblastic leukemia. The models were developed by some of the authors and introduced in [7].

The stability of all the equilibrium points was thoroughly investigated. When possible, parameter conditions for stability were determined. Some provocative and interesting conclusions resulted from the study of critical cases. For this purpose, a theorem designed by some of the authors and introduced in [8] was used. Numerical simulations were also used to validate and extend the study in some cases where the characteristic equation was too complicated.

From a biological and medical point of view, these models were validated first by the fact that the equilibrium points make sense biologically, as they correspond to existing states in which a patient may be found. Secondly, the numerical simulations, which were obtained using parameter values from the relevant literature (thus as close to reality as we could manage), depicted genuine developments of a potential patient's health.

In the case of the erythropoiesis model, the steady-states correspond to either the death of the patient or a chronic phase of the disease. As neither is the desired state of a patient, it would be preferable that these equilibrium points be unstable. The numerical analysis yielded this favorable result.

In the leukopoiesis model, the equilibrium points depicted either the death of the patient or a healthy state of the patient. For the considered configuration of the parameters, both these steady-states were unstable, but presented a positive evolution of the patient's condition.

Mathematical models that capture complicated situations in some diseases can aid in designing an adequate treatment. After a comprehensive mathematical study, these models can help to determine the correct dose of drug that needs to be administered for each individual patient.

**Author Contributions:** Conceptualization, I.B.; formal analysis, R.M.; methodology, I.B. and R.M.; resources, R.M.; software, A.-D.H.; validation, A.-D.H.; visualization, A.-D.H.; writing—original draft, I.B. and R.M.; writing—review and editing, I.B. and A.-D.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

