Effects of immune response and time delays in models of acute myeloid leukemia

In this paper, we propose a general acute myeloid leukemia (AML) model and introduce an immune response and time delays into this model to investigate their effects on the dynamics. Based on the existence, stability and local bifurcation of three types of equilibria, we show that the immune response is a best strategy for the control of the AML on the condition that the rates of proliferation and differentiation of the hematopoietic lineage exceed a threshold. In particular, a powerful immune response leads to a bistability feature meaning there exist the leukemia cells and healthy cells in the bone marrow or only the healthy cells. In addition, we further reveal that the time delays in the feedback regulation and immune response process induce a series of oscillations around the steady state, which shows that the leukemia cells are hardly eliminated. Our work in this paper aims to investigate the complex dynamics of this AML model with the immune response and time delays on the basis of mathematical models and numerical simulations, which may provide a theoretical guidance for the treatments of the AML.


Introduction
Acute myeloid leukemia (AML) is a blood cancer which is aroused by the malignant clone of hematopoietic stem cells or hematopoietic progenitor cells. Those anomalous cells mutate as leukemia original cells inhibit the normal proliferation and differentiation of hematopoietic cells and destroy the hematopoietic function of the bone marrow [1]. AML is a common myeloid disease with a great mortality rate and has attracted considerable attention in the clinical pathogenesis and treatment [2]. Based on a lot of experiments for the mouse, the leukemia stem cells have been found and defined [3,4], which shows a new way to the diagnosis and treatment of the AML [5,6]. With the development of molecular biology, there are more and more treatments of the AML, such as transplantation of bone marrow stem cells [7], targeted chemotherapy [8], combination chemotherapy [9] and resonance chemotherapy [10,11]. However, the treatments face a greater challenge because of the complicated interactions among the cells in the bone marrow. Therefore, it is of interest to model these interactions in the mathematical manner to provide a theoretical guidance to the AML [12][13][14].
There are two multistage lineages in the bone marrow, i.e., the hematopoietic lineage and leukemia lineage. The hematopoietic stem cell (HSC) of the hematopoietic lineage is multi-potent, which maintains the population by proliferation and differentiates into different mature cells. The hematopoietic progenitor cell (HPC) is mono-potent and has the ability to proliferate and differentiate into the functional cells (e.g., red blood cells and platelets). The leukemia stem cell (LSC) in the leukemia lineage has a similar feature with the HPC and proliferates into itself and differentiates into the leukemia cells. On the basis of cancer stem cells hypothesis [15], it is of significance to investigate the characteristics (e.g., proliferation, differentiation and death) of the cells in the hematopoietic lineage and the leukemia lineage, and interactions (e.g., competition and feedback) between them [16][17][18].
Presently mathematical models have been proposed to characterize the competition, feedback and immune response in the AML model [19]. Goardon et al. found that the LSC shares more characteristic with the HPC [20,21]. The LSC and HPC can proliferate and differentiate into only one type of cells by consuming the energy in the bone marrow, but they have to compete with each other, which results in the coupling between the hematopoietic lineage and leukemia lineage, further leads to complex dynamics of the system [22][23][24]. Rodriguez-Brenes et al. showed that a key of the carcinogenic process is that the proliferation and differentiation of the LSC break away the inhibition feedback from terminally differentiated cells in the leukemia lineage [25]. To control the output of multistage cell lineages, Lander et al. studied different types of feedback configurations, feedback sensitivities in the linages [26,27]. For the treatment and prevention of the AML, Crowell et al. regulated the feedback's strength among of the HSC, HPC and LSC [28,29]. To derive optimal control strategies, Sharp et al. introduced an immune function into the AML model and showed the effectiveness of treatment by improving the strength of immune response [30].
In biological systems, time delays play an important role in producing complex dynamics, small time delays often are harmless while larger delays lead to oscillations through a Hopf bifurcation [31][32][33][34]. Stépán firstly introduced a time delay into a predatorprey model and found that oscillations are induced by the delay [35]. Niculescu et al. studied the delayed effects on the asymptotic stability in modeling the post-transplantation dynamics of the immune response to chronic myelogenous leukemia [36]. In our previous work, the introduced delays induce a Hopf-born rotation in a time-delayed three-gene auto-regulated and mutually-repressed core genetic regulation net-work [37]. Inspired by the abovementioned discussions, we aim to investigate the dynamical behaviors of this AML model where the introduced delays characterize the lags between the feedback regulations in the hematopoietic lineage and the immune response in the leukemia lineage. We provide a lot of mathematical analysis of this model and do not focus on the biological relevance of delays in reality.
The remainder of this paper is organized as follows. In Sect. 2, we introduce an immune response and time delays into an AML model. Section 3 provides the analysis of the existence, stability and local bifurcation of three types of steady states for the non-delayed model. In Sect. 4, we study the counterparts of the model with time delays and compare them with the non-delayed model. Finally, conclusions are given in Sect. 5.

Crowell's model
Crowell et al. proposed a classical acute myeloid leukemia (AML) model [28], which is given by where S, A, D, L , T represent haematopoietic stem cells (HSCs), haematopoietic progenitor cells (HPCs), terminally differentiated blood cells, leukemia stem cells (LSCs) and fully differentiated leukemia cells, respectively. The quantities 0 < ρ S , ρ A , ρ L < 1 and 0 < δ S , δ A , δ L < 1 are the proliferation rates and differentiation rates of HSCs (S), HPCs (A) and LSCs (L). 0 < μ D < 1 and 0 < μ T < 1 are the death rates of terminally differentiated blood cells (D) and leukemia cells (T ). K 1 and K 2 are the population sizes of cells within the bone marrow which determine the carrying capacities. In the following studies, we scale the populations sizes K 1 = K 2 = 1. Z 1 ≡ S and Z 2 ≡ A + L are the total niche size, where A and L are coupled because of the competition between LSCs (L) and HPCs (A), which is resulted by the hypothesis that LSCs (L) and HPCs (A) occupy the same niche in the bone marrow [9,28,38].

Sharp's model
On the basis of the model (1), Sharp et al. introduced an immune response to apply the optimal control [30].
The feature of this model considers a small population of leukemia cells may be overcome by the healthy cells without the intervention [19]. Then, model (1) is written as where α > 0 and γ > 0 are the parameters controlling the intensity of immune response.

Jiao's model
To investigate the effects of feedback regulation, Jiao et al. introduced the Hill function to characterize the feedback regulation in the haematopoietic lineage [29]. The model is described in the form below: where 1+h i D n , i = 1, 2 and n is the Hill exponent. The feedback coefficients are g i , h i ≥ 0, i = 1, 2. p i0 and v i0 , (i = 1, 2, 3) are the maximal proliferation rates and division rates of S, A, L. Based on the biological significance and mathematical principles, it is assumed that 0 < p i0 < 1, v i0 > 0.

Our model involving both immune response and feedback regulation
Based on the models (2) and (3), we propose a new model with the immune response and feedback regulation to investigate the effect of intensity of the immune response on the dynamics. The model reads as 2.5 Our model introducing delays in the feedback regulation and immune response To investigate the effects of delays on this AML model, we introduce time delays in the feedback regulation and immune process. In details, τ 1 , τ 2 > 0 describe the time lags during feedback processes from terminally differentiated blood cells (D) to HSCs (S) and HPCs (A), respectively. The delay τ 3 > 0 shows the lag during the immune response of LSCs (L). Then, model (4) becomes A schematic diagram introducing the immune response and time delays is shown in Fig. 1; the definitions and basal values for all the parameters are given in Appendix B, see Table 1.

Effects of intensity of immune response
In this section, we discuss the existence, stability and local bifurcations of positive equilibria of the model (4) with the intensity of immune response.

Existence of the equilibria
It is obvious that the origin is a trivial equilibrium of the system (4), which means all the cells in the bone marrow die out. Similar with the analysis in Ref. [29], we divide the equilibria into three types, denoted by E l (0, 0, 0, where E l is a purely leukemic steady state with S = A = D = 0, which means that there are only LSCs (L) and terminally differentiated leukemia cells (T ). E h is a healthy steady state when there are only healthy cells in the bone marrow, i.e., L = T = 0. E c is a coexisting steady state Proof It is from Eq. (4) that The number of positive roots of Eq. (6) depends on the values of Δ, On the other hand, when Δ > 0 and α < p 30 v 30 (2 − 1 p 30 )γ or Δ = 0 and γ < (2 − 1 p 30 ), the system (4) has only one purely leukemic steady state E l (0, 0, 0, L * l1 , T * l1 ). As shown in Fig. 2, there is no purely leukemia steady state (III) under the condition α > p 30 v 30 (2 − 1/ p 30 )γ . By adjusting the appropriate immune response α and differentiation probability p 30 , there will not exist a purely leukemia steady state.
, and D * h is a root of the following equation: If p 10 > 0.5, it gives a 0 < 0 and F 0 (0) < 0. When a i > 0, there exists a positive root of Eq.
According to Eq. (7), the existence of E h is mainly determined by p 10 and p 20 , not related to the immune response (α) and the proliferation rate p 30 . When the rates p 10 and p 20 reach a threshold, the HPCs compete with the LSCs meaning that the increase of the immune response inhibits the leukemia stem cells.

Theorem 3 There are at most two coexisting steady states E c determined by the equations G(L
Proof From the system (4), the coexisting steady state E c satisfies the following equations:  (4) will have two or one or zero coexisting steady state The number of the coexisting steady state is same as the number of intersections between the curves H (L , D) = 0 and G(L , D) = 0.
There exists one or two coexisting steady states when the proliferation rate p 30 of the LSCs is bigger than the rate p 10 of the HSCs and the rate p 20 of the HPCs (see Fig. 3).

Stability of the equilibria
In this subsection, we discuss the stability of the equilibria E l , E h and E c and bifurcation mechanisms regarding the intensity parameter α of the immune response.
The Jacobian matrix J (S, A, D, L , T ) of the system (4) is given by where Substituting the origin in the Jacobian matrix (8), we obtain their eigenvalues regarding the origin, which is given by Obviously, the origin is a saddle when p 1 > 0.5 or p 2 > 0.5 or α < γ (2 p 30 − 1)v 30 .

Theorem 4
The purely leukemic steady state E l is locally stable when p 10 < 0.5, p 20 < Proof The Jacobian matrix at the state E l is Then, the eigenvalues are as follows: When As shown in Fig. 4, there are two purely leukemia steady states E l1 (0, 0, 0, 0.71295, 0.66542) and E l2 (0, 0, 0, 0.027051, 0.025247) when p 10 = 0.45, p 20 = 0.68, p 30 = 0.8. According to Theorem 4, it follows that E l1 is stable and E l2 is unstable. Since only the eigenvalue λ l4 is related to the intensity α, the stable equilibrium becomes unstable with the increasing values of the intensity.
According to the above analysis, it is an efficient way to treat the AML because increasing the immune response transforms the stable state E l into unstable, and further makes the orbits in the stable coexisting steady state E c or the healthy steady state E h .

Theorem 5
The steady states E c and E h are locally stable under the conditions H i > 0, i = 1, · · · , 5.

Proof The characteristic equation is
Based on the Routh-Hurwitz stability criterion, all the roots of Eq. (10) have negative real parts if and only if Therefore, the coexisting steady state E c and healthy steady state E h are locally stable when the conditions H i > 0, (i = 1, · · · , 5) hold.
From Fig. 5, there have no stable coexisting and healthy equilibria when the immune response is very weak (α → 0). With an increase of the immune response, there appears a bistability state, even a monostability state, which shows that the leukemia cells will be eliminated by a powerful immune response.

Local bifurcation
In this subsection, we discuss the effects of the intensity of immune response on the local bifurcation around three types of equilibria.

Transcritical bifurcation
A transcritical bifurcation is characterized by the exchange of stability of two equilibria. For this bifur- Then, its characteristic equation is Clearly, Eq. (11) has a simple zero eigenvalue at the threshold α [T h] = γ v 30 Therefore, the transversality conditions for the transcritical bifurcation are Based on Sotomayor's theorem for the transcritical bifurcation [39], we conclude that the system (4) experiences a transcritical bifurcation around E h at the threshold α [T h] . The proof of the bifurcation at E l is similar with the abovementioned proof, and it is omitted.
The system (4) has a stable healthy steady state E h (0.42047, 0.92976, 0.8898, 0, 0) with p 20 = 0.75, p 30 = 0.95 and μ D = 0.365. Then, the threshold is α [T h] = 0.000117, which is consistent with the numerical value in Fig. 6. After the transcritical bifurcation, the healthy steady state becomes stable, and the stable coexisting steady state E c loses its stability. In addition, the state E c has no biological significance because of the negative population of leukemia cells. In the same way, the theoretical value L(α [T l] ) = 2 − 1/ p 20 = 0.5294 with α = 0.0665 and p 10 = 0.45 approximates the numerical value L = 0.5294122 (see Fig. 7).
In Fig. 8, the coexisting steady state E c undergoes a Hopf bifurcation with p 10 = 0.51, p 20 = 0.85 and p 30 = 0.98, and induced an unstable limit cycle around E c .

Bistability of equilibria
Based on the results of the system (3), there exists a monostability. When this system is imposed by the immune response, there appears bistability, i.e., the healthy and coexisting equilibria E h and E c are stable. When the proliferation rate of haematopoietic cells is smaller than the rate of leukemia cells with a proper immune response, there will be a stable coexisting equilibrium and a stable healthy equilibrium (see Fig. 9).
When the ratio of cells of the hematopoietic lineage to the cells of the leukemia lineage in the bone marrow is small, the trajectory of the system (4) tends to a stable coexisting equilibrium (E c ). If the immune capacity α decreases and the hematopoietic function p 1 of the bone marrow is destroyed, the trajectory reaches a purely leukemia steady state E l in Fig. 7. On the other hand, the trajectory converges to a healthy equilibrium E h when the ratio of hematopoietic cells and leukemia cells in the bone marrow is large. In this case, increasing the immunity is an effective way to treat the AML.

Effects of time delays
In this section, we study transcritical and Hopf bifurcations of the delayed system (5), compare the corresponding bifurcations with the non-delayed system (4) and further explore the effects of time delays on the dynamics.

The characteristic equation of delayed model
The system (5) are rewritten in the matrix form, that is, Then, the characteristic equation of Eq. (5) is To obtain theoretical results for the time delays, we suppose τ 1 = τ 2 = τ ,τ 3 = 2τ . Then, the characteristic equation (12) is simplified as

Transcritical bifurcation
Firstly, we study the transcritical bifurcation around the purely leukemia steady state E l (0, 0, 0, L * l , T * l ). The characteristic equation (13) becomes . Specifically, the eigenvalues are Clearly, the eigenvalue λ l4 is related to the time delays, and other eigenvalues are same as the counterparts of non-delayed model. Therefore, it gives the following theorem.
Since the transcritical bifurcation takes place around the equilibrium with zero eigenvalues, the critical values of the bifurcation parameter α are not affected by the time delays. Comparing Figs. 10 and 6, Figs. 11 and 7, the delayed system (5) has the same critical values with the non-delayed system (4), which shows that the delays have no effects on the transcritical bifurcation.

Remark 1
The time delays have no effects on the transcritical bifurcation around three types of equilibria in the delayed system (5) with τ 1 = τ 2 = τ and τ 3 = 2τ .

Existence of Hopf bifurcation
In this section, we study the existence of Hopf bifurcation around the coexisting steady state E c and the healthy steady state E h . Suppose λ = ±iω (ω > 0) is a root of Eq. (13), then it follows When the equation H (v) = 0 has at least one positive root v * , the characteristic equation (13) has a pair of purely imaginary eigenvalues λ * = ±iω * with ω * = √ v * . Then, the critical delays for the Hopf bifurcation around the steady state E c are given by Taking the derivative of λ(τ ) with respect to τ in Eq. There is a Hopf bifurcation around E c and six number of Hopf bifurcations around E h , where the black dots denote a stable limit cycle, and the green, magenta, red, blue, yellow dots represent unstable limit cycles with one to five Floquet multipliers whose modulus is greater than 1 Then, the sign of the real part of Eq.
When the transversality condition T R F I −T I F R = 0 holds, the system (5) experiences a Hopf bifurcation at the time delays τ k around the coexisting steady state E c . Compared to the non-delayed system (4), the delays in the system (5) induce a series of Hopf bifurcations around the healthy steady state E h (see Fig. 12). The Hopf bifurcation induces a stable limit cycle which results in the oscillations of the leukemia cells (L) showing that it is impossible to eliminate all the leukemia cells (see Fig. 13). Due to the existence of time delays, the purely healthy steady state (L = T = 0) is just an ideal state because there always exist a small amount of leukemia stem cells and terminally differentiated leukemia cells in the bone marrow.

Conclusions
In the present study, we have proposed a new acute myeloid leukemia (AML) model incorporating with feedback regulation, immune response and time delays, which describes the competition between the haematopoietic progenitor cells and the leukemia stem cells. With the controlled bifurcation parameters, this model causes two bifurcation types-Hopf bifurcation and transcritical bifurcation. In the biological sense, the appearance of a Hopf bifurcation shows that the leukemia cells always exist. When the model undergoes a transcritical bifurcation, the concentration of the leukemia cells declines.
Regarding this model, we have shown that the existence and stability of three types of equilibria are determined by the proliferation rates and the strength of immune response. We also have found that enhancing immunity is an efficient way to prevent and treat the AML. With an increase of the immune response, the leukemia steady state loses its stability and the healthy steady state becomes stable at the different transcritical bifurcation points. When the immune response is introduced into this model, it is interesting that the model displays a bistability trait.