In this paper, a three-dimensional nonlinear delay differential system including Tumour cells, cytotoxic-T lymphocytes, T-helper cells is constructed to investigate the effects 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 effect of chemotherapy on depletion of the tumour cells is strong enough or the side effect 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.
Citation: Qingfeng Tang, Guohong Zhang. Stability and Hopf bifurcations in a competitive tumour-immune system with intrinsic recruitment delay and chemotherapy[J]. Mathematical Biosciences and Engineering, 2021, 18(3): 1941-1965. doi: 10.3934/mbe.2021101
Related Papers:
[1]
Dan Liu, Shigui Ruan, Deming Zhu .
Stable periodic oscillations in a two-stage cancer model of tumor and immune system interactions. Mathematical Biosciences and Engineering, 2012, 9(2): 347-368.
doi: 10.3934/mbe.2012.9.347
[2]
Yuting Ding, Gaoyang Liu, Yong An .
Stability and bifurcation analysis of a tumor-immune system with two delays and diffusion. Mathematical Biosciences and Engineering, 2022, 19(2): 1154-1173.
doi: 10.3934/mbe.2022053
[3]
Marek Bodnar, Monika Joanna Piotrowska, Urszula Foryś .
Gompertz model with delays and treatment: Mathematical analysis. Mathematical Biosciences and Engineering, 2013, 10(3): 551-563.
doi: 10.3934/mbe.2013.10.551
[4]
Jiawei Deng, Ping Jiang, Hongying Shu .
Viral infection dynamics with mitosis, intracellular delays and immune response. Mathematical Biosciences and Engineering, 2023, 20(2): 2937-2963.
doi: 10.3934/mbe.2023139
[5]
Huan Kong, Guohong Zhang, Kaifa Wang .
Stability and Hopf bifurcation in a virus model with self-proliferation and delayed activation of immune cells. Mathematical Biosciences and Engineering, 2020, 17(5): 4384-4405.
doi: 10.3934/mbe.2020242
[6]
Shishi Wang, Yuting Ding, Hongfan Lu, Silin Gong .
Stability and bifurcation analysis of $ SIQR $ for the COVID-19 epidemic model with time delay. Mathematical Biosciences and Engineering, 2021, 18(5): 5505-5524.
doi: 10.3934/mbe.2021278
[7]
Hongying Shu, Wanxiao Xu, Zenghui Hao .
Global dynamics of an immunosuppressive infection model with stage structure. Mathematical Biosciences and Engineering, 2020, 17(3): 2082-2102.
doi: 10.3934/mbe.2020111
[8]
Yu Yang, Gang Huang, Yueping Dong .
Stability and Hopf bifurcation of an HIV infection model with two time delays. Mathematical Biosciences and Engineering, 2023, 20(2): 1938-1959.
doi: 10.3934/mbe.2023089
[9]
Yuan Ma, Yunxian Dai .
Stability and Hopf bifurcation analysis of a fractional-order ring-hub structure neural network with delays under parameters delay feedback control. Mathematical Biosciences and Engineering, 2023, 20(11): 20093-20115.
doi: 10.3934/mbe.2023890
[10]
Zhichao Jiang, Xiaohua Bi, Tongqian Zhang, B.G. Sampath Aruna Pradeep .
Global Hopf bifurcation of a delayed phytoplankton-zooplankton system considering toxin producing effect and delay dependent coefficient. Mathematical Biosciences and Engineering, 2019, 16(5): 3807-3829.
doi: 10.3934/mbe.2019188
Abstract
In this paper, a three-dimensional nonlinear delay differential system including Tumour cells, cytotoxic-T lymphocytes, T-helper cells is constructed to investigate the effects 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 effect of chemotherapy on depletion of the tumour cells is strong enough or the side effect 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.
1.
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
Here M, N and Z represent the densities of tumor cells, hunting predator cells (cytotoxic T-lymphocytes) and resting predator cells (T-Helper cells), respectively; q is the conversion of normal cells to malignant ones (fixed input), r1 is the growth rate of tumor cells, k1 and k2 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
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 c1 reflects the chemotherapy effect on the depletion rate of tumour cells, c2 and c3 represent the side effect of chemotherapy on the hunting predator cells and resting predator cells, respectively. The term r1M(1−M(t−τ)k1) 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 −α2MN 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 re-scaling
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.
2.
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],R3+) denote the Banach space of continuous function mapping the interval [−τ,0] into R3+ with the topology of uniform convergence. The initial conditions are given by
x(ξ)=ϕ1(ξ),y(ξ)=ϕ2(ξ),z(ξ)=ϕ3(ξ)
(2.1)
with ϕi(ξ)≥0, ξ∈[−τ,0] and ϕi(0)>0(i=1,2,3).
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
where K(t)=(K1(t),K2(t),K3(t))T and K1(t)=x, K2(t)=y, K3(t)=z. 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 of Eq (1.4), we first prove the non-negativity over the time interval [0,τ]. Considering the right-hand side functions of Eq (1.4) over the time interval [0,τ], we have
It follows that the solutions of Eq (1.4) remain nonnegative in [0,τ]. Similarly, we can repeat the process over [τ,2τ] and so on by using the method of steps, it can be proved that for any finite interval [0,t], the solutions of Eq (1.4) remains non-negative.
Now we consider the boundeness of the solutions. It follows from the non-negativity of the solutions and the first equation of Eq (1.4) that
dxdt≤μ1x(1−x(t−τ)).
(2.2)
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+m2,m3}. By the non-negativity of the solutions of Eq (1.4), we have
Theorem 3.2.If m1>μ1 and m3>μ2, then Exyz0(0,0,0) is globally asymptotically stable.
Proof. It follows from Eq (3.8) that the characteristic equation of (3.1) at extinction equilibrium Exyz0(0,0,0) is
(λ−μ1+m1)(λ+k+m2)(λ−μ2+m3)=0.
It is easy to see that the eigenvalues are
λ1=μ1−m1,λ2=−k−m2<0,λ3=μ2−m3.
So if m1>μ1, and m3>μ2, equilibrium Exyz0(0,0,0) is locally asymptotically stable.
Now we prove the global stability of Exyz0(0,0,0). Note that m1>μ1 and m3>μ2. It follows that dxdt<0,dzdt<0, which implies that
limt→∞x(t)→0,limt→∞z(t)→0.
(3.9)
Then from Eq (3.9) and the second equation of (3.1) we have
limt→∞y(t)→0.
(3.10)
It follows from Eqs (3.9) and (3.10) that Exyz0(0,0,0) is globally attractive. The proof is completed.
Theorem 3.3.Assume that m1<μ1. If m3>μ2, then Eyz0(μ1−m1μ1,0,0) is globally asymptotically stable.
Proof. It follows from Eq (3.8) that the characteristic equation of (3.1) at Eyz0(μ1−m1μ1,0,0) is
(λ−m1+μ1)(λ+k+m2+α(μ1−m1)μ1)(λ−μ2+m3)=0.
Clearly the eigenvalues
λ1=m1−μ1<0,λ2=−k−m2−α(μ1−m1)μ1<0.
So if μ2<m3, we have λ3=μ2−m3<0. We know that Eyz0(μ1−m1μ1,0,0) is locally asymptotically stable.
Now we prove the global stability of Eyz0(μ1−m1μ1,0,0). It follows from m3>μ2 that dzdt<0, which implies that
z(t)→0ast→∞.
(3.11)
Then from Eq (3.11) and the second equation of (3.1), we have dzdt<0 for large t, which indicates that
y(t)→0ast→∞.
(3.12)
There exists T1(ϵ)>0 such that 0≤y(t)≤ϵ as t>T1(ϵ). Then from the first equation of (3.1), we have
(μ1−m1−ϵ)x−μ1x2≤dxdt≤(μ1−m1)x−μ1x2.
(3.13)
It follow from Eq (3.13), comparison principle and the arbitrariness of ϵ that
limt→∞x(t)=μ1−m1μ1.
(3.14)
We know from Eqs (3.11), (3.12) and (3.14) that Eyz0(μ1−m1μ1,0,0) is globally attractive. The proof is completed.
Theorem 3.4.Assume that m3<μ2. If m1>μ1 and μ2m2+ρm3>μ2(ρ−k), then Exy0(0,0,μ2−m3μ2) is globally asymptotically stable.
Proof. It follows from Eq (3.8) that the characteristic equation of (3.1) at equilibrium Exy0 is
(λ−m3+μ2)(λ−μ1+m1)(λ−ρμ2−ρm3−μ2k−μ2m2μ2)=0.
The eigenvalues are
λ1=m3−μ2,λ2=μ1−m1λ3=ρμ2−ρm3−μ2k−μ2m2μ2.
So if m1>μ1 and μ2m2+ρm3>μ2(ρ−k), Exy0 is locally asymptotically stable.
Now we prove the global stability of Exy0(0,0,μ2−m3μ2). It follows from m1>μ1 that dxdt<0, which implies that
limt→∞x(t)=0.
(3.15)
From the third equation of (3.1) we have
dzdt≤(μ2−m3)z−μ2z2,
which indicates that
lim supt→∞z(t)≤μ2−m3μ2.
(3.16)
Then from the second equation of (3.1) we have
dydt≤(ρz−k−m2)y≤(ρ(μ2−m3)μ2−k−m2)y
for large t. Since μ2m2+ρm3>μ2(ρ−k), we have dydt<0, which indicates that
limt→∞y(t)=0.
(3.17)
There exists T2(ϵ)>0 such that 0≤y(t)≤ϵ as t>T2(ϵ). It follows from the third equation of (3.1) that
dzdt≥(μ2−m3−ϵ)z−μ2z2,
which implies that
lim inft→∞z(t)≥μ2−m3−ϵμ2.
(3.18)
It follow from Eqs (3.16) and (3.18), the arbitrariness of ϵ that
limt↦∞z(t)=μ2−m3μ2.
(3.19)
Then from Eqs (3.15), (3.17) and (3.19), we know that Exy0(0,0,μ2−m3μ2) is globally attractive. The proof is completed.
Theorem 3.5.Assume that m1<μ1 and m3<μ2. If μ2m2+ρm3>μ2(ρ−k)+μ2αu1(m1−μ1), Ey0(μ1−m1μ1,0,μ2−m3μ2) is locally asymptotically stable.
Proof. It follow from Eq (3.8) that the characteristic equation of (3.1) at equilibrium Ey0 is
(λ−m3+μ2)(λ−m1+μ1)(λ−Aμ1μ2)=0,
where
A=μ2(ρ−k)+μ2αμ1(m1−μ1)−μ2m2−ρm3.
Note that eigenvalues
λ1=m3−μ2<0,λ2=m1−μ1<0
since m1<μ1 and m3<μ2. Thus if eigenvalue λ3=Aμ1μ2<0, that is
μ2m2+ρm3>μ2(ρ−k)+μ2αu1(m1−μ1),
the equilibrium Ey0 is locally asymptotically stable.
Theorem 3.6.Assume that μ2m2+ρm3<μ2(ρ−k). If μ2m2+ρm3<μ2(ρ−k)+ργ(m1−μ1), then Ex0(0,y1,z1) is locally asymptotically stable.
Proof. It follows from Eq (3.8) that the characteristic equation of (3.1) at tumor-free equilibrium Ex0(0,y1,z1) is
(λ−μ1+m1+y1)(λ2+μ2z1λ+γρy1z1)=0.
One of the eigenvalues is
λ1=μ1−m1−y1=μ1−m1−μ2ρ−μ2k−μ2m2−ρm3ργ.
The other two eigenvalues λ2 and λ3 are determined by the following equation
λ2+μ2z1λ+γρy1z1=0.
Note that
λ2+λ3=−μ2z1<0,λ2λ3=γρy1z1>0.
Then the real parts of eigenvalues λ2 and λ3 are negative. So if λ1<0, that is
μ2m2+ρm3<μ2(ρ−k)+ργ(m1−μ1),
the tumor-free equilibrium Ex0(0,y1,z1) 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 Ex0 and Ey0 are locally asymptotically stable if m1<μ1, m3<μ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:
Clearly a1>0. Then by Routh-Hurwitz Criterion, if a3>0 and a1a2−a3>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>m1, μ2>m3, which indicates that equilibria Exyz0, Exy0 and Exy0 can not be stable. It follows from Eq (3.6) that
which indicates that Ey0 can not be stable. It follows from the stability of E∗ and Eq (3.7) that
μ2m2+ρm3>μ2(ρ−k)+ργ(m1−μ1),
which indicates that Ex0 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.
Theorem 3.12.If the axial and planer singular points exist and
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)
ψ(x,y,z)=xξ11yξ22zξ33,
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 R3+. 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:
Γ(Exyz0)=ξ11(μ1−m1)+ξ22(−k−m2)+ξ33(μ2−m3)>0,
(3.20)
Γ(Eyz0)=ξ22(−k−m2−αμ1−m1μ1)+ξ33(μ2−m3)>0,
(3.21)
Γ(Exy0)=ξ11(μ1−m1)+ξ22(ρμ2−m3μ2−k−m2)>0,
(3.22)
Γ(Ey0)=ξ22(ρμ2−m3μ2−αμ1−m1μ1−k−m2)>0,
(3.23)
and
Γ(Ex0)=ξ11(μ1−m1−μ2ρ−μ2k−μ2m2−ρm3ργ)>0.
(3.24)
Noting μ1−m1>0 and μ2−m3>0, we can take large ξ11 and ξ33 and small ξ22 such that Eqs (3.20), (3.21) and (3.22) are satisfied. It follows that if
μ2(ρ−k)+ργ(m1−μ1)<μ2m2+ρm3<μ2(ρ−k)+μ2αμ1(m1−μ1),
Equations (3.23) and (3.24) hold. The proof is completed.
Remark 3.13. It is found that if the effect of chemotherapy on depletion of tumor cells is strong enough (m1>μ1), tumor cells would be cleared. While if the effect of chemotherapy is relatively weak (m1<μ1), the dynamics of Eq (3.1) rely on the side effects of chemotherapy on the hunting predator cells (m2) and resting predator cells (m3), respectively. Then to control or clear the tumor, it is important to increase the effect of chemotherapy on tumor cells and let the side effect of chemotherapy on the predator cells be under a threshold.
4.
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.
4.1. 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
Putting λ=iω into the characteristic Eq (4.1) and separating the real and imaginary parts, we have
ω3−b2ω=b5ωcos(ωτ)+(b4ω2−b6)sin(ωτ)
(4.2)
and
b1ω2−b3=b5ωsin(ωτ)+(−b4ω2+b6)cos(ωτ).
(4.3)
Adding up the squares of the corresponding sides of Eqs (4.2) and (4.3) yields the following algebra equation with respect to ω
ω6+p1ω4+p2ω2+p3=0,
(4.4)
where
p1=b21−2b2−b24,p2=b22−2b1b3+2b4b6−b25,p3=b23−b26.
Let u=ω2. Then Eq (4.4) becomes
Q(u)≡u3+p1u2+p2u+p3=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
p3=b23−b26=(αμ2−μ1ργ)(αμ2+μ1ργ)x∗2y∗2z∗2.
We have p3<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 p2<0, by Descartes sign rule, we know that Eq (4.5) has a unique positive real root u0. In addition it can be shown that Q′(u0)>0.
Let ω0=√u0. It follows from Eqs (4.2) and (4.3) that
Equation (4.1) has a pair of purely imaginary roots ±iω0.
Now, it remain to check the transversality condition required by the Hopf bifurcation theorem. Substituting λ(τ) into Eq (4.1) and differentiating the resulting equation in τ, we obtain
We know that at each τj, j=0,1,2,⋅⋅⋅, a pair of characteristic roots of Eq (4.1) cross the imaginary axis to the right. Then the transversality condition required by the Hopf bifurcation theorem is satisfied. We thus obtain the following result.
Theorem 4.1.Assume that the conditions in Theorem 3.8 hold.
(i) If p2>0, E∗ is locally asymptotically stable for τ∈[0,τ0] and unstable for τ>τ0.
(ii) Equation (1.4) undergoes a Hopf bifurcation at E∗ when τ=τ0. That is, Eq (1.4) has a branch of nonconstant periodic solutions bifurcating from E∗ near τ=τ0.
4.2. 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.
Let μ=τ−τ0 and then μ=0 is the Hopf bifurcation value of Eq (1.4) at E∗. Introduce the change of variables
x1(t)=x−x∗,x2(t)=y−y∗,x3(t)=z−z∗.
Let ˉxi(t)=xi(τt) and drop the bars for simplicity of notation. Define operators L(μ):C⟶R3 and F(μ,ϕ):R×C⟶R3 by
For φ∈C([0,1],(R3)∗) and ϕ∈C([−1,0],R3), we define a bilinear form as
<φ,ϕ>=ˉφT(0)ϕ(0)−∫0−1∫θ0ˉφT(ξ−θ)dη(θ)ϕ(ξ)dξ,
where η(θ)=η(θ,0). Then it is not difficult to show that A∗(0) is the adjoint operator of A(0). By the discussion at the beginning of this section, we know that ±iω0τ0 are eigenvalues of A(0). Thus, they are also eigenvalues of A∗(0).
Let q(θ)=(1,q2,q3)Teiω0τ0θ be the eigenvectors of matric A corresponding to eigenvalue iω0τ0. Then we have
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:
It is known that ν2 and β2 will determine the direction and stability of the Hopf bifurcation, and T2 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 T2>0(T2<0).
5.
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.
Then we obtain the boundary equilibria Ex0=(0,0.49,0.8750), Ey0=(0.95,0,0.9567). It is easy to check that the conditions in Theorems 3.5 and 3.6 are satisfied. Then both Ex0 and Ey0 are locally 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 Ex0 or Ey0 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.
Figure 1.
(a) Basins of attraction of Ex0 and Ey0. (b) (c)Trajectories are shown to converge to either Ex0 or Ey0 depending on their initial positions. The initial values in 1(b) is [x(0),y(0),z(0)]=[0.9,0.6,0.2] and in 1(c) is [x(0),y(0),z(0)]=[0.9,0.5,0.2]. 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.
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 Ex0 and Ey0 are unstable. One can see that the trajectories around E∗ converge to it (see Figure 2).
Figure 2.
The coexistence equilibrium E∗=(0.6,0.07,0.8) is locally asymptotically stable when time delay is absent. Normalized parameters are given in Eq (5.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, T2=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 3.
Time series and phase diagrams of the solutions of Eq (1.4) with different time delays τ. (a) (d):the positive equilibrium E∗=(0.6,0.07,0.8) of Eq (1.4) is local asymptotically stable when τ=14<τ0=14.5117. (b) (c) (e) (f): the bifurcated periodic solutions occur through Hopf bifurcations when τ=15,20,25>τ0. Normalized parameters are given in Equation (5.2).
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 4.
Amplitude and period of tumor cells population for the bifurcated periodic solutions of Eq (1.4) with respect to time delay τ. One can note that larger intrinsic recruitment delay lead to larger amplitude and longer period of the bifurcated periodic solution. Normalized parameters are given in Eq (5.2).
Figure 5 shows the sensitivity of state variable x(t) (tumour cells) to the parameters m1, m2, m3 and α. The oscillation behaviour indicates that the tumour cells population is very sensitive in the early time intervals. Note that ∂x∂m1<0. We know that the increase of m1 is beneficial to control and clear the tumor (see Figure 5(a)). On the other hand, the positivity of ∂x∂mi0 (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(b), (d)). The positivity of ∂x∂α implies that larger destructive influence of tumor cells 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]).
Figure 5.
Sensitivity functions ∂x∂m1, ∂x∂m2, ∂x∂m3 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 m1 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.
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 three-dimensional nonlinear delay differential system including Tumor cells, cytotoxic-T lymphocytes, T-helper 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).
Case
Conditions
Equilibria and stability
m1>μ1
m3>μ2
Exyz0 is GAS
m3<μ2, μ2m2+ρm3≥μ2(ρ−k)
Exy0 is GAS
μ2m2+ρm3<μ2(ρ−k)
Ex0 is LAS
m1<μ1
m3>μ2
Eyz0 is GAS
m3<μ2, μ2m2+ρm3>μ2(ρ−k)+μ2α(m1−μ1)μ1
Ey0 is LAS
μ2m2+ρm3<μ2(ρ−k)+ργ(m1−μ1)
Ex0 is LAS
Note: GAS means globally asymptotically stable; LAS means locally asymptotically stable.
It is found that there exist bistable phenomenon in Eq (1.4). Both tumor-free equilibrium Ex0 and Ey0 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 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 m1, m2, m3 and α are also investigated. It follows that tumour cells is negatively proportion with increasing the parameter m1 and is positively proportion with increasing the parameter m2, m3 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.
Acknowledgments
The authors thank the referees for their valuable comments and suggestions. This work is supported by the National Natural Science Foundation of P.R. China (11871403, 11701472), Fundamental Research Funds for the Central Universities (XDJK2020B050).
Conflict of interest
The authors declare no conflict of interest in this paper.
References
[1]
L. M. F. Merlo, J. W. Pepper, B. J. Reid, C. C. Maley, Cancer as an evolutionary and ecological process, Nat. Rev. Cancer, 6 (2006), 924–935. doi: 10.1038/nrc2013
[2]
A. Jemal, F. Bray, M. M. Center, Global cancer statistics, Ca. Cancer. J. Clin., 61 (2011), 69–90.
[3]
R. R. Sarkar, S. Banerjee, Cancer self remission and tumor stability-a stochastic approach, Math. Biosci., 196 (2005), 65–81. doi: 10.1016/j.mbs.2005.04.001
[4]
K. Subhas, J. J. Nieto, Mathematical modeling of tumor-immune competitive system, considering the role of time delay, Appl. Math. Comput., 340 (2019), 180–205.
[5]
D. Kirschner, J. Panetta, Modelling immunotherapy of the tumor-immune interaction, J. Math. Biol., 37 (1998), 235–252. doi: 10.1007/s002850050127
[6]
L. G. de Pillis, A. Radunskaya, The dynamics of an optimally controlled tumor model: A case study, Math. Comput. Model., 37 (2003), 1221–1244.
[7]
N. Tsur, Y. Kogana, M. Rehmc, Z. Agur, Response of patients with melanoma to immune checkpoint blockade-insights gleaned from analysis of a new mathematical mechanistic model, J. Theor. Biol., 485 (2020), 110033. doi: 10.1016/j.jtbi.2019.110033
[8]
M. A. Owen, J. A. Sherratt, Mathematical modelling macrophage dynamics in tumors, Math. Model. Meth. Appl. Sci., 9 (1999), 513–539. doi: 10.1142/S0218202599000270
[9]
J. L. Yu, S. R. J. Jang, A mathematical model of tumor-immune interactions with an immune checkpoint inhibitor, Appl. Math. Comput., 362 (2019), 124532.
[10]
R. A. Ku-Carrillo, S. E. Delgadillo, B. M. Chen-Charpentier, A mathematical model for the effect of obesity on cancer growth and on the immune system response, Appl. Math. Model., 40 (2016), 4908–4920. doi: 10.1016/j.apm.2015.12.018
[11]
N. Bellomo, L. Preziosi, Modelling and mathematical problems related to tumor evolution and its interaction with the immune system, Math. Comput. Model., 32 (2000), 413–452. doi: 10.1016/S0895-7177(00)00143-6
[12]
V. A. Kuznetsov, I. A. Makalkin, M. A. Taylor, A. S. Perelson, Nonlinear dynamics of immunogenic tumors: Parameter estimation and global bifurcation analysis, Bull. Math. Biol., 56 (1994), 295–321. doi: 10.1016/S0092-8240(05)80260-5
[13]
C. Letellier, F. Denis, L. A. Aguirre, What can be learned from a chaotic cancer model?, J. Theor. Biol., 322 (2013), 7–16. doi: 10.1016/j.jtbi.2013.01.003
[14]
A. Zazoua, W. D. Wang, Analysis of mathematical model of prostate cancer with androgen deprivation therapy, Commun. Nonlinear. Sci. Numer. Simul., 66 (2019), 41–60. doi: 10.1016/j.cnsns.2018.06.004
[15]
L. Y. Pang, Z. Zhao, X. Y. Song, Cost-effectiveness analysis of optimal strategy for tumor treatment, Chaos Solitons Fractals, 87 (2016), 293–301. doi: 10.1016/j.chaos.2016.03.032
[16]
S. Khajanchi, S. Banerjee, Quantifying the role of immunotherapeutic drug T11 target structure in progression of malignant gliomas: Mathematical modeling and dynamical perspective, Math. Biosci., 289 (2017), 69–77. doi: 10.1016/j.mbs.2017.04.006
[17]
S. Banerjee, S. Khajanchi, S. Chaudhur, A mathematical model to elucidate brain tumor abrogation by immunotherapy with T11 target structure, PLoS One, 10 (2015), 1–21. doi: 10.1371/journal.pone.0116884
[18]
J. Arciero, T. Jackson, D. Kirschner, A mathematical model of tumor-immune evasion and siRNA treatment, Discret. Contin. Dyn. Syst. Ser. B., 4 (2004), 39–58.
[19]
L. G. depillis, A. Eladdadi, A. E. Radunskaya, Modeling cancer-immune responses to therapy, J. Pharmacokinet. Pharmacodyn., 41 (2014), 461–0478. doi: 10.1007/s10928-014-9386-9
[20]
M. Itik, M. U. Salamci, S. P. Banks, Optimal control of drug therapy in cancer treatment, Nonlinear Anal., 71 (2009), 1473–1486. doi: 10.1016/j.na.2009.01.214
[21]
X. D. Liu, Q. Z. Li, J. X. Pan, A deterministic and stochastic model for the system dynamics of tumor-immune responses to chemotherapy, Phys. A, 500 (2018), 162–176. doi: 10.1016/j.physa.2018.02.118
[22]
L. G. depillis, W. Gu, K. R. Fister, T. Head, K. Maples, A. Murugan, et al., Chemotherapy for tumors: An analysis of the dynamics and a study of quadratic and linear optimal controls, Math. Biosci., 209 (2007), 292–315. doi: 10.1016/j.mbs.2006.05.003
[23]
P. Rokhforoz, A. A. Jamshidi, N. N. Sarvestani, Adaptive robust control of cancer chemotherapy with extended Kalman filter observer, Inform. Med. Unlock., 8 (2017), 1–7. doi: 10.1016/j.imu.2017.03.002
[24]
L. G. de Pillis, W. Gu, A. E. Radunskaya, Mixed immunotherapy and chemotherapy of tumors: modeling, applications and biological interpretations, J. Theor. Biol., 238 (2006), 841–862. doi: 10.1016/j.jtbi.2005.06.037
[25]
R. Eftimie, J. J. Gillard, D. A. Cantrell, Mathematical models for immunology: Current state of the art and future research directions, Bull. Math. Biol., 78 (2016), 2091–2134. doi: 10.1007/s11538-016-0214-9
[26]
M. Villasana, A. Radunskaya, A delay differential equation model for tumor growth, J. Math. Biol., 47 (2003), 270–294. doi: 10.1007/s00285-003-0211-0
[27]
Y. Radouane, Hopf bifurcation in differential equations with delay for tumor-immune system competition model, SIAM J. Appl. Math., 67 (2007), 1693–1703. doi: 10.1137/060657947
[28]
F. A. Rihan, D. H. Abdel Rahman, S. Lakshmanan, A. S. Alkhajeh, A time delay model of tumour-immune system interactions: Global dynamics, parameter estimation, sensitivity analysis, Appl. Math. Comput., 232 (2014), 606–623.
[29]
M. Yu, Y. P. Dong, Y. Takeuchi, Dual role of delay effects in a tumour-immune system, J. Biol. Dyn., 11 (2017), 334–347. doi: 10.1080/17513758.2016.1231347
[30]
S. Khajanchi, S. Banerjee, Stability and bifurcation analysis of delay induced tumor immune interaction model, Appl. Math. Comput., 248 (2014), 652–671.
[31]
S. Khajanchi, Bifurcation analysis of a delayed mathematical model for tumor growth, Chaos Solitons Fractals, 77 (2015), 264–276. doi: 10.1016/j.chaos.2015.06.001
[32]
S. Khajanchi, S. Banerjee, Influence of multiple delays in brain tumor and immune system interaction with T11 target structure as a potent stimulator, Math. Biosci., 302 (2018), 116–130. doi: 10.1016/j.mbs.2018.06.001
[33]
M. J. Piotrowska, M. Bodnar, Influence of distributed delays on the dynamics of a generalized immune system cancerous cells interactions model, Commun. Nonlinear. Sci. Numer. Simul., 54 (2018), 389–415. doi: 10.1016/j.cnsns.2017.06.003
[34]
S. Banerjee, R. R. Sarkar, Delay-induced model for tumor-immune interaction and control of malignant tumor growth, BioSystems, 91 (2008), 268–288. doi: 10.1016/j.biosystems.2007.10.002
[35]
A. E. Gohary, Chaos and optimal control of cancer self-remission and tumor system steady states, Chaos Solitons Fractals, 37 (2008), 1305–1316. doi: 10.1016/j.chaos.2006.10.060
[36]
H. Mayer, K. Zaenker, U. an der Heiden, A basic mathematical model of the immune response, Chaos, 5 (1995), 155–161. doi: 10.1063/1.166098
[37]
H. M. Byrne, The effect of time delay on the dynamics of avascular tumor growth, Math. Biosci., 144 (1997), 83–117. doi: 10.1016/S0025-5564(97)00023-0
[38]
P. Bi, S. G. Ruan, Bifurcations in delay differential equations and applications to tumor and immune system interaction models, SIAM J. Appl. Dyn., 12 (2014), 1847–1888.
[39]
S. G. Ruan, Nonlinear dynamics in tumor-immune system interaction models with delays, Discret. Contin. Dyn. Syst. Ser. B., 26 (2021), 541–602.
[40]
F. A. Rihan, Sensitivity analysis of dynamic systems with time lags, J. Comput. Appl. Math., 151 (2003), 445–462. doi: 10.1016/S0377-0427(02)00659-3
[41]
K. Gopalsamy, Nonoscillation in a delay-logistic equation, Quart. Appl. Math., 43 (1985), 189–197. doi: 10.1090/qam/793526
[42]
D. Mukherjee, P. C. Bhakta, A. B. Roy, Uniform persistence in Kolmogorov models with convex growth functions, Nonlinear Anal., 34 (1998), 427–432. doi: 10.1016/S0362-546X(97)00605-6
[43]
J. J. Wei, S. G. Ruan, Stability and bifurcation in a neural network model with two delays, Phys. D, 130 (1998), 255–272.
[44]
B. D. Hassard, N. D. Kazarinoff, Y. H. Wan, Theory and Application of Hopf Bifurcation, London Mathematical Society Lecture Note Series, Cambridge University Press, 1981.
This article has been cited by:
1.
Larysa Dzyubak, Oleksandr Dzyubak, Jan Awrejcewicz,
Nonlinear multiscale diffusion cancer invasion model with memory of states,
2023,
168,
09600779,
113091,
10.1016/j.chaos.2022.113091
2.
Yuliana Jao, Nur Erawaty,
Dynamic study of the pathogen-immune system interaction with natural delaying effects and protein therapy,
2022,
7,
2473-6988,
7471,
10.3934/math.2022419
3.
Ting Yu, Qinglong Wang, Shuqi Zhai,
Exploration on dynamics in a ratio-dependent predator-prey bioeconomic model with time delay and additional food supply,
2023,
20,
1551-0018,
15094,
10.3934/mbe.2023676
4.
Jianping Li, Nan Liu, Danni Wang, Hongli Yang,
Mathematical modeling and Hopf bifurcation analysis of tumor macrophage interaction with polarization delay,
2025,
19,
1751-3758,
10.1080/17513758.2025.2508240
Qingfeng Tang, Guohong Zhang. Stability and Hopf bifurcations in a competitive tumour-immune system with intrinsic recruitment delay and chemotherapy[J]. Mathematical Biosciences and Engineering, 2021, 18(3): 1941-1965. doi: 10.3934/mbe.2021101
Qingfeng Tang, Guohong Zhang. Stability and Hopf bifurcations in a competitive tumour-immune system with intrinsic recruitment delay and chemotherapy[J]. Mathematical Biosciences and Engineering, 2021, 18(3): 1941-1965. doi: 10.3934/mbe.2021101
Note: GAS means globally asymptotically stable; LAS means locally asymptotically stable.
Case
Conditions
Equilibria and stability
m1>μ1
m3>μ2
Exyz0 is GAS
m3<μ2, μ2m2+ρm3≥μ2(ρ−k)
Exy0 is GAS
μ2m2+ρm3<μ2(ρ−k)
Ex0 is LAS
m1<μ1
m3>μ2
Eyz0 is GAS
m3<μ2, μ2m2+ρm3>μ2(ρ−k)+μ2α(m1−μ1)μ1
Ey0 is LAS
μ2m2+ρm3<μ2(ρ−k)+ργ(m1−μ1)
Ex0 is LAS
Note: GAS means globally asymptotically stable; LAS means locally asymptotically stable.
Figure 1. (a) Basins of attraction of Ex0 and Ey0. (b) (c)Trajectories are shown to converge to either Ex0 or Ey0 depending on their initial positions. The initial values in 1(b) is [x(0),y(0),z(0)]=[0.9,0.6,0.2] and in 1(c) is [x(0),y(0),z(0)]=[0.9,0.5,0.2]. 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
Figure 2. The coexistence equilibrium E∗=(0.6,0.07,0.8) is locally asymptotically stable when time delay is absent. Normalized parameters are given in Eq (5.2)
Figure 3. Time series and phase diagrams of the solutions of Eq (1.4) with different time delays τ. (a) (d):the positive equilibrium E∗=(0.6,0.07,0.8) of Eq (1.4) is local asymptotically stable when τ=14<τ0=14.5117. (b) (c) (e) (f): the bifurcated periodic solutions occur through Hopf bifurcations when τ=15,20,25>τ0. Normalized parameters are given in Equation (5.2)
Figure 4. Amplitude and period of tumor cells population for the bifurcated periodic solutions of Eq (1.4) with respect to time delay τ. One can note that larger intrinsic recruitment delay lead to larger amplitude and longer period of the bifurcated periodic solution. Normalized parameters are given in Eq (5.2)
Figure 5. Sensitivity functions ∂x∂m1, ∂x∂m2, ∂x∂m3 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 m1 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