
Citation: Hongying Shu, Wanxiao Xu, Zenghui Hao. Global dynamics of an immunosuppressive infection model with stage structure[J]. Mathematical Biosciences and Engineering, 2020, 17(3): 2082-2102. doi: 10.3934/mbe.2020111
[1] | Jinhu Xu, Yicang Zhou . Bifurcation analysis of HIV-1 infection model with cell-to-cell transmission and immune response delay. Mathematical Biosciences and Engineering, 2016, 13(2): 343-367. doi: 10.3934/mbe.2015006 |
[2] | 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 |
[3] | Shunyi Li . Hopf bifurcation, stability switches and chaos in a prey-predator system with three stage structure and two time delays. Mathematical Biosciences and Engineering, 2019, 16(6): 6934-6961. doi: 10.3934/mbe.2019348 |
[4] | 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 |
[5] | 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 |
[6] | Kaushik Dehingia, Anusmita Das, Evren Hincal, Kamyar Hosseini, Sayed M. El Din . Within-host delay differential model for SARS-CoV-2 kinetics with saturated antiviral responses. Mathematical Biosciences and Engineering, 2023, 20(11): 20025-20049. doi: 10.3934/mbe.2023887 |
[7] | Yan Wang, Minmin Lu, Daqing Jiang . Viral dynamics of a latent HIV infection model with Beddington-DeAngelis incidence function, B-cell immune response and multiple delays. Mathematical Biosciences and Engineering, 2021, 18(1): 274-299. doi: 10.3934/mbe.2021014 |
[8] | Wanxiao Xu, Ping Jiang, Hongying Shu, Shanshan Tong . Modeling the fear effect in the predator-prey dynamics with an age structure in the predators. Mathematical Biosciences and Engineering, 2023, 20(7): 12625-12648. doi: 10.3934/mbe.2023562 |
[9] | 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 |
[10] | Juan Wang, Chunyang Qin, Yuming Chen, Xia Wang . Hopf bifurcation in a CTL-inclusive HIV-1 infection model with two time delays. Mathematical Biosciences and Engineering, 2019, 16(4): 2587-2612. doi: 10.3934/mbe.2019130 |
According to clinical data, drug treatment against persistent human infections sometimes fails to consistently eradicate the infection from the host [1,2,3], such as human immunodeficiency virus (HIV), hepatitis B virus (HBV), and hepatitis C virus (HCV). For all these viruses, drug treatment is not effective, and for HIV, lifelong therapy is commonly required to control viral replication [4,5]. The main reason is that these viruses are able to suppress and impair immune responses, resulting in persistent infection [6,7]. An alternative strategy is to use drug therapies to boost virus-specific immune response and then induce sustained viral suppression in the absence of drugs. Antiviral therapy aimed at boosting specific immune responses has attracted more and more attentions recently. In 2003, Komarova et al. [6] used mathematical models to study immune response dynamics during therapy in the context of immunosuppressive infections. They showed that a single phase of antiviral drug therapy is also possible to establish sustained immunity. They assumed that the degree of immune expansion depends on virus load and the response inhibits virus growth. This assumption is natural for any branch of the adaptive immune system, such as CD4 T cells, CD8 T cells, or antibodies. Denote by y(t) and z(t), respectively, the virus population and immune cell population at time t. Komarova et al. [6] proposed the following system of ordinary differential equations:
y′(t)=ry(t)(1−y(t)K)−py(t)z(t)−by(t),z′(t)=cy(t)z(t)1+dy(t)−qy(t)z(t)−μz(t). | (1.1) |
The virus population grows at a rate described by the logistic function. The viral replication rate is a linear decreasing function of viral loads with a maximum value of r and it vanishes at a viral load K. The treatment is modeled as a reduction in r. Immune cells kill virus at a rate pyz. b and μ are the natural decay rates of viral particles and immune cells, respectively. Immune expansion is modeled by the growth function cyz1+dy. When the virus load is low, the level of immune response is simply proportional to both the viral load and the immune cells. The immune response saturates when the virus load is sufficiently high. Immune cells are inhibited by the virus at a rate qyz. Komarova et al. presented a simple relationship between the timing of therapy and efficacy of the drugs. In the presence of strong viral suppression, they showed that therapy should be stopped relatively early because a longer duration of treatment may lead to a failure. On the other hand, in the presence of weaker viral suppression, stopping treatment too early is detrimental, and therapy has to be continued beyond a time threshold. Essentially, model (1.1) has two stable equilibria: a virus dominant equilibrium with no sustained immunity and an immune control equilibrium with sustained immunity. This bistability allows a solution from the basin of the attraction of the virus dominant equilibrium to be lifted to that of the immune control equilibrium via a single phase of therapy.
Model (1.1) assumes that the immune response is instantaneous. Nevertheless, it has been established that the immune response process involves a sequence of events such as antigenic activation, selection, and proliferation of the immune cells [8]. Moreover, oscillatory viral loads and immune cells were observed from clinical data [9], while the oscillatory behavior is not exhibited in model (1.1). A natural question is what is the cause of the oscillations. In 2014, Shu et al. [10] incorporated the time lag during the immune response process into model (1.1). By studying the dynamics between an immunosuppressive infection and antiviral immune response, they demonstrated that the time lag is the main factor causing oscillations. Here, τ is denoted as the time lag for the immune system to trigger a sequence of events. The activation rate of immune cells at time t depends on the virus load and the number of immune cells at time t−τ. The above assumptions lead to the following system:
y′(t)=ry(t)(1−y(t)K)−py(t)z(t)−by(t),z′(t)=cy(t−τ)z(t−τ)1+dy(t−τ)−qy(t)z(t)−μz(t). | (1.2) |
It was shown that the delayed antiviral immune response may induce sustained periodic oscillations, transient oscillations and even sustained aperiodic oscillations (chaos), and the model can admit two stable equilibria and can also allow a stable equilibrium to coexist with a stable periodic solution.
In the aforementioned work, the mortality of the population during the time lag has been ignored. Consideration of the survival probability during the time lag requires an additional delay-dependent multiplier in the nonlinear delayed feedback term. The death of the free viral particles depends on the mature immune cells instead of all the newly generated immune cells. We denote z(t) as the population of mature immune cells at time t, and consider the following delay differential equation with a delay-dependent coefficient:
y′(t)=f(y(t))−py(t)z(t)−by(t),z′(t)=e−δτh(y(t−τ))z(t−τ)−qy(t)z(t)−μz(t). | (1.3) |
Here, f(y) describes the virus population growth function; δ>0 is the death rate of the immature immune cells during the time lag; e−δτ represents the survival probability of the immature immune cells surviving τ time units before becoming mature; h(y)z denotes the immune cell growth function; the other parameters have the same meaning as in model (1.2). In addition to more general virus and immune cell growth functions, model (1.3) differs from model (1.2) in the sense that it introduces an additional delay-dependent survival probability e−δτ. As we shall see later, this quantity will make a significant difference in the bifurcation analysis. For instance, model (1.3) does not possess any positive equilibrium for sufficiently large τ and its global Hopf branches always have bounded τ components. However, for model (1.2), the existence condition of positive equilibrium is independent of τ and the global Hopf branches always have unbounded τ components. Furthermore, there exists only one stability switch of positive equilibrium for model (1.2), but the stability of positive equilibrium for our model (1.3) may switch multiple times.
Actually, the model (1.3) can be derived from a stage structured population model for u(t,a) as below
y′(t)=f(y(t))−py(t)z(t)−by(t),∂tu(t,a)+∂au(t,a)=−(μ(a)+q(a)y(t))u(t,a), |
with the stage-specific decay rate of immune cells, μ(a), inhibition rate of immune cells, q(a),
μ(a)={μ,a≥τ,δ,a<τ, q(a)={q,a≥τ,0,a<τ, |
where u(t,a) is the population of immune cells at age a and time t. Note that the population of mature immune cells is z(t)=∫∞τu(t,a)da. We integrate along characteristic lines to find
z′(t)=u(t,τ)−u(t,∞)−qy(t)z(t)−μz(t)=u(t−τ,0)e−δτ−u(t,∞)−qy(t)z(t)−μz(t). |
Note that u(t,∞)=0 and u(t−τ,0)=h(y(t−τ))z(t−τ). The equation for z(t) follows.
Throughout this paper, we assume that the virus population growth function f(y) and the immune expansion rate h(y) satisfy the following conditions:
(H1) f(y) is continuously differentiable, f(0)=0 and f′(0)>b; there exists ˉy>0 such that f(ˉy)=0, and (y−ˉy)f(y)<0 for positive y≠ˉy; f″(y)<0 for y∈[0,ˉy].
(H2) h(y) is continuously differentiable, h(0)=0 and h′(0)>q; h′(y)>0 and h″(y)<0 for all y≥0.
The assumptions (H1)-(H2) are biological conditions for infection dynamics. For instance, f(0)=0 means no new infection in the absence of virus; f′(0)>b indicates that the intrinsic growth rate is greater than the decay rate. Similarly, h(0)=0 implies that immune cells cannot reproduce without the presence of virus; h′(0)>q guarantees that the initial expansion rate is greater than the inhibition rate.
The rest of this paper is organized as follows. We first present some preliminary results concerning the well-posedness of model (1.3), and describe global dynamics of the ordinary differential equation system (2.4) in Section 2. Local and global stability analysis of equilibria as well as local and global Hopf bifurcation analysis are given in Section 3. In Section 4, numerical simulations based on the bifurcation analysis are reported and discussed. Finally, we give a brief summary in Section 5.
To establish the well-posedness of model (1.3), we choose the phase space C×C, where C is the Banach space of continuous functions on [−τ,0] defined by C:={ϕ:[−τ,0]→R is continuous}, and the norm is defined by ‖ϕ‖=sup−τ≤θ≤0ϕ(θ). The nonnegative cone of C is denoted by C+. As usual, ϕt∈C is defined by ϕt(θ)=ϕ(t+θ) for θ∈[−τ,0]. For biological applications, the initial condition of (1.3) is given as
(y0,z0)∈X:=C+×C+. |
The existence and uniqueness of the solution of model (1.3) follows from the standard theory of functional differential equations [11]. Using the same manner as in proving [12,Proposition 2.1], we can easily show that the solution of (1.3) with initial conditions in X is nonnegative, which implies that X is positively invariant under the solution map of (1.3).
Next, we shall prove that all solutions of (1.3) are ultimately bounded. It follows from the first equation in (1.3) that y′(t)≤f(y), which yields lim supt→∞y(t)≤ˉy. By the nonnegativity of the solution of (1.3) and (H1)-(H2), we obtain
(y(t)+ph′(0)z(t+τ))′≤Mf+(e−δτh(y)h′(0)y−1)py(t)z(t)−by(t)−μph′(0)z(t+τ)≤Mf−γ(y(t)+ph′(0)z(t+τ)), |
where Mf=maxy≥0f(y) and γ=min{b,μ}. Thus,
lim supt→∞(y(t)+ph′(0)z(t+τ))≤Mfγ, |
which implies that y(t) and z(t) are ultimately bounded in X. To summarize, we obtain the following proposition.
Proposition 2.1. The solutions of model (1.3) with the initial conditions in X are nonnegative, and the region
Γ={(ϕ,ψ)∈X:‖ϕ‖≤ˉy, ϕ(−τ)+ph′(0)ψ(0)≤Mfγ} |
is positively invariant and absorbing in X; namely, all solutions ultimately enter Γ.
Clearly, model (1.3) always has a trivial equilibrium E0=(0,0). Since f′(0)>b, there also exists a boundary equilibrium Eb=(˜y,0), where ˜y is the unique positive root of f(y)=by. We call Eb the virus dominant equilibrium (VDE). Assume that E∗=(y∗,z∗) is a positive equilibrium with y∗>0 and z∗>0, then it must satisfy the following equilibrium equations
e−δτh(y∗)=qy∗+μ, z∗=f(y∗)−by∗py∗. | (2.1) |
Obviously, z∗>0 if and only if y∗<˜y. Hence, E∗ exists if and only if y∗ is a positive root of
h(y)y+μ/q=qeδτ | (2.2) |
with y∗<˜y. It is easily to calculate that sign(h(y)y+μ/q)′=sign g(y), where g(y)=h′(y)(y+μq)−h(y). Note that g(0)>0,g(∞)<0 and g′(y)=h″(y)(y+μ/q)<0. There exists a unique yc>0 such that g(yc)=0, and (y−yc)g(y)<0 for y≠yc. Therefore, h(y)y+μ/q is strictly increasing on [0,yc), and strictly decreasing on (yc,∞). We denote
s:=h(yc)yc+μ/q=supy≥0h(y)y+μ/q. | (2.3) |
If s<qeδτ, (2.2) has no positive root. If s≥qeδτ, (2.2) has exactly two positive roots (counting multiplicity) y∗1 and y∗2 such that y∗1≤yc≤y∗2. These two roots coincide (i.e., y∗1=y∗2=yc) if and only if s=qeδτ. Throughout this paper, we make the following assumption to ensure the existence of positive roots for (2.2).
(H3) s>q and τ≤τmax, where τmax=1δlnsq.
If (H3) is violated, there does not exist any positive equilibrium for our model and the boundary equilibrium is locally asymptotically stable. Summarizing the above analysis, we obtain the following result.
Proposition 2.2. Assume that (H1)-(H3) are satisfied.
(i) If y∗1≥˜y holds, then there are two equilibria: E0=(0,0) and Eb=(˜y,0).
(ii) If y∗1<˜y≤y∗2 or y∗1=y∗2<˜y holds, then there are three equilibria: E0, Eb and E∗1=(y∗1,z∗1).
(iii) If y∗1<y∗2<˜y holds, then there are four equilibria: E0, Eb, E∗1 and E∗2=(y∗2,z∗2).
When τ=0, model (1.3) reduces to an ODE system which generalizes (1.1)
y′(t)=f(y(t))−py(t)z(t)−by(t),z′(t)=h(y(t))z(t)−qy(t)z(t)−μz(t). | (2.4) |
We next provide a complete description about its global dynamics.
Lemma 2.3. System (2.4) has no closed orbits.
Proof. Using an argument similar to Proposition 2.1, we can show that
Γ0={(ϕ,ψ)∈R2+:ϕ≤ˉy, ϕ+ph′(0)ψ≤Mfγ} |
is positively invariant and absorbing in R2+. Thus, it suffices to show that System (2.4) has no closed orbits in Γ0.
Let (P(y,z),Q(y,z)) denote the vector field of (2.4), then we have
∂(BP)∂y+∂(BQ)∂z=1z(t)(f(y)y)′, where B(y,z)=1yz. |
Note that (H1) implied that (f(y)y)′<0 for all y≤ˉy. Thus, ∂(BP)∂y+∂(BQ)∂z<0 in Γ0. The nonexistence of closed orbits for (2.4) follows from the classical Dulac-Bendixson criterion. This ends the proof.
We linearize system (2.4) at E0=(0,0) and calculate that the two eigenvalues are f′(0)−b>0 and −μ<0. Thus, E0 is a saddle point. For the linearized system at Eb=(˜y,0), the two eigenvalues are
λ1=f′(˜y)−b<0 and λ2=h(˜y)−q˜y−μ. |
From the proof of Proposition 2.2, we know h(˜y)−q˜y−μ<0 if ˜y<y∗1 or ˜y>y∗2, and h(˜y)−q˜y−μ>0 provided that y∗1<˜y<y∗2. This shows that Eb is stable if y∗1>˜y or y∗2<˜y, and Eb is a saddle point if y∗1<˜y<y∗2. For the linearized system at E∗1=(y∗1,z∗1), the associated characteristic equation is
λ2+(−f′(y∗1)+pz∗1+b)λ+py∗1z∗1(h′(y∗1)−q)=0. |
It follows from (2.1) and (H1) that −f′(y∗1)+pz∗1+b=−f′(y∗1)+f(y∗1)/y∗1>0. This, together with h′(y∗1)−q>0, implies that all eigenvalues must have negative real parts, and hence E∗1 is asymptotically stable. At E∗2=(y∗2,z∗2), the two eigenvalues λ1 and λ2 satisfy
λ1+λ2=f′(y∗2)−f(y∗2)y∗2<0, λ1λ2=py∗2z∗2(h′(y∗2)−q)<0. |
This implies that one eigenvalue must be positive and the other negative. Thus, E∗2 is a saddle point whenever it exists. Summarizing the above analysis, we obtain the following global stability results on model (2.4).
Theorem 2.4. Assume that (H1)-(H2) are satisfied and s>q.
(i) If y∗1>˜y holds, then E0 is a saddle point and Eb attracts all solutions with initial conditions in {(y(t),z(t))∈R2+:y(0)>0,z(0)≥0}.
(ii) If y∗1<˜y<y∗2 holds, then both E0 and Eb are saddle points and E∗1 attracts all solutions with initial conditions in {(y(t),z(t))∈R2+:y(0)>0,z(0)>0}.
(iii) If y∗1<y∗2<˜y holds, then both E0 and E∗2 are saddle points, and for any initial condition (y(0),z(0)) with y(0)>0, the solution of (2.4) approaches either Eb or E∗1, the stable manifold of E∗2 separate the basins of attraction of Eb and E∗1.
In this section, we study the dynamics of model (1.3) with τ>0.
By linearizing (1.3) at E0=(0,0), we find two eigenvalues f′(0)−b>0 and −μ<0. Hence, E0 is a saddle point. The characteristic equation associated with the linearization of model (1.3) at Eb is
(λ−f′(˜y)+b)(λ+q˜y+μ−e−δτh(˜y)e−λτ)=0. |
One eigenvalue is λ1=f′(˜y)−b<0. Hence Eb is locally asymptotically stable if all zeros of
Δb(λ)=λ+q˜y+μ−e−δτh(˜y)e−λτ |
have negative real parts. By [13,Lemma 6], there exists at least one positive eigenvalue if q˜y+μ<e−δτh(˜y), or equivalently, y∗1<˜y<y∗2. On the other hand, if q˜y+μ>e−δτh(˜y), that is, either ˜y<y∗1 or ˜y>y∗2, then all eigenvalues have negative real parts. For the critical case q˜y+μ=e−δτh(˜y); namely, either ˜y=y∗1 or ˜y=y∗2, one eigenvalue is 0, and all other eigenvalues have negative real parts. By using the method of normal forms, we obtain that Eb is locally asymptotically stable. To summarize, we have the following result on the local stability of E0 and Eb.
Theorem 3.1. Consider model (1.3) under the assumptions (H1)-(H3). E0 is a saddle point. If either ˜y≤y∗1 or ˜y≥y∗2, then Eb is locally asymptotically stable; while Eb is unstable if y∗1<˜y<y∗2.
To establish the global stability of Eb, we first show that the infection is persistent.
Lemma 3.2. lim inft→∞y(t)>0 for any solution of (1.3) with the initial condition in X1={(ϕ,ψ)∈X:ϕ(0)>0,ψ(0)≥0}.
Proof. We claim that lim supt→∞y(t)>0. If not, then y(t)→0 as t→∞. It then follows from (1.3) that z(t)→0 as t→∞. This contradicts to the fact that E0 is unstable. Thus, y(t) is weakly persistent. By using [14,Theorem 3.4.6] and Proposition 2.1, there exists a global attractor for the solution semiflow of (1.3). This, together with [15,Theorem 2.3], implies that y(t) is actually strongly persistent; that is, lim inft→∞y(t)>0.
Our next result is concerned with global stability of Eb.
Theorem 3.3. Assume that (H1)-(H3) hold. Then Eb of model (1.3) is globally asymptotically stable in the region X1 provided that τ≥τb, where τb:=1δlnh′(0)˜yq˜y+μ. Moreover, the condition τ≥τb implies that ˜y≤y∗1.
Proof. By Proposition 2.1 and Theorem 3.1, it suffices to show that Eb is globally attractive in X1∩Γ, which is a positively invariant set of (1.3). Let (y(t),z(t)) be a solution of (1.3) with the initial condition in X1, it follows from Lemma 3.2 that y(t)>0 for t>0. Motivated by the earlier work in [16], we construct a Lyapunov functional L:X1∩Γ→R by
L(yt,zt)=yt(0)−˜yln(yt(0))+czt(0)+ce−δτ∫0−τh(yt(θ))zt(θ)dθ, |
where c is a constant to be determined later. Calculating the time derivative of L along the solution, and using b=f(˜y)/˜y, h(y)≤h′(0)y for all y≥0, we obtain
dLdt=f(y)y(y−˜y)−b(y−˜y)+(p˜y−cμ)z+ce−δτh(y)z−(cq+p)yz,≤(f(y)y−f(˜y)˜y)(y−˜y)+(p˜y−cμ)z+(ce−δτh′(0)−cq−p)yz. |
Assumption (H1) implies that (f(y)y−f(˜y)˜y)(y−˜y)≤0 for all y≤ˉy. Since h′(0)e−δτ−q>0, the condition τ≥τb is the same as p˜yμ≤pe−δτh′(0)−q. We may choose a constant c∈[p˜yμ,pe−δτh′(0)−q] such that dL/dt≤0. Moreover, dL/dt=0 if and only if y(t)=˜y and z(t)=0. Thus the maximal compact invariant set in {dL/dt=0} is the singleton {Eb}. By the LaSalle invariance principle [11], Eb is globally attractive in X1∩Γ. Since X1∩Γ is absorbing in X1, we conclude that Eb is globally attractive in X1. By Theorem 3.2, Eb is globally asymptotically stable in X1 if τ≥τb.
If τ≥τb, then e−δτh′(0)˜y−q˜y−μ≤0=e−δτh(y∗1)−qy∗1−μ≤e−δτh′(0)y∗1−qy∗1−μ. Since e−δτh′(0)>q, we obtain ˜y≤y∗1.
Here, we remark that if y∗1=y∗2<˜y, then two positive equilibria E∗1 and E∗2 coincide, and Eb is locally asymptotically stable.
In this subsection, we investigate the stability of E∗1 and identify parameter regions in which the time delay can destabilize E∗1, lead to Hopf bifurcation and induce sustained oscillations.
The characteristic equation associated with the linearization of system (1.3) at equilibrium E∗i (i=1,2) is
Δi(λ)=λ2+a1,iλ+a0,i+(b1,iλ+b0,i)e−λτ=0, | (3.1) |
where
a0,i=(f(y∗i)y∗i−f′(y∗i))h(y∗i)e−δτ−pqy∗iz∗i, a1,i=f(y∗i)y∗i−f′(y∗i)+h(y∗i)e−δτ>0,b0,i=(f′(y∗i)−f(y∗i)y∗i)h(y∗i)e−δτ+py∗iz∗ih′(y∗i)e−δτ, b1,i=−h(y∗i)e−δτ<0. |
We first investigate the stability of E∗1, regarding τ as the bifurcation parameter. In the proof of Theorem 2.4, we have demonstrated that, when τ=0, all eigenvalues of (3.1) with i=1 lie to the left of the imaginary axis and E∗1 is locally asymptotically stable. As τ increases, E∗1 may lose its stability only when some eigenvalues cross the imaginary axis to the right. In view of
a0,1+b0,1=py∗1z∗1(h′(y∗1)e−δτ−q)>0, |
0 is not an eigenvalue for any τ≥0. For simplicity, we drop the subscript 1 in the following arguments. Substituting λ=iω(ω>0) into (3.1) and separating the real and imaginary parts, we obtain
ω2−a0=b0cos(ωτ)+b1ωsin(ωτ),−a1ω=b1ωcos(ωτ)−b0sin(ωτ). |
Squaring and adding both the above equations lead to
G(ω,τ):=ω4+c1(τ)ω2+c0(τ)=0, | (3.2) |
where
c1(τ)=(f(y∗1)y∗1−f′(y∗1))2+2pqy∗1z∗1>0, c0(τ)=a20−b20. |
Then iω(ω>0) is an imaginary root of (3.1) if and only if G(ω,τ)=0. Since c1(τ)>0 for all 0≤τ≤τmax. we know that G(ω,τ)=0 has exactly one positive root if and only if c0(τ)<0, or equivalently, a0<b0.
If a0<b0, the implicit function theorem implies that there exists a unique C1 function ω=ω∗(τ) such that G(ω∗(τ),τ)=0 for 0<τ≤τmax. For iω∗(τ) to be a root of (3.1), ω∗(τ) needs to satisfy the system
sin(ω(τ)τ)=b1ω(τ)3+(a1b0−a0b1)ω(τ)b21ω(τ)2+b20:=g1(τ),cos(ω(τ)τ)=(b0−a1b1)ω(τ)2−a0b0b21ω(τ)2+b20:=g2(τ). | (3.3) |
Set
I={τ:0≤τ≤τmax satisfies a0<b0}. | (3.4) |
For τ∈I, let θ(τ) be the unique solution of sinθ=g1 and cosθ=g2 in (0,2π]; that is,
θ(τ)={arccosg2(τ), if ω(τ)2≤a0−a1b0b1,2π−arccosg2(τ), if ω(τ)2>a0−a1b0b1. |
Following Beretta and Kuang [17], we define
Sn(τ)=τ−θ(τ)+2nπω(τ) for τ∈I with n∈N. | (3.5) |
One can check that ±iω∗(τ∗) are a pair of purely imaginary eigenvalues of (3.1) if and only if τ∗ is a zero of function Sn(τ) for some n∈N. From [17,Theorem 2.2], we have
Sign(dReλ(τ)dτ|τ=τ∗)=Sign(∂G∂ω(ω∗(τ∗),τ∗))SignS′n(τ∗). |
Note that ∂G∂ω(ω∗(τ∗),τ∗)>0, thus we have the following result concerning the transversality condition.
Lemma 3.4. Assume that Sn(τ) has a positive root τ∗∈I for some n∈N, then a pair of simple purely imaginary roots ±iω∗(τ∗) of (3.1) exist at τ=τ∗, and
Sign(dReλ(τ)dτ|τ=τ∗)=SignS′n(τ∗). |
Moreover, this pair of simple purely imaginary roots ±iω∗(τ∗) cross the imaginary axis from left to right at τ=τ∗ if S′n(τ∗)>0, and from right to left if S′n(τ∗)<0.
It is easily seen that Sn(τ)>Sn+1(τ) for all τ∈I,n∈N. When τ=0, the asymptotic stability of E∗1 implies that S0(0)<0. Denote
ˆτ=supI:=sup{τ:0≤τ≤τmax satisfies a0<b0}. | (3.6) |
If ˆτ<τmax, then c0(ˆτ)=0, and thus ω(τ)→0 as τ→ˆτ−. This, together with (3.3), yields limτ→ˆτ−sinθ(τ)=0 and limτ→ˆτ−cosθ(τ)=−1. Therefore, limτ→ˆτ−θ(τ)=π, limτ→ˆτ−Sn(τ)=−∞.
If supτ∈IS0(τ)<0, Sn(τ) has no zeros in I for all n∈N, which excludes the existence of purely imaginary eigenvalues and thus implies that E∗1 is locally asymptotically stable for all τ∈[0,τmax].
If supτ∈IS0(τ)=0, S0(τ) has a double zero in I, denoted by τ∗, and S′0(τ∗)=0. This, together with Lemma 3.4, implies that the transversality condition at τ∗ is not satisfied and all eigenvalues remain to the left of the pure imaginary axis. Thus, E∗1 is still locally asymptotically stable for τ∈[0,τmax].
If supτ∈IS0(τ)>0, it then follows from S0(0)<0 and limτ→ˆτ−Sn(τ)=−∞ that S0(τ) has at least two zeros in I. For simplicity, we assume that
(H4) supτ∈IS0(τ)>0 and all zeros of Sn(τ) with odd multiplicity are simple.
Under assumption (H4), each Sn(τ) has even number of simple zeros. Now, we collect all simple zeros of Sn(τ) with n≥0 and list them in increasing order: 0<τ0<τ1<⋯<τ2K−1<ˆτ (K∈N+). For each integer 0≤i≤K−1, we apply Lemma 3.4 to obtain S′m(τi)>0 and S′m(τ2K−i−1)<0 for some m∈N. Therefore, the pair of simple conjugate purely imaginary eigenvalues ±iω(τi) crosses the imaginary axis from left to right, and the pair of simple conjugate purely imaginary eigenvalues ±iω(τ2K−i−1) crosses the imaginary axis from right to left. Thus, system (1.3) undergoes a Hopf bifurcation at E∗1 when τ=τj (0≤j≤2K−1). Moreover, E∗1 is asymptotically stable for τ∈[0,τ0)∪(τ2K−1,τmax], and unstable for τ∈(τ0,τ2K−1).
For each k=0,⋯,2K−1, there exists nk such that τk is a simple zero of Snk(τ). Let Tk be the period of periodic solution bifurcated at τk. Applying the Hopf bifurcation theorem for delay differential equations [11,18], we have
Tk=2πω(τk)=2πτkθ(τk)+2nkπ. |
This, together with θ(τk)∈(0,2π], implies that Tk>τk if nk=0, and τknk+1≤Tk<τknk if nk>0. To summarize, we have the following results on stability of E∗1 and Hopf bifurcation.
Theorem 3.5. Assume that (H1)-(H3) hold, y∗1<˜y and ˆτ<τmax. Let I,Sn(τ),τmax,ˆτ be defined in (3.4), (3.5), (H3) and (3.6), respectively.
(i) If either I=∅ or supτ∈IS0(τ)≤0, then E∗1 is locally asymptotically stable for all τ∈[0,τmax].
(ii) If (H4) holds, then there exist exactly 2K local Hopf bifurcation values, namely, 0<τ0<τ1<⋯<τ2K−1<ˆτ such that model (3.1) undergoes a Hopf bifurcation at E∗1 when τ=τj for 0≤j≤2K−1. E∗1 is locally asymptotically stable for τ∈[0,τ0)∪(τ2K−1,τmax], and unstable for τ∈(τ0,τ2K−1). Moreover, the period Tk of periodic solution bifurcated at τk satisfies Tk>τk if τk is a simple zero of S0(τ), and τknk+1≤Tk<τknk if τk is a simple zero of Snk(τ) with nk>0.
Theorem 3.6. Assume that (H1)-(H3) are satisfied. If y∗1<y∗2<˜y holds, then E∗2 is unstable for all τ≥0. Moreover, if a0,2≤b0,2, then the characteristic equation (3.1) with i=2 has no purely imaginary eigenvalues.
Proof. Note that in (3.1) with i=2, Δ2(0)=py∗2z∗2(h′(y∗2)e−δτ−q)<0. For any τ≥0, we know that Δ2(∞)>0, thus as long as E∗2 exists, the associated characteristic equation must have at least one positive real root and hence E∗2 is always unstable for τ≥0.
Using a similar argument in the study of root distribution for (3.1) at E∗1, we note that c1,2(τ)>0,c0,2(τ)=a20,2−b20,2 and a0,2+b0,2=py∗2z∗2(h′(y∗2)e−δτ−q)<0. Thus, G(ω,τ) in (3.2) has no positive zeros if a0,2≤b0,2. Therefore, (3.1) with i=2 has no purely imaginary eigenvalues if a0,2≤b0,2.
Since E∗2 is unstable (and thus biologically irrelevant) whenever it exists, we name E∗1, instead of E∗2, the immune control equilibrium (ICE).
Note that Theorem 3.5 states that if y∗1<˜y and (H4) holds, then periodic solutions can bifurcate from E∗1 when τ is near the local Hopf bifurcation values τk,k=0,1,2,…,2K−1. For integer n≥0, we denote
J:={τ0,τ1,⋯,τ2K−1}, Jn:={τ∈J:Sn(τ)=0}, J0=J∖J0. | (3.7) |
According to Theorem 3.5, Jn has even number of bifurcation values, and the period of any periodic solution bifurcated near a bifurcation point in J0 is bounded from both below and above, while the periodic solution bifurcated near a bifurcation point in J0 has a lower bound for its period. As we shall see later in the numerical simulation, it seems impossible to find an upper bound for such period. In what follows, we will restrict our investigation on the set J0 and assume that J0≠∅. Especially, we will discuss the global continuation of periodic solutions bifurcated from the point (E∗1,τk) with τk∈J0 as the bifurcation parameter τ varies. We shall use the global Hopf bifurcation theorem for delay differential equations [19] and show that model (1.3) admits periodic solutions globally for all τ∈(τ_,ˉτ), where τ_=minJ0 and ˉτ=maxJ0.
Let x(t)=(y(τt),z(τt))T, model (1.3) can be rewritten as a general functional differential equation:
x′(t)=F(xt,τ,T),(t,τ,T)∈R+×I×R+, | (3.8) |
where xt(θ)=x(t+θ) for θ∈[−1,0], and xt∈X:=C([−1,0],R2+), and
F(xt,τ,T)=(τf(x1t(0))−pτx1t(0)x2t(0)−bτx1t(0)τe−δτh(x1t(−1))x2t(−1)−qτx1t(0)x2t(0)−μτx2t(0)) | (3.9) |
with xt=(x1t,x2t)∈X. Identifying the subspace of X consisting of all constant functions with R2+, we obtain a restricted function
˜F(x,τ,T):=F∣R2+×I×R+=(τf(x1)−pτx1x2−bτx1τe−δτh(x1)x2−qτx1x2−μτx2). |
Obviously, ˜F is twice continuously differentiable, i.e., the condition (A1) in [19] holds. We denote the set of stationary solutions of system (3.8) by E(F)={(˜x,˜τ,˜T)∈R2+×I×R+:˜F(˜x,˜τ,˜T)=0}. It follows from Proposition 2.2 that (ⅰ) if y∗1<˜y<y∗2 holds,
E(F)={(E0,τ,T),(Eb,τ,T),(E∗1,τ,T);(τ,T)∈I×R+}; |
(ⅱ) if y∗1<y∗2<˜y holds,
E(F)={(E0,τ,T),(Eb,τ,T),(E∗1,τ,T),(E∗2,τ,T);(τ,T)∈I×R+}. |
For any stationary solution (˜x,τ,T), the characteristic function is
Δ(˜x,τ,T)(λ)=λId−DF(˜x,τ,T)(eλ⋅Id). |
Note that if either y∗1<˜y<y∗2 or y∗1<y∗2<˜y holds, 0 is not a eigenvalue of any stationary solution of (3.8) and hence the condition (A2) in [19] holds. It can be checked easily from (3.9) that the smoothness condition (A3) in [19] is satisfied.
Theorem 3.5 implies that if y∗1<˜y and (H4) holds, then for each integer 0≤k≤2K−1 the stationary solution (E∗1,τk,2π/(ωkτk)) is an isolated center of (3.8), where ωk=ω(τk) is the unique positive zero of G(ω,τ) in (3.2), and there is only one purely imaginary characteristic value of the form im(2π/˜T) with m=1 and ˜T=2π/(ωkτk). Moreover, if follows from Lemma 3.4 that the crossing number ζ1(E∗1,τk,2π/(ωkτk)) at each of these centers is
ζ1(E∗1,τk,2πωkτk)=−Sign(dReλ(τ)dτ|τ=τk)={−1, 0≤k≤K−1,1, K≤k≤2K−1. | (3.10) |
Thus the condition (A4) in [19] holds. We then define a closed subset Σ(F) of X×I×R+ by
Σ(F)=Cl{(x,τ,T)∈X×I×R+:x is a nontrivial T-periodic soltuion of (3.8)}, |
and for each integer 0≤k≤2K−1, let C(E∗1,τk,2π/(ωkτk)) denote the connected component of (E∗1,τk,2π/(ωkτk)) in Σ(F). By Theorem 3.5, C(E∗1,τk,2π/(ωkτk)) is a nonempty subset of Σ(F).
To find the interval of τ in which periodic solutions exist, we shall further investigate the properties of periodic solutions of (3.8).
Lemma 3.7. All nonconstant and nonnegative periodic solutions of (3.8) are uniformly bounded. Actually, we have 0<y(t),z(t)≤M for any t∈R, where M=max{ˉy,h′(0)Mf/(pγ)}.
Proof. Since y(t) is nonconstant and nonnegative, there exists t0>0 such that y(t0)>0. Integrating the equation for y′(t) gives
e∫tt0[pz(s)+b]dsy(t)=y(t0)+∫tt0e∫rt0[pz(s)+b]dsf(y(r))dr>0,t>t0. |
Hence, y(t)>0 for all t>t0. Since y is periodic, we have y(t)>0 for all t>0. Similarly, if z(t0)>0 for some t0>0, we integrate the equation for z′(t) to obtain
e∫tt0[qy(s)+μ]dsz(t)=z(t0)+∫tt0e∫rt0[qy(s)+μ]dse−δτh(y(r−τ))z(r−τ)dr>0 |
for all t>t0. This together with periodicity of z implies that z(t)>0 for all t>0.
By Proposition 2.1, we have lim supt→∞y(t)≤ˉy. We claim y(t)≤M for all t≥0. Otherwise, if y(t1)>M for some t1≥0, then
limn→∞y(t1+nT)=y(t1)>M, |
where T is the period of the nonconstant and nonnegative periodic solutions. This contradicts the fact that y(t) is eventually bounded above by M as t→∞. Thus, M is the upper bound of y(t). Similarly, from Proposition 2.1, we have lim supt→∞z(t)≤h′(0)Mf/(pγ)≤M. This implies that z(t) is uniformly bounded by M.
Lemma 3.8. System (3.8) has no nontrivial periodic solution of period 1.
Proof. Note that any nontrivial 1-periodic solution x(t) of system (3.8) is also a nontrivial periodic solution of model (2.4), which does not admit any nontrivial periodic solution via Lemma 2.3.
We are now in the position to state our result concerning the properties of the global Hopf branches.
Theorem 3.9. Assume that (H1)-(H4), either y∗1<˜y<y∗2 or y∗1<y∗2<˜y holds, ˆτ<τmax,J0≠∅ and a0,2≤b0,2. Then for system (3.8), we have the following results:
(i) The connected component C(E∗1,τk,2π/(ωkτk)) is bounded for all τk∈J0.
(ii) Let n≥1. If Jn has 2kn bifurcation values, ordered as τn,1<⋯<τn,2kn, then for each i=1,⋯,kn, there exists j>kn such that τn,i and τn,j are connected by a global Hopf branch. Similarly, for each j=kn+1,⋯,2kn, there exists i≤kn such that τn,i and τn,j are connected by a global Hopf branch. Especially, if kn=1, then the two bifurcation values of Jn are connected.
(iii) For each τ∈(minJn,maxJn) with integer n≥1, system (3.8) has at least one periodic solution with period in (1/(n+1),1/n).
Proof. Lemma 3.7 implies that the projection of C(E∗1,τk,2π/(ωkτk)) onto X is bounded. Note that system (3.8) has no periodic solutions of period 1, and thus no periodic solutions of period 1/nk or 1/(nk+1) for any positive integer nk. It follows from Lemma 3.8 that the period T of a periodic solution on the connected component C(E∗1,τk,2π/(ωkτk)) satisfies
1nk+1<T<1nk |
with τk∈Jnk and nk≥1. Hence, the projection of C(E∗1,τk,2π/(ωkτk)) onto the T-space is bounded for τk∈J0. Note that τ∈I and I is a bounded interval. Therefore, C(E∗1,τk,2π/(ωkτk)) is bounded in X×I×R+ for any τk∈J0. This proves (ⅰ).
From the proof of Theorems 3.1 and 3.6, we know that the stationary solutions (E0,τ,T), (Eb,τ,T) and (E∗2,τ,T) cannot be a center for any τ and T if either y∗1<˜y<y∗2 or y∗1<y∗2<˜y holds and a0,2≤b0,2. It then follows from the global bifurcation theorem [19,Theorem 3.3] that E:=C(E∗1,τk,2π/(ωkτk))∩E(F) is finite and
∑(˜x,τ,T)∈Eζ1(˜x,τ,T)=0. |
If Jn has 2kn bifurcation values, ordered as τn,1<⋯<τn,2kn, then by (3.10), ζ1(˜x,τn,i,T)=−1 for each i=1,⋯,kn, and ζ1(˜x,τn,j,T)=1 for each j=kn+1,⋯,2kn. This together with the above equation implies (ⅱ).
Note that τn,1=minJn and τn,2kn=maxJn. The τ projection of the global bifurcation branch bifurcated from τn,1 includes (τn,1,τn,kn+1), while the τ projection of the global bifurcation branch bifurcated from τn,2kn includes (τn,kn,τn,2kn). Thus, for any τ∈(τn,1,τn,2k), system (3.8) has at least one periodic solution. According to the proof of (ⅰ), the period of this periodic solution lies in (1/(n+1),1/n). Thus, (ⅲ) is proved.
In this section, we present some numerical simulations to demonstrate our theoretical results. We will use the Matlab package DDE-BIFTOOL developed by Engelborghs et al. [20,21] to sketch the global Hopf branches C(E∗1,τk,2π/(ωkτk)) for 0≤k≤2K−1. Following the study of ODE model in [6], we choose
f(y)=ry(1−yK), h(y)=cy1+dy, |
and set the parameter values as follows:
r=6day−1,K=3virusmm−3,p=1mm3 cells−1 day−1,b=3day−1,c=4mm3 virus−1day−1,d=0.5mm3 virus−1,δ=0.03day−1,q=1mm3 virus−1day−1,μ=0.5day−1. | (4.1) |
Here, the units for viral loads y and immune cell density z are virus per mm3 and cells per mm3, respectively. The delay τ has the unit in days. It is easy to calculate
ˆτ=17.1939<τmax=19.1788<τb=36.6204. |
From Proposition 2.2, there are exactly two equilibria E0=(0,0) and Eb=(1.5,0) when τ>τmax, and, if τ≤τmax, there exists a third equilibrium E∗1=(y∗1,z∗1), where y∗1,z∗1 depend on τ. This is the unique positive equilibrium when τ∈[0,τe] and there exists another positive equilibrium E∗2=(y∗2,z∗2) when τ∈(τe,τmax); see Figure 1. Here, τe is the critical value when y∗2=˜y; that is,
τe=1δlnc˜y(q˜y+μ)(1+d˜y)=17.9666. |
At τ=τmax, the two positive equilibria coincide: y∗1=y∗2.
Theorem 3.1 implies that E0 is a saddle point, and Eb is locally asymptotically stable for τ≥τe, and unstable for τ<τe, Moreover, by Theorem 3.3, Eb is globally asymptotically stable in the region X1 provided that τ≥τb. It can be verified that the conditions of case (ⅱ) in Theorem 3.5 are satisfied only if 0≤τ<ˆτ. By Theorem 3.5, we know that there are exactly 4 local Hopf bifurcation values with K=2, namely,
τ0≈0.2994<τ1≈8.9274<τ2≈13.9287<τ3≈17.0322, |
as shown in Figure 2. Correspondingly,
ω0≈0.9537>ω1≈0.7540>ω2≈0.4994>ω3≈0.1069. |
Clearly, E∗1 is stable for τ∈[0,τ0)∪(τ3,τmax] and unstable for τ∈(τ0,τ3), which implies that τ0 and τ3 are stability switches. Remark that τ3<τe.
In Figure 3, we plot the four global Hopf branches C(E∗1,τk,2π/(ωkτk)) with 0≤k≤3. It is observed that the branches C(E∗1,τ1,2π/(ω1τ1)) and C(E∗1,τ2,2π/(ω2τ2)) are bounded and connected, which agrees with Theorem 3.9. The periodic solutions on these global Hopf branches are plotted on the phase plane of y and z; see Figures 4–6. As we know, the stability of the periodic solution is completely determined by the associated principal Floquet multiplier: if the principal Floquet multiplier is larger than 1, then the corresponding periodic solution is unstable, otherwise, the bifurcated periodic solution is stable [11].
By using DDE-BIFTOOL, we can further calculate the associated principal Floquet multipliers; see Figure 7. Note that the periodic solution on the first branch C(E∗1,τ0,2π/(ω0τ0)) is stable for small τ, and becomes unstable as τ increases; the periodic solution on the second branch C(E∗1,τ1,2π/(ω1τ1))–or equivalently, the third branch C(E∗1,τ2,2π/(ω2τ2))–is always unstable; the periodic solution on the fourth branch C(E∗1,τ3,2π/(ω3τ3)) is stable for τ near τ3, and becomes unstable as τ decreases. Moreover, for any τ∈(τ1,τ2), model (1.3) has at least one periodic solution on the second (or third) global Hopf branch with period in the interval (τ/2,τ), as shown in Figure 8.
It is interesting to note from Figure 3 that the first and the fourth Hopf branches do not connect. This phenomenon is quite different from that in the existing literature. We have proved that both the periodic solutions and delays on any global Hopf branch are uniformly bounded. According to the global Hopf bifurcation theory [19], either the first and the fourth Hopf branches connect, or these two branches have unbounded periods. From Figure 8, we observe that the periods of periodic oscillations on the fourth branch C(E∗1,τ3,2π/(ω3τ3)) seem to be unbounded. To further explore this interesting phenomenon, we numerically sketch in Figure 9 a bifurcation diagram for model (1.3). It is suggested that model (1.3) can have chaotic solutions and thus allow aperiodic oscillations when τ lies in the interval where the first and fourth Hopf branches disappear. Figure 9 also confirms that τ0 and τ3 are stability switches: (ⅰ) when τ crosses τ0 from left to right, the positive equilibrium E∗1 becomes unstable, while a stable periodic solution bifurcates; and (ⅱ) when τ crosses τ3 from left to right, E∗1 regains its stability while there exists a stable periodic solution when τ is close to τ3 from the left. It should be mentioned that the bifurcation diagram has a jump at τmax, which is due to the fact that there exists no positive equilibrium when τ>τmax but the positive equilibrium at τ=τmax is nonzero.
Considering the fact the immune response process involves a sequence of events such as antigenic activation, selection, and proliferation of the immune cells and the mortality of immune cells during the time lag, we derive an immunosuppressive infection model from a stage structured population model. We find that the time lag τ plays a key role in the infection and immune response dynamics. As τ increases, the model dynamics shifts through four possible outcomes: (ⅰ) virus persists without immune response when τ>τmax; (ⅱ) the outcome is initial condition dependent and delay dependent, either the immunity can be completely destroyed, or the immune mediated control can be established when τ∈(τe,τmax]; (ⅲ) immune response controls the virus when τ∈[0,τ0)∪(τ2K−1,τe); and (ⅳ) sustained oscillation appears or chaotic dynamics occurs when τ∈(τ0,τ2K−1). Our theoretical results coincide with the phenomenon of oscillatory viral loads and immune cells observed in clinical data [9].
From numerical simulation, we also observe an unbounded global Hopf bifurcation branch with bounded periodic solutions, bounded delays, but unbounded periods. To the best of our knowledge, this phenomenon is new, because, in the existing literature, either the global Hopf branch is bounded or it is unbounded with unbounded delays. A detailed theoretical study of this interesting numerical phenomenon should enrich the global Hopf bifurcation theory for delay differential equations, and we leave it as a future work.
We are very grateful to the anonymous referees for their valuable comments and suggestions, which help to improve the presentation of this manuscript. H. Shu was partially supported by the National Natural Science Foundation of China (No.11971285, No.11601392), and the Fundamental Research Funds for the Central Universities.
The authors declare there is no conflict of interest.
[1] | F. C. Bekkering, C. Stalgis, J. G. McHutchison, J. T. Brouwer, A. S. Perelson, Estimation of early hepatitis C viral clearance in patients receiving daily interferon and ribavirin therapy using a mathematical model, Hepatology, 33 (2001), 419-423. |
[2] | A. U. Neumann, N. P. Lam, H. Dahari, D. R. Gretch, T. E. Wiley, T. J. Layden, et al., Hepatitis C viral dynamics in vivo and the antiviral efficacy of interferon-α therapy, Science, 282 (1998), 103-107. |
[3] | S. A. Whalley, J. M. Murray, D. Brown, G. J. M. Webster, V. C. Emery, G. M. Dusheiko, et al., Kinetics of acute hepatitis B virus infection in humans, J. Exp. Med., 193 (2001), 847-854. |
[4] | T. W. Chun, L. Stuyver, S. B. Mizell, L. A. Ehler, J. A. Mican, M. Baseler, et al., Presence of an inducible HIV-1 latent reservoir during highly active antiretroviral therapy, Proc. Natl. Acad. Sci. USA, 94 (1997), 13193-13197. |
[5] | A. S. Perelson, Modelling viral and immune system dynamics, Nat. Rev. Immunol., 2 (2002), 28-36. |
[6] | N. L. Komarova, E. Barnes,P. Klenerman, D. Wodarz, Boosting immunity by antiviral drug therapy: a simple relationship among timing, efficacy, and success, Proc. Natl. Acad. Sci. USA, 100 (2003), 1855-1860. |
[7] | E. S. Rosenberg, M. Altfeld, S. H. Poon, M. N. Phillips, B. M. Wilkes, R. L. Eldridge, et al., Immune control of HIV-1 after early treatment of acute infection, Nature, 407 (2000), 523-526. |
[8] | A. Fenton, J. Lello, M. B. Bonsall, Pathogen responses to host immunity: the impact of time delays and memory on the evolution of virulence, Proc. R. Soc. B, 273 (2006), 2083-2090. |
[9] | G. M. Ortiz, J. Hu, J. A. Goldwitz, R. Chandwani, M. Larsson, N. Bhardwaj, et al., Residual viral replication during antiretroviral therapy boosts human immunodeficiency virus type 1-specific CD8+ T-cell responses in subjects treated early after infection, J. Virol., 76 (2002), 411-415. |
[10] | H. Shu, L. Wang, J. Watmough, Sustained and transient oscillations and chaos induced by delayed antiviral immune response in an immunosuppressive infection model, J. Math. Biol., 68 (2014), 477-503. |
[11] | J. K. Hale, S. M. V. Lunel, Introduction to Functional Differential Equations, Springer-Verlag, New York, 1993. |
[12] | X. Pan, Y. Chen, H. Shu, Rich dynamics in a delayed HTLV-I infection model: stability switch, multiple stable cycles, and torus, J. Math. Anal. Appl., 479 (2019), 2214-2235. |
[13] | H. Shu, X. Hu, L. Wang, J. Watmough, Delay induced stability switch, multitype bistability and chaos in an intraguild predation model, J. Math. Biol., 71 (2015), 1269-1298. |
[14] | J. K. Hale, Asymptotic Behavior of Dissipative Systems, AMS, Providence, 1988. |
[15] | H. Thieme, Uniform persistence and permanence for non-autonomous semiflows in population biology, Math. Biosci., 166 (2000), 173-201. |
[16] | H. Shu, L. Wang, J. Watmough, Global stability of a nonlinear viral infection model with infinitely distributed intracellular delays and CTL immune responses, SIAM J. Appl. Math., 73 (2013), 1280-1302. |
[17] | E. Beretta, Y. Kuang, Geometric stability switch criteria in delay differential systems with delay dependent parameters, SIAM J. Math. Anal., 33 (2002), 1144-1165. |
[18] | B. D. Hassard, N. D. Kazarinoff, Y. H. Wan, Theory and Applications of Hopf Bifurcation, Cambridge University Press, Cambridge, 1981. |
[19] | J. Wu, Symmetric functional differential equations and neural networks with memory, Trans. Amer. Math. Soc., 350 (1998), 4799-4838. |
[20] | K. Engelborghs, T. Luzyanina, G. Samaey, DDE-BIFTOOL v. 2.00: A Matlab package for bifurcation analysis of delay differential equations, Technical Report TW-330, K. U. Leuven, Belgium, 2001. |
[21] | K. Engelborghs, T. Luzyanina, D. Roose, Numerical bifurcation analysis of delay differential equations using DDE-BIFTOOL, ACM Trans. Math. Software, 28 (2002), 1-21. |