Stability and Hopf bifurcations in a competitive tumour-immune system with intrinsic recruitment delay and chemotherapy

: In this paper, a three-dimensional nonlinear delay di ﬀ erential system including Tumour cells, cytotoxic-T lymphocytes, T-helper cells is constructed to investigate the e ﬀ ects of intrinsic recruitment delay and chemotherapy. It is found that the introduction of chemotherapy and time delay can generate richer dynamics in tumor-immune system. In particular, there exists bistable phenomenon and the tumour cells would be cleared if the e ﬀ ect of chemotherapy on depletion of the tumour cells is strong enough or the side e ﬀ ect of chemotherapy on the hunting predator cells is under a threshold. It is also shown that a branch of stable periodic solutions bifurcates from the coexistence equilibrium when the intrinsic recruitment delay of tumor crosses the threshold which is new mechanism, which can help understand the short-term oscillations in tumour sizes as well as long-term tumour relapse. Numerical simulations are presented to illustrate that larger intrinsic recruitment delay of tumor leads to larger amplitude and longer period of the bifurcated periodic solution, which indicates that there exists longer relapse time and then contributes to the control of tumour growth.


Introduction
Cancer or malignant tumor results from uncontrolled growth of normal cells [1]. According to the report of World Health Organization, there are about 14.1 million new cancer cases and 8.2 million deaths in the world and it is expected that annual cancer cases will rise from 14 million in 2012 up to 22 million in the following two decades [2].
The immune system can monitor and protect the host from tumor evasion. However, the mechanisms of tumor evasion and immune response are still not completely understood. One of the adaptive immune responses is cell-mediated immunity(CMI). Cell-mediated immunity involves the production of cytotoxic T-lymphocytes (CTLs) and release of various cytokines in response to an antigen mediated by T-lymphocytes. CTLs are able to recognize and destroy tumor cells. The immune system can be classified into two subclasses, namely, the hunting cells (cytotoxic T lymphocytes) and the resting cells (T Helper cells). Most CTLs require cytokines from helper (resting) T-cells in order to be activated efficiently [3,4].
Different mathematical models have been constructed to understand the complex interaction between tumor and anti-tumor elements [5][6][7][8][9][10][11]. Kuznetsov et al. [12] proposed a classical tumor-immune model to explain the phenomena of "sneaking through" of the tumor and formation of a tumor "dormant state". Letellier et al. [13] found there exist chaotic attractor in a simple model of three competing cell populations (host, immune and tumour cells). Assia and Wang [14] investigated the effects of stochastic noises on tumors dynamics under treatment.
In [3], Sarkar et al. classified the immune system into two subclasses and proposed the following system (1.1) Here M, N and Z represent the densities of tumor cells, hunting predator cells (cytotoxic Tlymphocytes) and resting predator cells (T-Helper cells), respectively; q is the conversion of normal cells to malignant ones (fixed input), r 1 is the growth rate of tumor cells, k 1 and k 2 are the maximum capacity of tumor cells and resting cells, respectively; α is the predation/destruction rate of tumor cells by the hunting cells, β is the conversion rate of resting cells to hunting cells, d is the natural death rate of hunting cells. The authors obtained some thresholds to control the malignant tumor growth. In order to prolong the survival of patients it is essential to initiate a specific treatment regimen such as surgery, radiotherapy and chemotherapy [15][16][17][18][19][20]. Among various treatment methods, chemotherapy is very important and widely used in practice. However chemotherapy also has some side effects and it may kill the normal cells and immune cells [21,22]. Some mathematical models have been proposed to evaluate the effects of chemotherapy on tumor-immune dynamical behavior [23,24]. We also refer to [25] and the references therein for related works.
To derive more realistic models, recent papers [26][27][28][29][30][31][32][33] focus on delayed-induced oscillations in tumor-immune system dynamics to explain the short-term oscillations as well as long-term tumour relapse during the progress of the disease. Particularly, following [3], Subhas et al. [4] proposed the following tumor-immune competitive system with time delay (1.2) They neglect the fixed input term q of tumor in original Eq (1.1) by assuming that the tumour cells are not benign but malignant. Another modification to original Eq (1.1) is that a discrete time delay is added to prescribe the convention time from resting predator cells to hunting predator cells. It is found in Eq (1.2) that the introduction of time delay can lead to appearance of periodic solution. Some other tumor immune models with time delay were investigated in [34,35].
Considering the necessary times for molecule production and proliferation, the authors of [36][37][38][39] introduced the intrinsic recruitment delay of tumor. In fact, as a type of special population, the growth of tumor also need supply of nutrition and food and thus there exist intrinsic recruitment delay. In this paper, based on the models of Sarkar et al. [3], when considering the intrinsic recruitment delay and the effects of chemotherapy, we propose a new deterministic model describing tumor-immune responses as Here c 1 reflects the chemotherapy effect on the depletion rate of tumour cells, c 2 and c 3 represent the side effect of chemotherapy on the hunting predator cells and resting predator cells, respectively. The term r 1 M(1 − M(t−τ) k 1 ) describes the intrinsic recruitment of tumor and the positive time delay constant τ describes the necessary times for molecule production, proliferation, etc. We also modify the Eq (1.1) by adding the term −α 2 MN in the second equation, which represents the loss of hunting predator cells due to competition and destructive influence of tumour cells and is an important characteristic of cancer [12].
In order to simplify the analysis, we non-dimensionalize the Eq (1.3) by using the following rescaling After dropping the over bar notation for notational clarity, we obtain the following renormalization model (1.4) The rest of this paper is organized as follows. In Section 2, we obtain the non-negativity and boundeness of solutions of Eq (1.4). In Section 3, we study the existence and stability of possible steady states of the model. In Section 4, we study the changes in the stability of positive equilibrium though the Hopf bifurcation analysis. The direction and stability of the bifurcated periodic solution is also obtained. In Section 5, we carry out some numerical simulations and provide some biological interpretations. Conclusions and discussions are then presented in Section 6.

Non-negativity and boundeness
In this section, we discuss the non-negativity and boundedness of the solutions of Eq (1.4). For τ > 0, let C = C([−τ, 0], R 3 + ) denote the Banach space of continuous function mapping the interval [−τ, 0] into R 3 + with the topology of uniform convergence. The initial conditions are given by Theorem 2.1. The solutions of Eq (1.4) satisfying the initial conditions Eq (2.1) are non-negative and ultimately bounded.
Proof. We first define the right-hand side function of Eq (1.4) as It is obvious that the function G(t, K(t)) is locally Lipschitz and by the standard theory of functional differential equation, we know that there exists a unique solution for a given initial conditions. To prove the non-negativity of solutions Note that the solution is bounded for delay-logistic equation (refer to [41]). We know from Eq (2.2) and the comparison principle that x(t) is ultimately bounded for any nonnegative initial conditions. Define a new variable U(t) = y(t) + ρ γ z(t), and let d = min{k + m 2 , m 3 }. By the non-negativity of the solutions of Eq (1.4), we have

Existence and stability of equilibria
In this section, we investigate the existence and stability of equilibria for Eq (1.4) as delay is absent, that is Clearly, some of equilibria of Eq (3.1) can be obtained easily as follows: • There always exists extinction equilibrium E xyz 0 (0, 0, 0). • If m 1 < µ 1 , the Eq (3.1) exists immunity-free equilibrium E yz For the existence of coexistence equilibrium, we have the following theorem.
Proof. The coexistence equilibrium E * (x * , y * , z * ) of Eq (3.1) satisfies the following equations: From the first equation of (3.2), we have Substituting Eq (3.3) into the third equation of (3.2) comes to Substituting Eqs (3.3) and (3.4) into the second equation of (3.2) yields It is easy to see from Eqs (3.3) and (3.4) that y * > 0, z * > 0 if and only if It follows from Eq (3.5) that x * > 0 if and only if To study the stability of the equilibria above when τ = 0, we evaluate the Jacobian matrix of Theorem 3.2. If m 1 > µ 1 and m 3 > µ 2 , then E xyz 0 (0, 0, 0) is globally asymptotically stable. Proof. It follows from Eq (3.8) that the characteristic equation of (3.1) at extinction equilibrium It is easy to see that the eigenvalues are So if m 1 > µ 1 , and m 3 > µ 2 , equilibrium E xyz 0 (0, 0, 0) is locally asymptotically stable. Now we prove the global stability of E xyz 0 (0, 0, 0). Note that m 1 > µ 1 and m 3 > µ 2 . It follows that dx dt < 0, dz dt < 0, which implies that Then from Eq (3.9) and the second equation of (3.1) we have It follows from Eqs (3.9) and (3.10) that E xyz 0 (0, 0, 0) is globally attractive. The proof is completed.
Clearly the eigenvalues Then from Eq (3.11) and the second equation of (3.1), we have dz dt < 0 for large t, which indicates that There exists T 1 ( ) > 0 such that 0 ≤ y(t) ≤ as t > T 1 ( ). Then from the first equation of (3.1), we have It follow from Eq (3.13), comparison principle and the arbitrariness of that We know from Eqs (3.11), (3.12) and (3.14) that E yz 0 ( µ 1 −m 1 µ 1 , 0, 0) is globally attractive. The proof is completed.
Proof. It follows from Eq (3.8) that the characteristic equation of (3.1) at equilibrium E xy 0 is The eigenvalues are So if m 1 > µ 1 and µ 2 m 2 + ρm 3 > µ 2 (ρ − k), E xy 0 is locally asymptotically stable. Now we prove the global stability of E xy From the third equation of (3.1) we have Then from the second equation of (3.1) we have There exists T 2 ( ) > 0 such that 0 ≤ y(t) ≤ as t > T 2 ( ). It follows from the third equation of (3.1) It follow from Eqs (3.16) and (3.18), the arbitrariness of that Then from Eqs (3.15), (3.17) and (3.19), we know that E xy 0 (0, 0, µ 2 −m 3 µ 2 ) is globally attractive. The proof is completed.
is locally asymptotically stable.
Proof. It follows from Eq (3.8) that the characteristic equation of (3.1) at tumor-free equilibrium One of the eigenvalues is The other two eigenvalues λ 2 and λ 3 are determined by the following equation Then the real parts of eigenvalues λ 2 and λ 3 are negative. So if λ 1 < 0, that is the tumor-free equilibrium E x 0 (0, y 1 , z 1 ) is locally asymptotically stable. Remark 3.7. It can be seen from Theorem 3.5 and Theorem 3.6 that there exist bistable phenomenon in the Eq (3.1). Both E x 0 and E y 0 are locally asymptotically stable if m 1 < µ 1 , m 3 < µ 2 and which indicates that if the tumor can be cleared or not depends on the initial densities of different types of cells.
Theorem 3.8. Assume that Eqs (3.6) and (3.7) hold. The coexistence equilibrium E * is locally asymptotically stable if and only if the following conditions hold: Proof. It follows from Eq (3.8) that the characteristic equation of (3.1) at coexistence equilibrium E * is Clearly a 1 > 0. Then by Routh-Hurwitz Criterion, if a 3 > 0 and a 1 a 2 − a 3 > 0, that is the conditions (H1) and (H2) hold, E * is locally asymptotically stable.
Remark 3.9. It can be seen from Theorem 3.8 that if the destructive influence of tumor cells is small, that is α is small enough, the conditions (H1) and (H2) hold, and the coexistence equilibrium E * is locally asymptotically stable.
Theorem 3.10. If the coexistence equilibrium E * is asymptotically stable, all the other boundary equilibria are unstable.
Proof. If E * exist, we have µ 1 > m 1 , µ 2 > m 3 , which indicates that equilibria E xyz 0 , E xy 0 and E xy 0 can not be stable. It follows from Eq (3.6) that which indicates that E y 0 can not be stable. It follows from the stability of E * and Eq (3.7) that which indicates that E x 0 can not be stable. Remark 3.11. It is believed that, although it is difficulty to complete the proof, the coexistence equilibrium E * is globally asymptotically stable when it is locally asymptotically stable. However, we can obtain the result about uniform persistence of Eq (3.1), which indicates that three types of cells can coexist.
Proof. In order to prove the theorem we shall employ the method of average Lyapunov function [42]. Let us consider the following average Lyapunov function for the Eq (3.1) where ξ 11 , ξ 22 , ξ 33 are all nonnegative constants to be chosen suitably in the subsequent steps. It can be noted that ψ(x, y, z) is a positive continuous function defined in the positive quadrant R 3 + . Taking the time derivative along the solution of Eq (3.1), we have It is easy to show that the ω-limit sets for Eq (3.1) on the boundary of the positive cone consists of fixed points. Thus to prove the uniform persistence of the Eq (3.1), it is sufficient to verify that Γ(x, y, z) is positive for all singular points and the suitable choice of ξ 11 , ξ 22 , ξ 33 > 0, that is, the following conditions must be satisfied: and Equations (3.23) and (3.24) hold. The proof is completed.

Hopf bifurcation
In this section, we take the discrete time delay τ as a bifurcation parameter and show that when E * loses its stability and a Hopf bifurcation occurs when the time delay passes through a critical value.

Existence of Hopf bifurcation
In order to study the stability of the coexistence equilibrium E * of Eq (1.4), we first compute the Jacobian matrix as following Then the characteristic equation at E * is Putting λ = iω into the characteristic Eq (4.1) and separating the real and imaginary parts, we have Adding up the squares of the corresponding sides of Eqs (4.2) and (4.3) yields the following algebra equation with respect to ω ω 6 + p 1 ω 4 + p 2 ω 2 + p 3 = 0, (4.4) where Then Eq (4.4) becomes Q(u) ≡ u 3 + p 1 u 2 + p 2 u + p 3 = 0. (4.5) If Eq (4.5) has a positive real root u, the characteristic Eq (4.1) has a purely imaginary root iω = √ u; otherwise, Eq (4.1) has no purely imaginary root. Note that p 3 = b 2 3 − b 2 6 = (αµ 2 − µ 1 ργ)(αµ 2 + µ 1 ργ)x * 2 y * 2 z * 2 . We have p 3 < 0 if and only if (H1) is satisfied. We know that the stability of coexistence equilibrium E * as τ = 0 implies that Eq (4.5) has at least one positive real root. If we further assume that p 2 < 0, by Descartes sign rule, we know that Eq (4.5) has a unique positive real root u 0 . In addition it can be shown that Q (u 0 ) > 0.
We thus obtain the following result.

Direction and stability of the Hopf bifurcation
In the previous section, we studied the stability of the coexistence equilibrium E * of Eq (1.4) and the existence of Hopf bifurcations. We know Eq (1.4) has a branch of nonconstant periodic solutions bifurcating from E * near τ = τ 0 . It is interesting to further determine the direction, stability of bifurcated periodic solutions.

Then we can chooseD
Now following the standard notations and algorithms given in [44] and using a computation process similar to that of Wei and Ruan [43], we can obtain the following coefficients which will be used for determining the important qualities: where W (1) 11 (−1), W (i) 20 (0) and W (i) 11 (0) are given, respectively, by Furthermore the constant vector E 20 satisfy the following equation It follows that The constant vector E 20 satisfy the following equation It follows that So far, we have calculated g 20 , g 11 , g 02 , g 21 and then we can obtain ν 2 = − Re(c 1 (0)) Re(λ (τ 0 )) , β 2 = 2Re(c 1 (0)), It is known that ν 2 and β 2 will determine the direction and stability of the Hopf bifurcation, and T 2 determines the period of the bifurcated periodic solutions. In particular, if ν 2 > 0(ν 2 < 0), the Hopf bifurcation is supercritical (subcritical) and the bifurcated periodic solutions exist for τ > τ 0 (τ < τ 0 ). If β 2 < 0(β 2 > 0), the bifurcated periodic solutions are stable(unstable). The period will become longer(shorter) if T 2 > 0(T 2 < 0).

Numerical simulations
In this section, we carry out some numerical simulations to display the qualitative behaviours of Eq (1.4) and confirm the theoretical results obtained in the above sections. We perform simulations by using MATLAB software with realistic parameter values given in Table 1. According to the ranges of parameters in Table 1, we first take the following normalized parameter values Normalized parameters are given in Eq (5.1). One can note that when the initial density of hunting predator cells is small, the tumour cells tend to be persistent; when the initial density of hunting predator cells is relatively large, the tumour cells will be cleared. asymptotically stable and E * is unstable. Figure 1 presents a phase diagram and time series of the solutions of Eq (1.4). One can see that the trajectories converge to either E x 0 or E y 0 with different initial positions. Specially, when the initial density of hunting predator cells is small, the tumour cells tend to be persistent; when the initial density of hunting predator cells is relatively large, the tumour cells will be cleared.
To investigate the stability of coexistence equilibrium E * , we choose a slightly different set of the normalized parameter values as following The coexistence equilibrium is E * = (0.6, 0.07, 0.8). One can show that the conditions in Theorem 3.8 are satisfied and then E * is locally asymptotically stable when τ = 0, which implies that all the boundary equilibria including E x 0 and E y 0 are unstable. One can see that the trajectories around E * converge to it (see Figure 2).  From Eq (4.6), we can obtain the Hopf bifurcation point τ 0 = 14.5117. By Eq (4.8), we obtain the Hopf bifurcation parameters ν 2 = 40.3056, β 2 = −0.3561, T 2 = 5.8354. Then we know that the Hopf bifurcation at E * is supercritical, the bifurcated periodic solutions is stable and the period will increase with the increase of time delay. Figure 3 gives the time series and phase diagrams of the solutions of Equation (1.4) with different time delays. One can see that when the time delay is under the threshold, that is τ = 14 < τ 0 , E * is a stable focus(see Figure 3(a),(e)), which can explain the short-term oscillations in tumour sizes; when τ = 15, 20, 25 > τ 0 , a branch of stable periodic solutions bifurcate from E * and the period become larger with the increase of τ (see Figure 3(b)-(d), 3(f)-(h)), which can explain the long-term tumour relapse. Figure 4 presents the amplitude and period of the bifurcated periodic solutions of Eq (1.4) with respect to different time delay τ. One can see that, larger intrinsic recruitment delay lead to larger amplitude and longer period of the bifurcated periodic solutions, which indicates that the period of tumor relapse become longer and then can provide us with a longer period to control or delete the tumor. Furthermore, with the increase of time delay, the minimum of the periodic tumor tends to zero, which indicates larger delay τ contributes to the control of tumor growth. Figure 5 shows the sensitivity of state variable x(t) (tumour cells) to the parameters m 1 , m 2 , m 3 and α. The oscillation behaviour indicates that the tumour cells population is very sensitive in the early time intervals. Note that ∂x ∂m 1 < 0. We know that the increase of m 1 is beneficial to control and clear the tumor (see Figure 5(a)). On the other hand, the positivity of ∂x ∂m i 0 (i = 1, 2) indicates that stronger side effect of chemotherapy on hunting predator cells and resting predator cells leads to increase of tumor cells (see Figure 5  is not conductive to the control of tumor. Here we have used the so called "direct approach" to find sensitivity functions of Eq (1.4) (refer to [40]).

Conclusions
The long-term relapse phenomenon of tumor, that is the appearance of periodic solution in a system including tumor, is interesting for the control or deletion of tumor. In this paper, we construct a threedimensional nonlinear delay differential system including Tumor cells, cytotoxic-T lymphocytes, Thelper cells with recruitment delay and chemotherapy to investigate long-term relapse phenomena of tumor.
We first show that the existence and stability of the equilibria are influenced by the effects of chemotherapy and destructive influence of tumor cells and the results are summarized in the following Table 2. Table 2. Equilibria and their stability of Eq (1.4).
It is found that there exist bistable phenomenon in Eq (1.4). Both tumor-free equilibrium E x 0 and E y 0 are locally asymptotically stable at the same time, which indicates that if the tumour cells can be eradicated or not rely on the initial density of the hunting predator cells (see Figure 1). We point out  Figure 5. Sensitivity functions ∂x ∂m 1 , ∂x ∂m 2 , ∂x ∂m 3 and ∂x ∂α for Eq (1.4). Normalized parameters are given in Eq (5.2). One can see that the tumour cells population is very sensitive to the perturbation of the chosen set of parameters in the early time intervals and the sensitivity decreases by time. The increase of m 1 is beneficial to control and clear the tumor. However stronger side effect of chemotherapy and destructive influence of tumor cells on predator cells leads to increase of tumor cells.
that the bistable phenomenon appearing in this paper is a new dynamical behavior compared with the results in [27] and [28], which indicates that the introduction of chemotherapy and destructive influence of tumor cells can generate richer dynamics in tumor-immune system.
It is found that the coexistence equilibrium is locally asymptotically stable if the destructive influence of tumor cells is small enough (see Figure 2). We also show that the introduction of intrinsic recruitment delay of tumor can leads to the appearance of periodic solutions through Hopf bifurcation, which is a different factor that explain the long-term relapse phenomena of tumor compared with the factor in [27] and [28] (see Figure 3). We further investigate the direction and stability of bifurcated periodic solutions by applying normal form method and center manifold theory. Finally, by numerical simulations, it is shown that, with the increase of delay τ, the amplitude and period of the bifurcated periodic solutions increase and the minimum of the periodic tumor cells tend to zero, which indicates that the increase of delay τ is beneficial to the control of the growth of tumor (see Figure 4). The sensitivity of tumour cells population to small perturbations in m 1 , m 2 , m 3 and α are also investigated. It follows that tumour cells is negatively proportion with increasing the parameter m 1 and is positively proportion with increasing the parameter m 2 , m 3 and α and they are very sensitive in the early time intervals and the sensitivity decreases by time to be insensitive in the steady state (see Figure 5).
Some aspects of the problem remain to be examined and studied in the future. For instance, we plan to extend our analysis to global Hopf bifurcation analysis. Some other factors that influence the dynamics of tumor including the biological heterogeneity of tumor and the resistance to immunotherapy or chemotherapy may be investigated.