# Analysis of a Viral Infection Model with Delayed Nonlytic Immune Response

Authors: Mengye Chen; Liang You; Jie Tang; Shasha Su; Ruiming Zhang
Journal: Discrete Dynamics in Nature and Society
 (2015)

---

## Abstract

We investigate the dynamical behavior of a virus infection model with delayed nonlytic immune response. By analyzing corresponding characteristic equations, the local stabilities of two boundary equilibria are established. By using suitable Lyapunov functional and LaSalle’s invariance principle, we establish the global stability of the infection-free equilibrium. We find that the infection free equilibriumE0 is globally asymptotically stable when R0⩽1, and the infected equilibrium without immunity E1 is local asymptotically stable when 1<R0⩽1+bβ/cd. Under the condition R0>1+bβ/cd we obtain the sufficient conditions to the local stability of the infected equilibrium with immunity E2. We show that the time delay can change the stability of E2 and lead to the existence of Hopf bifurcations. The stabilities of bifurcating periodic solutions are studied and numerical simulations to our theorems are provided.


---

## Body

## 1. Introduction

Mathematical models have been proven valuable in understanding the population dynamics of viral load in vivo. A proper model may play significant role in a better understanding of the disease and the various drug therapy strategies. Viral infection models have received great attention in recent years [1–7]. In most viral infections, cytotoxic T lymphocytes (CTLs) play a critical role in antiviral defense by providing a cell-mediated response to specific foreign antigens associated with cells. Therefore, the immune response after a viral infection is universal and necessary to eliminate or control the disease.Recently, there have been a lot of papers on virus dynamics within host; some include the immune response directly [8–13]. During viral infections, the host immune system reacts with innate and antigen-specific immune response. Both types of response can be subdivided broadly into lytic and nonlytic components. Lytic components kill infected cells, whereas the nonlytic components inhibit viral replication through soluble mediators. As a part of the innate response, cytotoxic T lymphocytes (CTLs) kill infected cells, while antibodies neutralize free virus particles and thus inhibit the infection of susceptible cells. In addition, CD4+ and CD8+ T cells can secrete cytokines that inhibit viral replication [12]. In order to investigate the role of direct lytic and nonlytic inhibition of viral replication by immune cells in viral infections, a mathematical model was constructed to describe the basic dynamics of the interaction among susceptible host cells, a virus population, and immune response, which is described by the following differential equations [14, 15]:
(1)x′(t)=s-dx-βxy1+qz,y′(t)=βxy1+qz-ay-pyz,z′(t)=cy-bz,
where x(t), y(t), and z(t) represent the densities of uninfected target cells, virus, and CTL cells at time t, respectively. Uninfected cells are produced at rate s, die at rate dx, and become infected by virus at rate βxy without the immune response; to model nonlytic antiviral, viral replication is inhibited by the immune response at rate 1+qz; infected cells die at rate ay and are removed at rate pyz by the CTL immune response. The virus-specific CTL cells proliferate at rate cy by contact with infected cells and die at rate bz. The parameter p expresses the strength of the lytic component; the parameter q expresses the efficacy of the nonlytic component.Time delays cannot be ignored in models for immune response. Antigenic stimulation generating CTLs may need a period of timeτ [16–18]; the CTL response at time t may depend on the population of antigen at a previous time t-τ. Under the assumption of retarded immune response, Wang et al. [18] studied the effects of the time delay for immune response and assumed the time evolution of the population of CTL cells z is governed by the delayed nonlinear differential equation z′(t)=cy(t-τ)-bz. Li and Shu [19] and Xie et al. [20] investigated the effects of a time delay on a three-dimensional system with z′(t)=cy(t-τ)z(t-τ). However, both studies do not consider the nonlytic component. In this paper we propose a more general model,
(2)x′(t)=s-dx-βxy1+qz,y′(t)=βxy1+qz-ay-pyz,z′(t)=cy(t-τ)z(t-τ)-bz,
with the initial conditions
(3)(x(0),y0(θ),z0(θ))∈R+×C+×C+,
where C+=C([-τ,0],R+) which is the Banach space of continuous functions mapping the interval [-τ,0] into R+ with the topology of uniform convergence. Clearly, the models studied in [18–20] correspond to the case q=0 in our general model (2). In this work we study (2) in the more general case q⩾0.This paper is organized as follows. In Section2, we study the local asymptotic stability of the infection-free equilibrium and the immune-exhausted equilibrium of system (2), and the global asymptotic stability of the infection-free equilibrium also is investigated. In Section 3, we analyze the local stability of the positive equilibrium and the existence of Hopf bifurcations. In Section 4, the direction and stability of Hopf bifurcation are analyzed by the normal form theory and center manifold approach. In Section 5, we present numerical simulations to illustrate our results. Finally, in Section 6 we provide concluding remarks.


## 2. The Stability Analysis of Equilibrium

Considering the existence of the three equilibria, then we have the following conclusions.Proposition 1.
LetR0=sβ/ad. Then, consider the following. (1)
IfR0<1, system (2) only has an infection-free equilibrium E0=(x0,0,0), where x0=s/d.(2)
If1<R0⩽1+bβ/cd, in addition to the infection-free equilibrium E0, system (2) has another unique equilibrium E1(x1,y1,0) where x1=a/β and y1=(sβ-ad)/aβ.(3)
IfR0>1+bβ/cd, in addition to the infection-free equilibrium E0, system (2) has another infected equilibrium E2=(x2,y2,z2) that corresponds to the survival of free virus and CTL, where
(4)x2=-1dc(bpz2+ab-sc),y2=bc,z2=((acdq+bβp+cdp)2-4cdpq(abβ+acd-βcs)(acdq+bβp+cdp)2-4cdpq(abβ+acd-βcs)-acdq-bβp-cdp)×(2cdpq)-1.By the similar method in [20], we have the following result.Theorem 2.
Under the above initial conditions  (3), all solutions of system (2) are positive for t>0 and there exists M>0, such that all the solutions satisfy x(t)<M, y(t)<M, and z(t)<M for all large t.Theorem 3.
IfR0<1, then the infection-free equilibrium E0 is locally asymptotically stable, and E0 is unstable if R0>1.Proof.
The characteristic equation aboutE0 is given by
(5)|d+λsβd00a+λ-sβd000b+λ|=0;
that is,
(6)(λ+b)(λ+d)(λ+a-βsd)=0.
It is clear that (6) has negative roots λ1=-d and λ2=-b. If R0<1, then λ3=(sβ-ad)/d<(ad-ad)/d=0. This shows that all roots of (6) have the negative real part; infection-free equilibrium E0 is locally asymptotically stable. If R0>1, then λ3>0. Equation (6) has a positive root; infection-free equilibrium E0 is unstable.Theorem 4.
If1<R0⩽1+bβ/cd, then immune-exhausted equilibrium E1(x1,y1,0) is locally asymptotically stable, and E1 is unstable if R0>1+bβ/cd.Proof.
The characteristic equation aboutE1 is given by
(7)|d+y1β+λa-aqy1-y1βλpy1+aqy100b-ce-λτy1+λ|=0;
that is,
(8)(aβy1+dλ+λ2+βλy1)(b-cy1e-λτ+λ)=0.
It is easy to see that
(9)λ1+λ2=-sβa<0,λ1λ2=sβ-ad>0.
If R0>1, then λ1 and λ2 are both negative real roots, and λ3 is determined by the following equation:
(10)f(λ)≜(b-cy1e-λτ+λ)=0.
If 1<R0⩽1+bβ/cd and τ=0, then λ3=cy1-b=(-abβ-acd+βcs)/aβ<0. Therefore, if τ=0 and 1<R0⩽1+bβ/cd, then E1 is locally asymptotically stable. If τ>0, let iω (ω>0) be a solution of (10). Substituting λ=iω into (10) and separating the real and imaginary parts, we have
(11)b=cy1cos⁡(τω)ω=-cy1sin(τω).
Squaring and adding the two equations of (11) together, we obtain b2+ω2=c2y12, which has solution
(12)ω2=(c(ad-sβ)-abβ)(abβ+acd-βcs)a2β2.
Since 1<R0⩽1+bβ/cd, then we have (ad-sβ)<0 and (abβ+acd-βcs)⩾0; thus we get ω2⩽0, which contradicts with our original assumption ω>0. Therefore, (10) must have no pure imaginary root. Hence, in the case of 1<R0⩽1+bβ/cd, E1 is locally asymptotically stable for all τ⩾0. If R0>1+bβ/cd, then f(0)=(abβ+acd-βcs)/aβ<0, and lim⁡λ→+∞f(λ)=+∞. By the intermediate value theorem of continuous functions, we see that equation f(λ)=0 has at least one positive real root; thus the immune-exhausted equilibrium E1 is unstable.Theorem 5.
IfR0⩽1, the infection-free equilibrium E0(x0,0,0) is globally asymptotically stable.Proof.
Define the Lyapunov functional(13)V1(t)=x-x0-x0ln⁡xx0+y+pcz+p∫t-τt‍y(θ)z(θ)dθ;
we calculate the derivative of V1(t) along positive solutions of system (2) then substitute s for s=dx0 to obtain
(14)dV1(t)dt=(1-x0x(t))[-d(x(t)-x0)-βx(t)y(t)qz(t)+1]+βx(t)y(t)qz(t)+1-ay(t)-pbcz(t)=-d[x0-x(t)]2x(t)-bpz(t)c+y(t)[-adqz(t)+βs-ad]d(qz(t)+1)=-d[x0-x(t)]2x(t)-bpz(t)c+y(t)[-adqz(t)+ad(R0-1)]d(qz(t)+1).
Clearly, in the case of R0⩽1, we have dV1(t)/dx≤0 since x(t), y(t), and z(t) are positive, and dV1(t)/dx=0 if and only if (x,y,z)=(x0,0,0). Thus E0 is globally asymptotically stable by Lyapunov-LaSalle invariance principle.


## 3. The Stability of Infected Equilibrium and Hopf Bifurcation

Lety2/(1+qz2)=A, apply βx2/(1+qz2)=a+pz2 to the characteristic equation of system (2) at the positive equilibrium E2, and we have
(15)|λ+d+Aβa+pz2-Aq(a+pz2)-AβλA(aq+p+2pqz2)0-ce-λτz2λ-e-λτb+b|=0;
that is,
(16)P(λ,τ)+Q(λ,τ)e-λτ=0,
where
(17)P(λ,τ)=λ3+a1λ2+a2λ+a3,Q(λ,τ)=a4λ2+a5λ+a6,a1=Aβ+b+d,a2=aAβ+Abβ+Aβpz2+bd,a3=aAbβ+Abβpz2,a4=-b,a5=aAcqz2-Abβ+2Acpqz22+Acpz2-bd,a6=aAcdqz2+A2βcpqz22+A2βcpz2-Abβpz2+2Acdpqz22+Acdpz2-aAbβ.If τ=0, (16) becomes
(18)λ3+(a1+a4)λ2+(a2+a5)λ+a3+a6=0,
where
(19)a1+a4=Aβ+d>0,a2+a5=aAβ+aAcqz2+2Acpqz22+Acpz2+Aβpz2>0,a3+a6=aAcdqz2+A2βcpqz22+A2βcpz2+2Acdpqz22+Acdpz2>0.
Then,
(20)Δ1=a1+a4>0,Δ2=(a1+a4)(a2+a5)-(a3+a6)Δ2=Aβ(a+pz2)(Aβ+Acqz2+d)>0.
By the Routh-Hurwitz criterion, we know that all roots of (18) have negative real parts. From the above analysis, the following theorem holds.Theorem 6.
IfR0>1+bβ/cd, the equilibrium E2 is locally asymptotically stable for τ=0.From Theorem6, when τ=0, all roots of (16) lie to the left of imaginary axis. But, as τ is increased from zero, some of the roots may cross the imaginary axis to the right. Then the equilibrium E2 becomes unstable. Now we suppose (16) has a pure imaginary root λ=iω (ω>0). Obviously, λ=iω (ω>0) is a root of (16) if and only if ω satisfies
(21)-iω3-a1ω2+ia2ω+a3+(-a4ω2+ia5ω+a6)×(cos⁡(τω)-isin(τω))=0.
Separating the real and imaginary parts, we have
(22)a3-a1ω2=(a4ω2-a6)cos⁡(τω)-a5ωsin(τω),ω3-a2ω=(a4ω2-a6)sin(τω)+a5ωcos⁡(τω).
From (21), we also obtain
(23)|-iω3-a1ω2+ia2ω+a3|2=|-a4ω2+ia5ω+a6|2|e-iωτ|2.
So we have
(24)ω6+b1ω4+b2ω2+b3=0,
where
(25)b1=a12-a42-2a2,b2=a22-a52-2a1a3+2a4a6,b3=a32-a62.
Denote
(26)G(x)=x3+b1x2+b2x+b3.
Therefore, if (16) has a pure imaginary root iω, then equation
(27)G(x)=0
has a positive real root ω2.Suppose that (27) has K,  1⩽K⩽3, positive real roots, which are xk,  1⩽k⩽K, respectively. From (22), we have
(28)cos⁡(τxk)=Qk=-a5xk(xk3/2-a2xk)-(a3-a1xk)(a6-a4xk)-(a6-a4xk)2-a52xk,(29)sin(τxk)=Pk=xk(a4xk2-a2a4xk+a1a5xk-a6xk-a3a5+a2a6)a52xk+a42xk2-2a4a6xk+a62.
Let
(30)τk(j)={1xk[arccos(Qk)+2jπ],ifPn⩾0,1xk[2π-arccos(Qk)+2jπ],ifPn<0,
where 1⩽k⩽K and j=0,1,2,…. Then ±ixk is a pair of pure imaginary roots of (16).Theorem 7 (see [21]).
Suppose the characteristic equation is the form(31)f0(s)+f1(s)e-sτ=0,
where f0 and f1 are continuously differentiable with respect to s. One of the roots to (31) is s(τ)=α(τ)+iω(τ), where s(τ) is continuously differentiable with respect to τ, satisfying α(τ0)=0 and ω(τ0) for a positive real number τ0. Denote Φ(ω)=|f0(iω)|2-|f1(iω)|2; then we have sign⁡[dRe(s)/dτ|τ=τ0]=sign⁡[((1/2ω)·(dΦ/dω))|ω=ω0].Letλ(τ)=ξ(τ)+iω(τ) be the root of (16) near τ=τk(j) satisfying ξ(τk(j))=0, ω(τk(j))=xk. Then, by Theorem 7, we have the following theorem.Theorem 8.
d
ξ
/
d
τ
|
τ
=
τ
k
(
j
) and dG/dx|x=xk have the same sign.We chooseτ0=min⁡{τk(j)∣1⩽k⩽K,j=0,1,2,…}, where τk(j) is defined by (30). Summarizing the above analysis and combining [22], the following theorem holds.Theorem 9.
(
1
) For system (2), its equilibrium E2 is asymptotically stable for τ∈[0,τ0) if (27) has some positive real roots;
(
2
) if xk0 is a simple positive real root of (27), then system (2) undergoes a Hopf bifurcation at the equilibrium E2 when τ=τ0.


## 4. Direction and Stability of Hopf Bifurcation

In this section, we will study the direction of the Hopf bifurcations and stability of bifurcating periodic solutions by applying the normal theory and the center manifold theorem introduced in [22]. We always assume that system (2) undergoes Hopf bifurcation at the positive equilibrium E2 for τ=τ0, and then ±iω are pure imaginary roots of the characteristic equation at the positive equilibrium E2.Letu1(t)=x(τt)-x*, u2(t)=y(τt)-y*, and u3(t)=z(τt)-z*, τ=τ0+ν; then ν=0 is the Hopf bifurcation value. The Taylor expansion of (2) at the equilibrium point is given by
(32)du1(t)dt=(τ0+ν)×{-(d+βy*qz*+1)u1-βx*qz*+1u2+βqx*y*(qz*+1)2u300000+βqu1u3y*(qz*+1)2-βu1u2qz*+1+βqx*u2u3(1+qz*)2+⋯},du2(t)dt=(τ0+ν)×{βy*qz*+1u1-(py*+βqx*y*(qz*+1)2)u3000000-(p+βqx*(qz*+1)2)u2u3000000-βqy*u1u3(qz*+1)2+βu1u2qz*+1+⋯},du3(t)dt=(τ0+ν){-bu3+cz*u2(t-1)+cy*u3(t-1)000000000000000+cu3(t-1)u2(t-1)}.
We denote the above system by
(33)u˙(t)=Lν(ut)+f(ν,ut),
where u(t)=(u1(t),u2(t),u3(t))⊤∈R3, and the maps Lν:C→R3, f:R×C→R3 are given by
(34)Lν(ϕ)=(τ0+ν)×(-(d+βy*qz*+1)-βx*qz*+1βqx*y*(qz*+1)2βy*qz*+10-(py*+βqx*y*(qz*+1)2)00-b)×(ϕ1(0)ϕ2(0)ϕ3(0))+(τ0+ν)(0000000cz*cy*)(ϕ1(-1)ϕ2(-1)ϕ3(-1)),f(ν,ϕ)=(τ0+ν){(-βϕ1(0)ϕ2(0)qz*+1βϕ1(0)ϕ2(0)qz*+10)000000000l00+(βqy*ϕ1(0)ϕ3(0)(qz*+1)2-βqy*ϕ1(0)ϕ3(0)(qz*+1)20)000000000l00+(βqx*ϕ2(0)ϕ3(0)(qz*+1)2-(βqx*+p)ϕ2(0)ϕ3(0)(qz*+1)2cϕ2(-1)ϕ3(-1))+⋯}.
By the Riesz representation theorem, there exists a function η(θ,ν) of bounded variation for θ∈[0,1] such that
(35)Lν(ϕ)=∫-10‍ϕ(θ)dη(θ,ν),ϕ∈C.
In fact, we can choose
(36)η(θ,ν)=(τ0+ν)×(-(d+βy*qz*+1)-βx*qz*+1βqx*y*(qz*+1)2βy*qz*+10-(py*+βqx*y*(qz*+1)2)00-b)×δ(θ)-(τ0+ν)(0000000cz*cy*)δ(θ+1),
where δ denote the Dirac delta function. For ϕ∈C([-1,0],R3), let
(37)A(ν)ϕ(θ)={dϕ(θ)dθ,ifθ∈[-1,0),∫-10dη(ξ,ν)ϕ(ξ),‍ifθ=0;R(ν)ϕ(θ)={0,ifθ∈[-1,0),f(ν,ϕ),ifθ=0,
then system (33) is equivalent to
(38)u˙(t)=A(ν)ut(θ)+R(ν)ut(θ),
where ut(θ)=u(t+θ) for θ∈[-1,0].Forφ∈C([-1,0],R3), the adjoint operator of A(0) is defined by
(39)A*(0)φ(s)={-dφ(s)ds,ifs∈(0,1],∫-10‍dη⊤(ξ,0)φ(-ξ),ifs=0,
and define a bilinear inner product
(40)〈φ(s),ϕ(θ)〉=φ¯⊤(0)ϕ(0)-∫-10‍∫ξ=0θ‍φ¯⊤(ξ-θ)dη(θ)ϕ(ξ)dξ,
where η(θ)=η(θ,0). By the discussion in Section 3 we know that ±iωτ0 are eigenvalues of A(0). Hence, they are also eigenvalues of A*. We define q(θ) and q*(s) to be the eigenvectors of A and A* corresponding to the eigenvalues iωτ0 and -iωτ0, respectively.Assume thatq(θ)=(1,q1,q2)⊤eiωτ0θ is the eigenvector of A(0) corresponding to iωτ0; then A(0)q(θ)=iωτ0q(θ). It follows from the definition of A(0) that we have
(41)τ0(iω+(d+βy*qz*+1)βx*qz*+1-βqx*y*(qz*+1)2-βy*qz*+1iω(py*+βqx*y*(qz*+1)2)0-cz*e-iωτ0iω+b-cy*e-iωτ0)q(0)=(000).
Then we obtain
(42)q1=-((qz*+1)(p(qz*+1)(βy*+(d+iω)(qz*+1))(βy*+(d+iω)(qz*+1))+βqx*(d+iω))×(βx*(βqx*+(qz*+1)(pqz*+p+iqω)))-1,q2=(i(qz*+1)(ω(d+iω)(qz*+1)2+βy*(qωz*-iβx*+ω))+βy*(qωz*-iβx*+ω))hhhhhhhhhhhh+βy*(qωz*-iβx*+ω)))×(βx*y*(βqx*+(qz*+1)(pqz*+p+iqω)))-1.
Similarly, let q*(s)=D(1,q1*,q2*)⊤eiωτ0s be the eigenvector of A* corresponding to -iωτ0; by the definition of A* we get
(43)q1*=βy*+(d-iω)(qz*+1)βy*,q2*=-(ie-iτ0ω(ω(d-iω)(qz*+1)2+βy*(βcy*z*(qz*+1))-1hhhhhhhhhhhhh×(qωz*+iβx*+ω)ie-iτ0ω)hhhhhlh×(βcy*z*(qz*+1))-1).
In order to assure 〈q*(s),q(θ)〉=1, we need to determine the value of D. From (40) we have
(44)〈q*(s),q(θ)〉=D¯(1,q1*¯,q2*¯)(1,q1,q2)⊤-∫-10‍∫ξ=0θ‍D¯(1,q1*¯,q2*¯)e-iωτ0(ξ-θ)dη(θ)000000000×(1,q1,q2)⊤eiωτ0ξdξ=D¯{1+q1q1*¯+q2q2*¯-∫-10‍(1,q1*¯,q2*¯)θeiωτ0θ000×∫-10‍[dη(θ)](1,q1,q2)⊤}=D¯{(0000000cz*cy*)(1+q1q1*¯+q2q2*¯)+τ0e-iωτ0(1,q1*¯,q2*¯)000×(0000000cz*cy*)(1,q1,q2)⊤}=D¯{(1+q1q1*¯+q2q2*¯)+cτ0e-iωτ0q2*¯000×(1+q1q1*¯+q2q2*¯)(z*q1+y*q2)}=1.
Hence
(45)D¯=[(1+q1q1*¯+q2q2*¯)+cτ0e-iωτ0q2*¯(z*q1+y*q2)]-1.
Using the method in [22] we obtain the following coefficients:
(46)g20=2τ0D¯(F21q1*¯+F31q2*¯+F11),g11=τ0D¯(F22q1*¯+F32q2*¯+F12),g02=2τ0D¯(F23q1*¯+F33q2*¯+F13),g21=2τ0D¯(F24q1*¯+F34q2*¯+F14),
where
(47)F11=β(-q2)q22Q3x*y*+Q2×(βqq1q2x*+βqq2y*)-βq1Q,F12=-2βq2q2Q3x*y*q2¯+Q2×(βqq1x*q2¯+βqq2x*q1¯+βqy*q2¯+βqq2y*)+Q(-βq1¯-βq1),F13=β(-q2)Q3x*y*q2¯2+Q2×(βqx*q1¯q2¯+βqy*q2¯)-βQq1¯,F14=3βq3q22Q4x*y*q2¯+Q3×(-βq22q2x*q1¯-2βq1q2q2x*q2¯-2βq2q2y*q2¯00000-βq22q2y*-βq2W320(0)x*y*q2¯00000-2βq2q2W311(0)x*y*)+Q2(12βqq1q2¯+βqq2q1¯+βqq2q10000000+12βqW120(0)y*q2¯0000000+βqq2W111(0)y*+12βqW220(0)x*q2¯0000000+βqq2W211(0)x*0000000+βqq1W311(0)x*+12βqW320(0)x*q1¯0000000+βqW311(0)y*+12βqW320(0)y*)+Q(12-βW211(0)-βq1W111(0)00000000-12βW220(0)-12βW120(0)q1¯),F21=-pq1q2+βq2q22Q3x*y*+Q2×(β(-q)q1q2x*-βqq2y*)+βq1Q,F22=-pq1q2¯-pq2q1¯+2βq2q2Q3x*y*q2¯+Q(βq1¯+βq1)+Q2(β(-q)q1x*q2¯-βqq2x*q1¯000000-βqy*q2¯-βqq2y*),F23=-pq1¯q2¯+βq2Q3x*y*q2¯2+Q2×(β(-q)x*q1¯q2¯-βqy*q2¯)+βQq1¯,F24=-12pW220(0)q2¯-3βq3q22Q4x*y*q2¯-pq2W211(0)-12pW320(0)q1¯-pq1W311(0)+Q(12βW120(0)q1¯+βq1W111(0)0000+βW211(0)+12βW220(0))+Q2(-βqq2q1-βqq2q1¯-12βqW120(0)y*q2¯000000-βqq2W111(0)y*-βqq2W211(0)x*000000-12βqW220(0)x*q2¯-12βqW320(0)x*q1¯000000-βqq1q2¯+β(-q)q1W311(0)x*000000-βqW311(0)y*-12βqW320(0)y*)+Q3(βq22q2x*q1¯+2βq1q2q2x*q2¯00000+2βq2q2y*q2¯+βq22q2y*00000+2βq2q2W311(0)x*y*00000+βq2W320(0)x*y*q2¯),F31=cq1q2e-2iτ0ω,F32=cq1q2¯+cq2q1¯,F33=cq1¯q2¯e2iτ0ω,F34=12cW220(-1)q2¯eiτ0ω+12cW320(-1)q1¯eiτ0ω+cq2W211(-1)e-iτ0ω+cq1W311(-1)e-iτ0ω,Q=11+qz*,W120=W20(1),W220=W20(2),W320=W20(3),W111=W11(1),W211=W11(2),W311=W11(3),W20(θ)=ig20ωτ0q(0)eiθωτ0+ig02¯3ωτ0q(0)¯e-iθωτ0+E1e2iθωτ0,W11(θ)=-ig11ωτ0q(0)eiθωτ0+ig11¯ωτ0q(0)¯e-iθωτ0+E2,E1=(d+2iω+βy*qz*+1βx*qz*+1-qβx*y*(qz*+1)2-βy*qz*+12iωpy*+qβx*y*(qz*+1)20-ce-2iωτ0z*b+2iω-ce-2iωτ0y*)-1×(F11F21F31),E2=(d+βy*qz*+1βx*qz*+1-qβx*y*(qz*+1)2-βy*qz*+10py*+qβx*y*(qz*+1)20-cz*b-cy*)-1(F12F22F32).
Then we get
(48)c1(0)=i2ωτ0(g20g11-2|g11|2-|g02|3)+g212,ν2=-Re{c1(0)}Re{λ′(τ0)},β2=2Re{c1(0)},T2=-Re{c1(0)}+ν2Re{λ′(τ0)}ωτ0,
which determine the quantities of bifurcating periodic solutions reduced on the central manifold at the critical value τ0. More precisely, ν2 determines the direction of Hopf bifurcation and β2 determines the stability of bifurcating periodic solution. Summarizing the above discussion, we have the following main result.Theorem 10.
Ifν2>0 (ν2<0), then the Hopf bifurcation is supercritical (subcritical) and the bifurcating periodic solutions exist for τ>τ0 (τ<τ0). Moreover, β2 determines the stability of the bifurcating periodic solutions. The bifurcating periodic solutions are stable (unstable) if β2<0 (β2>0).


## 5. Numerical Simulation

In order to demonstrate our results and find complex dynamic behavior of system (2), we provide numerical calculations for different birth rates of susceptible cells s and time delays τ. The parameter values are chosen from literatures [14, 15, 20]. If we choose the following data set: s=200, β=0.002, q=0.01, a=5, b=0.1, c=0.2, d=0.1, and p=0.1, then we get R0=0.8<1, E0=(2000,0,0). Hence, by Theorem 5, E0 is globally asymptotically stable. Figure 1 illustrates this fact.Figure 1
The phase diagram of system (2) with different initial values and different τ.If we choose the following data set:s=251, β=0.002, q=0.01, a=5, b=0.1, c=0.2, d=0.1, and p=0.1, then we get 1<R0⩽1.01, E1=(2500,0.3,0). Hence, by Theorem 4, E2 is local asymptotically stable. Figure 2 illustrates this fact.Figure 2
Local asymptotic stability of the immune-exhausted equilibriumE1.In addition, we choose a set of parameters:s=257, β=0.002, q=0.01, a=5, b=0.1, c=0.2, d=0.1, and p=0.1, then we get R0>1.01, E2=(2544.7,0.5,0.5937), the positive real roots of (27), x1=0.00607941 and x2=0.00194817. It satisfies the conditions (2) of Theorem 9. We can calculate that
(49)c1(0)=-0.00241449-0.003121i,τ1(0)=8.78303,τ2(0)=37.8851,τ1(1)=89.3671,τ1(2)=169.951.By Theorem9, we have that E2 is stable when τ<τ1(0) (Figure 3). Let τ0=τ1(0); we can compute c1(0)=-0.000199515-0.0139977i, G′(x1)=0.000059>0. Since sign⁡[Reλ′(τ0)]=sign⁡[G′(x1)] and Re{c1(0)}<0, then ν2>0; according to Theorem 10, there exist stable periodic solutions (Figure 4) for τ>τ1(0) in a τ1(0) neighborhood. When τ exceeds τ2(0), we can observe that another change of stability appears (Figure 5). Since G′(x2)=-0.00004<0 and Re{c1(0)}<0 then ν2>0. Applying Theorem 10, there exist stable periodic solutions for τ<τ2(0). By the same method, the stability of E2 is changed when τ exceeds τ1(1), τ1(2). We can also obtain that there exist stable periodic solutions for τ>τ1(1)(Figure 6). The bifurcation diagram (Figure 7) shows that the stability of E2 is changed at τ1(0), τ2(0), and τ1(1) clearly.Figure 3
Stability ofE2 for τ=5.Figure 4
A stable bifurcating periodic solution forτ=13.Figure 5
Stability ofE2 for τ=40.Figure 6
A stable bifurcating periodic solution forτ=95.Figure 7
A bifurcation diagram in whichτ is bifurcating parameter.As shown in Figures3–7, for different delays in CTL response, the system is stable at a positive equilibrium when the delay is small or it is at a stable periodic oscillation when the delay is large. Our simulations demonstrated that delays may affect cell-virus and virus-immune system dynamics. These results could provide insights into many important questions such as how to further improve our mathematical model, how to better implement clinical diagnosis, and how to design and administrate more efficient medical treatment.


## 6. Conclusion

In this paper, we have studied a virus infection model with delayed nonlytic immune response. We obtain the sufficient and necessary conditions for the existence of the equilibria. Global stability of the infection-free equilibrium has been given by the Lyapunov-LaSalle type theorem. We find that the infection-free equilibriumE0 is globally asymptotically stable if R0⩽1. In this case, the virus is finally cleared up. By analyzing corresponding characteristic equation, we obtain that the immune-exhausted equilibrium E1 is locally asymptotically stable if 1<R0⩽1+bβ/cd. In this case, the infection becomes chronic but without CTL immune response. By choosing time delay as a bifurcation parameter, we have shown that a limit cycle emerges via Hopf bifurcation when the delay passes through the critical value τ0. By applying normal form theory and center manifold theorem, the direction of Hopf bifurcation and stability of bifurcating periodic solutions have been investigated. Some numerical simulations are performed to verify the correctness of theoretical analysis. Clearly, compared to the earlier studies, our analysis shows that the new time delay we introduce can change the stability of the equilibrium E2. We would like to mention that whether the equilibrium E2 is globally asymptotically stable when R0>1+bβ/cd. We leave this for our future work.



---

*Source: 235420-2015-05-14.xml*