# Modelling the Drugs Therapy for HIV Infection with Discrete-Time Delay

Authors: Xueyong Zhou; Xiangyun Shi
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

A discrete-time-delay differential mathematical model that described HIV infection of CD4+ T cells with drugs therapy is analyzed. The stability of the two equilibria and the existence of Hopf bifurcation at the positive equilibrium are investigated. Using the normal form theory and center manifold argument, the explicit formulas which determine the stability, the direction, and the period of bifurcating periodic solutions are derived. Numerical simulations are carried out to explain the mathematical conclusions.


---

## Body

## 1. Introduction

Recently there has been a substantial effort in the mathematical modelling of virus dynamics [1–8]. These models focus on uninfected target cells, infected cells that are producing virus, and virus. A basic mathematical model describing HIV infection dynamic model is of the following form which has been studied in [5, 9]:
(1)dT(t)dt=s-μ1T(t)-kT(t)V(t),dI(t)dt=kT(t)V(t)-δI(t),dV(t)dt=NδI(t)-cV(t).In system (1), the following variables are includes: T(t) uninfected cells at time t (unit is cells mm-3), I(t) infected cells at time t (unit is cells mm-3), and V(t) virus at time t (unit is virions mm-3). Parameters μ1,δ, and c are the death rates of the uninfected T cells, the infected T cells, and the virus particles, respectively. k is the contact rate between uninfected T cells and the virus particles. N is the average number of virus particles produced by an infected T cell.Reverse transcriptase inhibitors (RTIs) are a class of antiretroviral drugs used to treat HIV infection. RTIs inhibitors work by inhibiting the action of reverse transcriptase. RTIs inhibit the activity of reverse transcriptase, a viral DNA polymerase enzyme that retroviruses need to reproduce. In [10], Srivastava et al. developed a mathematical model for primary infection with RTIs. They subdivided the infected cells class in two subclasses: pre-RT (denoted by I1(t)) and post-RT (denoted by I2(t)). They assumed that a virus enters a resting CD4+ T cell, the viral RNA may not be completely reverse transcribed into DNA, the unintegrated virus may decay with time and partial DNA transcripts are labile and degrade quickly [11, 12]. And they also assumed that a fraction of cells ηaI1(t) in pre-RT class reverts back to uninfected class and the remaining (1-η)aI1(t) proceeds to post-RT class and becomes productively infected due to presence of RT inhibitors. The model of Srivastava et al. is as follows
(2)dT(t)dt=s-μ1T(t)-kT(t)V(t)+(ηa+b)I1(t),dI1(t)dt=kT(t)V(t)-(d+a+b)I1(t),dI2(t)dt=(1-η)aI1(t)-δI2(t),dV(t)dt=NδI2(t)-cV(t),
where 0<η<1 is the efficacy of reverse transcriptase inhibitors (RTIs), a is the transition rate from pre-RT (i.e., I1(t)) infected CD4+ T cells class to productively post-RT (i.e., I2(t)) which is a productively infected class, and b is the reverting rate of infected cells to uninfected class due to noncompletion of reverse transcription [11, 12].Protease inhibitors (PIs) are a class of drugs used to treat or prevent infection by viruses, including HIV and hepatitis C. PIs prevent viral replication by inhibiting the activity of HIV-1 protease, an enzyme used by the viruses to cleave nascent proteins for final assembly of new  virus. The new virous are noninfectious. Virions that were created prior to drug treatment remain infectious. Thus, in the presence of a protease inhibitor, two types of virus particles (i.e., infectious virions and noninfectious virions) should be considered [5]. We need the drug to be highly effective if we use single drug to treat. Hence, combination anti-HIV therapy is now the standard of care for people with HIV. So far as we know, there are few mathematical models about the effects of combination anti-HIV therapy [7, 13]. Therefore, considering the effects of both RTIs and PIs, model (2) can be modified to
(3)dT(t)dt=s-μ1T(t)-kT(t)V1(t)+(ηa+b)I1(t),dI1(t)dt=kT(t)V1(t)-(d+a+b)I1(t),dI2(t)dt=(1-η)aI1(t)-δI2(t),dV1(t)dt=(1-p)NδI2(t)-cV1(t),dV2(t)dt=pNδI2(t)-cV2(t),
where variables V1(t) and V2(t) denote infectious and non-infectious virus at time t, respectively. And V(t)=V1(t)+V2(t) is the total virus concentration at time t. Parameter p∈[0,1] denotes the effectiveness of PIs with p=1 meaning that the therapy with PIs is 100% effective and no newly infectious virus particles will be produced [5].In the real situation, there may be a delay between the time target cells which are contacted by the virus particles and the time the contacted cells become actively affected meaning that the contacting virions enter cells. Hence, time delays of one type or another have been incorporated into viral dynamical models by many authors. The first model that included this type “intracellular” delay was developed by Herz et al. [14] and assumed that cells became productively infected time units after HIV initial infection. Nelson et al. [15] extend the development of delay models of HIV infection and treatment to the general case of combination antiviral therapy that is less than completely efficacious. Recently, in studying the viral clearance rates, Perelson et al. [9] assumed that there are two types of delays that occur between the administration of drug and the observed decline in viral load: a pharmacological delay that occurs between the ingestion of drug and its appearance within cells and an intracellular delay that is between initial infection of a cell by HIV and the release of new virions. Furthermore, the growth of CD4+ T cells in humans is not well understood.Recently, studies in various fields such as biology, control, economy, chemistry, and electrodynamics have shown that delay differential equations play an important role in explaining many different phenomena [16–20]. Srivastava et al. [10] proposed and analyzed a mathematical model for the effect of RTIs on the dynamics of HIV. In [21], Culshaw and Ruan have considered that the basic model of HIV infection in host was extended to incorporate logistic growth and an intracellular delay. However, none of these models have incorporated antiretroviral therapy, logistic growth of the CD4+ T cell, and intracellular delay. Here, we build on the basic model of HIV pathogenesis in host, adding the effects of antiretroviral therapy, logistic growth of the CD4+ T cell, and intracellular delay. Hence, we can obtain the following model:
(4)dT(t)dt=s+rT(t)(1-T(t)Tmax⁡)-μ1T(t)-kT(t)V1(t)+(ηa+b)I1(t),dI1(t)dt=kT(t-τ)V1(t-τ)-(d+a+b)I1(t),dI2(t)dt=(1-η)aI1(t)-δI2(t),dV1(t)dt=(1-p)NδI2(t)-cV1(t),dV2(t)dt=pNδI2(t)-cV2(t).
In model (4), T(t), I1(t), I2(t), V1(t), and V2(t) represent the density of susceptible CD4+ T cells, infected CD4+ T cells before reverse transcription (i.e., those infected cells which are in pre-RT class), infected CD4+ T cells in which reverse transcription is completed (post-RT class), infectious virus, and noninfectious virus at time t, respectively. The meaning of the parameters are as follows: s is the source term for uninfected CD4+ T cell, k is the rate at which CD4+ T cell becomes infected with virus, μ1 is the death rate of healthy CD4+ T cell, η is the efficacy of RTIs, a is the transition rate from pre-RT infected CD4+ T cells to productively post-RT, b is the reverting rate of infected cells to uninfected class, d is the death rate of infected T cells, δ is the death rate of actively infected T cells I2, N is the number of virions produced by infected CD4+ T cells, c is the clearance rate of virus, r is the maximum proliferation rate, Tmax⁡ is the T cell population density at which proliferation shuts off, p is the efficacy of protease inhibitor, and τ is the “intracellular” delay.Note that the non-infectious HIV virusV2(t) does not appear in the first four equations of system (4). Thus, we can consider the following subsystem of system (4):(5)dT(t)dt=s+rT(t)(1-T(t)Tmax⁡)-μ1T(t)-kT(t)V1(t)+(ηa+b)I1(t),dI1(t)dt=kT(t-τ)V1(t-τ)-(d+a+b)I1(t),dI2(t)dt=(1-η)aI1(t)-δI2(t),dV1(t)dt=(1-p)NδI2(t)-cV1(t).In this paper, we will discuss the dynamics of model (5). This paper is organized as follows. In Section 2, we present some preliminaries about system (5), for example, the positivity of solutions and the expression of equilibria. We discuss the local stability of the uninfected equilibrium in Section 3. In Section 4, we discuss the local stability and Hopf bifurcation at the infected equilibrium. In Section 5, the direction and stability of the local Hopf bifurcation are established. In Section 6, some numerical simulations are performed to illustrate the analytical results found. A brief discussion is presented in the last section.


## 2. Preliminaries

System (5) is a system of delay differential equations. For such a system, initial functions need to be specified and well-posedness needs to be addressed. We denote by ℂ the Banach space of continuous functions φ:[-τ,0]→ℝ4 with norm
(6)∥φ∥=sup⁡-τ≤ς≤0{|φ1(ς)|,|φ2(ς)|,|φ3(ς)|,|φ4(ς)|},
where φ=(φ1,φ2,φ3,φ4). As usual, the initial condition of (5) is given as
(7)T(ς)=φ1(ς),I1(0)=φ2(0),I2(0)=φ3(0),V1(ς)=φ4(ς),ς∈[-τ,0],
where the initial function φ=(φ1,φ2,φ3,φ4) belongs to the Banach space ℂ=ℂ([-τ,0],ℝ4) of continuous functions mapping the initial [-τ,0] into ℝ4. For biological reasons, the initial functions are assumed as
(8)T(ς)=φ1(ς)≥0,ς∈[-τ,0],φ1(0)>0;I1(0)=φ2(0)>0;I2(0)=φ3(0)>0;V1(ς)=φ4(ς)≥0,ς∈[-τ,0],φ4(0)>0.
In this paper, we will discuss the dynamical behavior of system (5) with the initial conditions in (8). By the fundamental theory of functional differential equations [22], we know that there is a unique solution (T(t),I1(t),I2(t),V1(t)) to system (5) with initial conditions (8).Firstly, we present the positivity of the solutions. System (5) can be put into the matrix form
(9)X˙(t)=G(X(t)),whereX(t)=(T(t),I1(t),I2(t),V1(t))⊤∈ℝ4 and G(X(t)) is given by(10)G(X(t))=(G1(X(t))G2(X(t))G3(X(t))G4(X(t)))=(s+rT(t)(1-T(t)Tmax⁡)-μ1T(t)-kT(t)V1(t)+(ηa+b)I1(t)kT(t-τ)V1(t-τ)-(d+a+b)I1(t)(1-η)aI1(t)-δI2(t)(1-p)NδI2(t)-cV1(t)).Letℝ+4=[0,+∞)×[0,+∞)×[0,+∞)×[0,+∞) be the nonnegative octant in ℝ4;  G:ℝ+4+1→ℝ4, G∈ℂ∞(ℝ4) (where G is a function of the variable X(t)∈ℝ+4) is locally Lipschitz and satisfies the condition
(11)Gi(X(t))|xi(t)=0,X(t)∈ℝ+4≥0,
where x1(t)=T(t),x2(t)=I1(t),x3(t)=I2(t), and x4(t)=V1(t).Due to lemma in [23] any solution of (9) with X(ς)∈ℂ+, say X(t)=X(t,X(ς)), is such that X(t)∈ℝ+4 for all t≥0.System (5) has an uninfected (boundary) equilibrium and an infected (positive) steady state. The uninfected equilibrium is E0(T0,0,0,0), where
(12)T0=Tmax⁡2r[r-μ1+(r-μ1)2+4rsTmax⁡].The infected equilibrium isE*(T*,I1*,I2*,V1*), where
(13)T*=c(d+a+b)(1-p)(1-η)kNa,I1*=1d+(1-η)a[s-dT*+rT*(1-T*Tmax⁡)],I2*=(1-η)aδ[d+(1-η)a][s-dT*+rT*(1-T*Tmax⁡)],V1*=(1-p)(1-η)Nac[d+(1-η)a][s-dT*+rT*(1-T*Tmax⁡)].The basic reproductive number is given asℛ0=T0/T*. The basic reproductive number ℛ0 measures the average number virus-producing target cells produced by an single virus-producing target cell during its entire infectious period in an entirely uninfected targeT cell population [24, 25]. It is easy to see that ℛ0>1 ensures the existence of the infected equilibrium E*.


## 3. Stability of Uninfected EquilibriumE0

In this section, we will discuss the stability of the uninfected equilibriumE0(T0,0,0,0).LetE-(T-,I-1,I-2,V-1) be any arbitrary equilibrium. To study the stability of the steady state E-, let us define
(14)x(t)=T(t)-T-,y1(t)=I1(t)-I-1,y2(t)=I2(t)-I-2,z(t)=V1(t)-V-1.
Then, the linearized system of (5) around the equilibrium E- is given by
(15)ddt(x(t)y1(t)y2(t)z(t))=A1(x(t)y1(t)y2(t)z(t))+A2(x(t-τ)y1(t-τ)y2(t-τ)z(t-τ)),
where A1 and A2 are 4×4 matrices given by(16)A1=(-μ1+r-2rT-Tmax⁡-kV-1ηa+b0-kT-0-a-b-d000(1-η)a-δ000(1-p)Nδ-c),A2=(0000kV-100kT-00000000).Hence, the characteristic equation of system (5) at E- is given by
(17)det⁡(A+Be-λτ-λℐ)=0,
where ℐ is a 4×4 identity matrix that is,(18)|-μ1+r-2rT-Tmax⁡-kV-1-ληa+b0-kT-kV-1e-λτ-a-b-d-λ0kT-e-λτ0(1-η)a-δ-λ000(1-p)Nδ-c-λ|=0.Theorem 1.
(1) Ifℛ0<1, E0 is locally asymptotically stable for any time delay τ≥0. (2) If ℛ0>1, E0 is unstable for any time delay τ≥0. (3) If ℛ0=1, it is a critical case.Proof.
For uninfected equilibriumE0, (18) reduces to
(19)(r-μ1-2rT0Tmax⁡-λ)[λ3+b1λ2+b2λ+b3+c3e-λτ]=0,
where
(20)b1=a+b+d+δ+c,b2=(a+b+d)(c+δ)+cδ,b3=(a+b+d)cδ,c3=-(1-p)(1-η)akNδT0.
It is clear that (19) has the characteristic root λ1=r-μ1-(2rT0/Tmax⁡)=-(r-μ1)2+4rs/Tmax⁡<0.
Next, we will consider the transcendental polynomial(21)λ3+b1λ2+b2λ+b3+c3e-λτ=0.
Forτ=0, we have that
(22)λ3+b1λ2+b2λ+b3+c3=0.
Obviously, b1>0, b2>0, and b3+c3>0 since ℛ0<1. We also get
(23)b1b2-(b3+c3)=(a+b+d)2(c+δ)+(a+b+d)(c+δ)2+cδ(c+δ)+(1-p)(1-η)akNδT0>0.
This shows that all the roots of (22) have negative real parts for τ=0 by using Routh-Hurwitz theorem.
In the following, we investigate the existence of purely imaginary rootsλ=iω, ω>0, of (21). If τ>0 and λ=iω with ω>0 is a solution of (21), then separating the real and imaginary parts gives
(24)ω3-b2ω=-c3sin⁡(ωτ),b1ω2-b3=c3cos⁡(ωτ).
Squaring and adding both equations of (24) yields
(25)f(ω,τ)=ω6+m1ω4+m2ω2+b32-c32=0,
where
(26)m1=(a+b+d)2+c2+δ2>0,m2=(a+b+d)2(c2+δ2)+(a+b+d)2cδ+(cδ)2+cδ(c+δ)(a+b+d)>0.
Lettingy=ω2 yields
(27)y3+m1y2+m2y+b32-c32=0.
Ifℛ0<1, then b32-c32>0. Therefore, by claim 1 in [21], it is evident that (27) has no positive real roots. This shows that (21) cannot have a purely imaginary root for any τ>0. Therefore, the uninfected equilibrium E0 is locally asymptotically stable for any τ≥0 provided that ℛ0<1.
Ifℛ0=1, the transcendental polynomial (21) becomes
(28)λ3+b1λ2+b2λ+b3-b3e-λτ=0.
It is clear that λ=0 is a simple root of (28). We further show that any root of (28) must have negative real part except λ=0.
In fact, if (28) has imaginary root u±iω for some u≥0, ω≥0, and τ≥0, from (28) we have
(29)u3-3uω2+b1u2-b1ω2+b2u+b3=b3e-uτcos⁡⁡(ωτ),-ω3+3u2ω+b2ω+2b1uω=-b3e-uτsin(ωτ),
which, together with u≥0, implies that
(30)[u3-3uω2+b1u2-b1ω2+b2u+b3]2+[-ω3+3u2ω+b2ω+2b1uω]2=b32e-2uτ≤b32.
However, it is easy to check that the previous inequality is not true. Hence, it shows that any root of (28) has negative real part except λ=0, which implies that the trivial solution of (5) is stable for any time delay τ≥0.
Ifℛ0>1, let
(31)f(λ)=λ3+b1λ2+b2λ+b3+c3e-λτ=0.
Note that f(0)=b3+c3<0 since R0>1 and lim⁡λ→+∞⁡f(λ)=+∞. It follows from the continuity of the function f(λ) on (-∞,+∞) that equation f(λ)=0 has at least one positive root. Hence, characteristic equation (19) has at least one positive. Thus, E0 is unstable. Therefore, our results in this theorem are proved.


## 4. Dynamical Behavior of Endemic EquilibriumE*

In general, the nonlinear delay system will undergo a Hopf bifurcation when the delay passes through a critical value of the delay, for which the stability of the existing equilibrium changes from stable status to unstable status and a self-excited limit cycle emerges at this moment. Under certain conditions, the existence of a Hopf bifurcation can be determined from linear stability analysis; it requires that at the bifurcation point, the characteristic function has exactly one pair of conjugate roots on the imaginary axis, and as the delay passes through the bifurcation point, this pair of characteristic roots cross from the left-half complex plane to the right-half complex plane or vice verse [19, 26]. The crossing direction is the same as that mentioned previously in linear stability analysis. Thus, the determination of the crossing direction is very important for both stability analysis and Hopf bifurcation. In this section, we will consider the dynamical behavior of endemic equilibrium E*. Some conditions for Hopf bifurcation around equilibrium E* to occur are obtained by using the time delay τ as a bifurcation parameter.For endemic equilibriumE*(T*,I1*,I2*,V1*), (18) reduces to(32)|-μ1+r-2rT*Tmax⁡-kV1*-ληa+b0-kT*kV1*e-λτ-a-b-d-λ0kT*e-λτ0a(1-η)-δ-λ000(1-p)Nδ-c-λ|=0;that is,
(33)λ4+p1λ3+p2λ2+p3λ+p4-(q2λ2+q3λ+q4)e-λτ=0,
where
(34)p1=a+b+d+c+δ-Ω,p2=-Ω(a+b+d)+(c+δ)(a+b+d-Ω)+cδ,p3=-(c+δ)Ω(a+b+d)+cδ(a+b+d-Ω),p4=-cδΩ(a+b+d),q2=(ηa+b)kV1*,q3=(c+δ)(ηa+b)kV1*+(1-p)(1-η)NδakT*,q4=cδ(ηa+b)kV1*-(1-p)(1-η)NδakT*(Ω+kV1*),Ω=-μ1+r-2rT*Tmax⁡-kV1*<0.Obviously,p1>0. In addition, in view of Routh-Hurwitz criteria, we can easily know that all roots of (33) with τ=0 have negative real parts if the following condition holds:
(35)(H):p4-q4>0,p3-q3>0,p1[(p3-q3)(p2-q2)-p1(p4-q4)]>(p1-q3)2.Let us considerτ≠0 and assume λ(τ)=ϕ(τ)+iψ(τ), where ϕ(τ),ψ(τ)∈R. Substituting λ(τ)=ϕ(τ)+iψ(τ) and rewriting (33) in terms of its real and imaginary parts, we obtain(36a)υ4+ω4-6υ2ω2+p1(υ3-3υω2)+p2(υ2-ω2)+p3υ+p4=e-τυ{q2[(υ2-ω2)cos⁡(τω)+2υωsin(τω)]+q3[υcos⁡(τω)+ωsin(τω)]+q4cos⁡(τω)(υ2-ω2)},(36b)4υω(υ2-ω2)+p1(3υ2ω-ω3)+p2(2υω)+p3ω=e-τυ{q2[(-υ2+ω2)sin(τω)+2υωcos⁡(τω)]+q3[-υsin(τω)+ωcos⁡(τω)]+q4[-sin(τω)](-υ2+ω2)}.Letτ1* be such that υ(τ1*)=0 and ω(τ1*)=ω(τ*); then (36a) and (36b) reduce to(37a)ω1*4-p2ω1*2+p4=(-q2ω1*2+q4)cos⁡(τ1*ω1*)+q3ω1*sin(τ1*ω1*),(37b)-p1ω1*3+p3ω1*=q3ω1*cos⁡(τ1*ω1*)+(q2ω1*2-q4)sin(τ1*ω1*).Eliminating τ, we have
(38)ω1*8+(p12-2p2)ω1*6+(p22+2p4-2p1p3-q22)ω1*4+(p32-2p2p4+2q2q4-q32)ω1*2+(p42-q42)=0.Suppose thatω1* is the last positive simple root of (38). We will now show that, with this value of ω1*, there is a τ1* such that υ(τ1*)=0 and ω(τ1*)=ω1*. Given ω1*, (37a) and (37b) can be written as(39a)U=Φcos⁡(τ1*ω1*)+Ψsin(τ1*ω1*),(39b)V=Ψcos⁡(τ1*ω1*)-Φsin(τ1*ω1*),where
(40)Φ=-q2ω1*2+q4,Ψ=q3ω1*,U=ω1*4-p2ω1*2,V=-p1ω1*3+p3ω1*,Φ2+Ψ2=U2+V2=H2,
where H>0.Equations(41)Φ=Hcos⁡θ,Ψ=Hsinθ
determine a unique θ∈[0,2π). With this value of θ,(42a)Hcos⁡(τ1*ω1*)cos⁡θ+Hsin(τ1*ω1*)sinθ=U,(42b)Hcos⁡(τ1*ω1*)sinθ-Hsin(τ1*ω1*)cos⁡θ=V.Hence,(43a)Hcos⁡(τ1*ω1*-θ)=U,(43b)Hsin(τ1*ω1*-θ)=-V.Equations (43a) and (43b) determine τ1*ω1*-θ uniquely in [0,2π) and hence τ1* uniquely in [θ/ω1*,(θ+2π)/ω1*). To apply the Hopf bifurcation theorem as stated in [27], we state the following lemma.Lemma 2 (see [28]).
Suppose (38) has at least one simple positive root and ω1* is the last such root. Then, iω(τ1*)=iω1* is a simple root of (33) and υ(τ)+iω(τ) is differentiable with respect to τ in a neighbourhood of τ=τ1*.
Next, to establish Hopf bifurcation atτ1=τ1*, we need to verify the transversality condition
(44)dυdτ|τ=τ1*≠0.
Differentiating equations (36a) and (36b) with respect to τ, setting υ=0 and ω=ω1*, solving for dυ/dτ|τ=τ1* and dω/dτ|τ=τ1*, and using (37a) and (37b), we obtain
(45)dυdτ|τ=τ1*=1Γ12+Γ22{ω1*2[4ω1*6+3ω1*4(p12-2p2)+2ω1*2(p22+2p4-2p1p3-q22)+p32-q32-2p2p4+2q2q4]}.
Here,
(46)Γ1=-4ω1*3+2p2ω1*+τ1*(-p1ω1*3+p3ω1*)+q3sin(τ1*ω1*)-2q2ω1*cos⁡(τ1*ω1*),Γ2=-3p1ω1*2+p3+τ1*(ω1*4-p2ω1*2+p4)+(-2q2ω1*)sin(τ1*ω1*)-q3cos⁡(τ1*ω1*),Γ12+Γ22>0,
as iω(τ1*) is a simple root of (33). Let ς=ω1*2; then (38) reduces to ν(ς)=0, where
(47)ν(ς)=ς4+(p12-2p2)ς3+(p22+2p4-2p1p3-q22)ς2+(p32-2p2p4+2q2q4-q32)ς+(p42-q42).
Hence,
(48)dνdς=4ς3+3ς2(p12-2p2)+2ς(p22+2p4-2p1p3-q22)+(p32-q32-2p2p4+2q2q4).
Ifω1*2 is the first positive simple root of (38), then
(49)dνdς|ς=ω1*2>0.Hence, using (45) and (48) we deduce that
(50)dυdτ|τ=τ1*>0.Theorem 3.
Suppose that (38) has at least one simple positive root and ω1* is the last such root. Then, there is a Hopf bifurcation for the system (5) as τ passes upwards through τ1* leading to a periodic solution that bifurcates from E*.Next, we will give the sensible conditions that the Hopf bifurcation occurs around equilibriumE*. Firstly, we need the following important lemma.Definef1=p12-2p2, f2=p22+2p4-2p1p3-q22, f3=p32-2p2p4+2q2q4-q32, f4=p42-q42, and ς=ω1*2, then (38) becomes
(51)ς4+f1ς3+f2ς2+f3ς+f4=0.Lemma 4 (see [28]).
Iff4<0, then the quartic equation
(52)ν(ς)=ς4+f1ς3+f2ς2+f3ς+f4=0
has a strictly positive triple root k1 if and only if(1)
3f12≥8f2;(2)
f1<0 or f2<0;(3)
β1 satisfies the equation 6β2+3f1β+f2=0;(4)
f3=β12(3f1+8β1);(5)
f4=β13(-f1-3β1).We also need the following mild condition.Condition 1.
Either(i)
8f2>3f12;(ii)
f1≥0 and f2≥0;(iii)
or if3f12≥8f2 and also either f1<0 or f2<0, then if β1 is a strictly positive root of the quadratic equation, 6β2+3f1β+f2=0; either f3≠β12(3f1+8β1) or f4≠β13(-f1-3β1).Equation (38) has at least one positive real root for ω1*2 if |p4|<|q4|. By Lemma 2, this is a simple root if Condition 1 is satisfied. Thus, from Lemma 2 and Theorem 3, we can get the following theorem.Theorem 5.
Suppose that(i)
ℛ0>1 and the unique endemic equilibrium E* exists; and(ii)
Condition1 holds and |p4|<|q4| so f4<0.
Then, there is a Hopf bifurcation for the system (5) as τ passes upwards through τ1* leading to a periodic solution that bifurcates from E*.Remark 6.
If (38) has a positive root ω1*, from (37a) and (37b) we can obtain
(53)τj*=1ω1*arcsin[(((q3*4ω1*)2+(q2ω1*2-q4)2)-1(-q3ω1*(ω1*4-p2ω1*2+p4)+(p1ω1*3-p3ω1*)(q2ω1*2-q4))×((q3*4ω1*)2+(q2ω1*2-q4)2)-1)+2jπ((q3*4ω1*)2+(q2ω1*2-q4)2)-1],j=0,1,2,….


## 5. Direction and Stability of the Hopf Bifurcation

In the previous section, we obtain the conditions under which a family of periodic solutions bifurcates from the positive equilibriumE* at the critical value of τ1*. As pointed out in Hassard et al. [29], it is interesting to determine the direction, stability, and period of the periodic solutions bifurcating from the positive equilibrium E*. Following the ideas of Hassard et al., we derive the explicit formulas for determining the properties of the Hopf bifurcation at the critical value of τ1* by using the normal form and the center manifold theory. Throughout this section, we always assume that system (5) undergoes Hopf bifurcation at the positive equilibrium E* for τ=τ1*, and then ±iω1* is corresponding purely imaginary roots of the characteristic equation at the positive equilibrium E*. In the this section, for convenience, we use τ* and ω* instead of τ1* and ω1*, respectively.Letx1(t)=T(t)-T*, x2(t)=I1(t)-I1*, x3(t)=I2(t)-I2*, x4(t)=V1(t)-V1*, x-i(t)=xi(τt), (i=1,2,3,4), and τ=τ*+μ; system (5) is transformed into an functional differential equation (FDE) in ℂ=ℂ([-1,0],ℝ4) as
(54)dxdt=Lμ(xt)+f(μ,xt),
where x(t)=(x1(t),x2(t),x3(t),x4(t))⊤∈ℝ4 and Lμ:ℂ→ℝ,f:ℝ×ℂ→ℝ are given, respectively, by
(55)Lμ(ϕ)=(τ*+μ)×(Ωηa+b0-kT*0-a-b-d000a(1-η)-δ000(1-p)Nδ-c)×(ϕ1(0)ϕ2(0)ϕ3(0)ϕ4(0))+(τ*+μ)(0000kV1*00kT*00000000)(ϕ1(-1)ϕ2(-1)ϕ3(-1)ϕ4(-1)),(56)f(μ,ϕ)=(τ*+μ)(-rTmax⁡ϕ12(0)-kϕ1(0)ϕ4(0)kϕ1(-1)ϕ4(-1)00).
By the Riesz representation theorem, there exists a function η(θ,μ) of bounded variation for θ∈[-1,0], such that
(57)Lμϕ=∫-10‍dη(θ,0)ϕ(θ),
for ϕ∈ℂ.In fact, we can choose(58)η(θ,μ)=(τ*+μ)(Ωηa+b0-kT*0-a-b-d000a(1-η)-δ000(1-p)Nδ-c)δ(θ)-(τ*+μ)(0000kV1*00kT*00000000)δ(θ+1),
where δ is the Dirac delta function. For ϕ∈ℂ′([-1,0],ℝ4), define
(59)A(μ)ϕ={dϕ(θ)dθ,θ∈[-1,0),∫-10‍dη(μ,s)ϕ(s),θ=0,R(μ)ϕ={0,θ∈[-1,0),f(μ,ϕ),θ=0.
Then, system (54) is equivalent to
(60)x˙t=A(μ)xt+R(μ)xt,
where xt(θ)=x(t+θ) for θ∈[-1,0].Forψ∈ℂ1([-1,0],(ℝ4)*), define
(61)A*ψ(s)={-dψ(s)ds,s∈(0,1],∫-10‍dηT(t,0)ψ(-t),s=0,
and a bilinear inner product
(62)〈ψ(s),ϕ(θ)〉=ψ¯(0)ϕ(0)-∫-10‍∫ζ-θθ‍ψ¯(ζ-θ)dη(θ)ϕ(ζ)dζ,
where η(θ)=η(θ,0). Then, A(0) and A* are adjoint operators. By the discussion in Section 4, we know that ±iω*τ* are eigenvalues of A(0). Thus, they are also eigenvalues of A*. We first need to compute the eigenvector of A(0) and A* corresponding to +iω*τ* and -iω*τ*, respectively.Suppose thatq(θ)=(1,a,β,γ)⊤eiω*τ*θ is the eigenvector of A(0) corresponding to +iω*τ*; then A(0)q(θ)=iω*τ*q(θ). It follows from the definition of A(0) and (55), (57), and (58) that(63)τ*(iω*-Ω-ηa-b0kT*-kV1*e-iω*τ*iω*+a+b+d0-kT*e-iω*τ*0-a(1-η)iω*+δ000-(1-p)Nδiω*+c)×q(0)=(0000).Thus, we can easily obtainq(0)=(1,α,β,γ)⊤, where
(64)α=(iω*-Ω-kV1*)e-iω*τ*(ηa+b)e-iω*τ*-(iω*+a+b+d),γ=1kT*[(ηa+b)α-iω*+Ω],β=(iω*+c)γ(1-p)Nδ.Similarly, letq*(s)=D(1,α*,β*,γ*)eiω*τ*s be the eigenvector of A* corresponding to -iω*τ*. By the definition of A* and (55)–(57), we can compute
(65)α*=iω*-ΩkV1*e-iω*τ*,β*=-ηa-b+(iω*+a+b+d)α*a(1-η),γ*=α*kT*e-iω*τ*-kT*iω*+c.In order to assure〈q*(s),q(θ)〉=1, we need to determine the value of D. From (62), we have〈q*(s),q(θ)〉=D¯(1,α*¯,β*¯,γ*¯)(1,α,β,γ)⊤-∫-10‍∫ζ=0θ‍D¯(1,α*¯,β*¯,γ*¯)e-iω*τ*(ζ-θ)dη(θ)(1,α,β,γ)⊤eiω*τ*ζdζ=D¯{∫-10‍(1,α*¯,β*¯,γ*¯)θeiω*τ*θ1+αα*¯+ββ*¯+γγ*¯-∫-10‍(1,α*¯,β*¯,γ*¯)θeiω*τ*θdη(θ)(1,α,β,γ)⊤}=D¯{1+αα*¯+ββ*¯+γγ*¯+τ*e-iω*τ*(α*kV1*+α*kT*γ-)}.Thus, we can choose D as
(67)D=11+αα*¯+ββ*¯+γγ*¯+τ*eiω*τ*(kV1*α*+kT*α*γ-).In the remainder of this section, we use the same notations as in [29]; we first compute the coordinates to describe the center manifold C0 at μ=0. Let xt be the solution of (60) when μ=0. Define
(68)Z(t)=〈q*,xt〉,W(t,θ)=xt(θ)-2Re{z(t)q(θ)}.
On the center manifold C0, we have
(69)W(t,θ)=W(z(t),z-(t),θ),
where
(70)W(z,z-,θ)=W20(θ)z22+W11(θ)zz-+W02(θ)z-22+W30(θ)z36+⋯,z and z- are local coordinates for center manifold C0 in the direction of q* and q*¯. Note that W is real if xt is real. We only consider real solutions. For solution xt∈C0 of (60), since μ=0, we have
(71)z˙(t)=iω*τ*z+q*¯(0)f(0,W(z,z-,θ)+2Re{zq(θ)})=iω*τ*z+q*¯(0)f0(z,z-).
We rewrite this equation as
(72)z˙(t)=iω*τ*z(t)+g(z,z-),
where
(73)g(z,z-)=q*¯(0)f0(z,z-)=g20z22+g11zz-+g02z-22+g21z2z-2+⋯.
It follows from (68) and (70) that
xt(θ)=W(t,θ)-2Re{z(t)q(t)}=W20(θ)z22+W11(θ)zz-+W02(θ)z-22+(1,α,β,γ)Teiω*τ*θz+(1,α-,β-,γ-)Te-iω*τ*θz-+⋯.
It follows together with (56) that
(75)g(z,z-)=q*¯(0)f0(z,z-)=q*¯(0)f(0,xt)=τkD¯(1,α*¯,β*¯,γ*¯)(-rTmax⁡ϕ12(0)-kϕ1(0)ϕ4(0)kϕ1(-1)ϕ4(-1)00)=-τ*D¯rTmax⁡[z+z-+W20(1)(0)z22+W11(1)(0)zz-+W02(1)(0)z-22+o(|(z,z-)3|)]2-τ*D¯k[z+z-+W20(1)(0)z22+W11(1)(0)zz-+W02(1)(0)z-22+o(|(z,z-)3|)]×[γz+γ-z-+W20(4)(0)z22+W11(4)(0)zz-+W02(4)(0)z-22+o(|z,z-|3)]+τ*D¯α-*k[e-iω*τ*z+eiω*τ*z-+W20(1)(-1)z22+W11(1)(-1)zz-+W02(1)(-1)z-22+o(|z,z-|3)]×[γe-iω*τ*z+γ-eiω*τ*z-+W20(4)(-1)z22+W11(4)(-1)zz-+W02(4)(-1)z-22+o(|z,z-|3)].Comparing the coefficients with (73), we have
(76)g20=2τ*D-[-rTmax⁡-kγ+kα-*γe-2iω*τ*],g11=2τ*D-[-rTmax⁡-kReγ+kα-*Reγ],g02=2τ*D-[-rTmax⁡-kγ-+kα-*γ-e2iω*τ*],g21=τ*D-kα-*×[-rTmax⁡(4W11(1)(0)+2W20(1)(0))-k(2W11(4)(0)+W20(4)(0)+γ-W20(1)(0)+2γW11(1)(0))+kα-*(2e-iω*τ*W11(4)(-1)+eiω*τ*W20(4)(-1)+γ-eiω*τ*W20(1)(-1)+2γe-iω*τ*W11(1)(-1))-rTmax⁡].Since there areW20(θ) and W11(θ) in g21, we still need to compute them. From (60) and (68), we have
(77)W˙=x˙t-z˙q-z-˙q-={AW-2Re{q*¯(0)f0q(θ)},θ∈[-1,0),AW-2Re{q*¯(0)f0q(θ)}+f0,θ=0,≜AW+H(z,z-,θ),
where(78)H(z,z-,θ)=H20(θ)z22+H11(θ)zz-+H02(θ)z-22+⋯.
Substituting the corresponding series into (77) and comparing the coefficients, we obtain
(79)(A-2iω*τ*)W20(θ)=-H20,AW11(θ)=-H11,….
From (77), we know that for θ∈[-1,0),
(80)H(z,z-,θ)=-q*¯(0)f0q(θ)-q*(0)f-0q-(θ)=-g(z,z-)q(θ)-g(z,z-)q-(θ).
Comparing the coefficients with (78) gives
(81)H20(θ)=-g20q(θ)-g-02q-(θ),(82)H11(θ)=-g11q(θ)-g-11q-(θ).
From (79), (81), and the definition of A, it follows that
(83)W˙20=2iω*τ*W20(θ)+g20q(θ)+g-02q-(θ).
Notice that q(θ)=(1,α,β,γ)⊤eiω*τ*θ; hence,
(84)W20(θ)=ig20ω*τ*q(0)eiω*τ*θ+ig-023ω*τ*q-(0)e-iω*τ*θ+E1e2iω*τ*θ,
where E1=(E1(1),E1(2),E1(3),E1(4))∈ℝ4 is a constant vector. Similarly, from (79) and (82), we obtain
(85)W11(θ)=-ig11ω*τ*q(0)eiω*τ*θ+ig-11ω*τ*q-(0)e-iω*τ*θ+E2,
where E2=(E2(1),E2(2),E2(3),E2(4))∈ℝ4 is also a constant vector.In what follows, we will seek appropriateE1 and E2. From the definition of A and (79), we obtain
(86)∫-10‍dη(θ)W20(θ)=2iω*τ*W20(θ)-H20(θ),(87)∫-10‍dη(θ)W11(θ)=-H11(θ),
where η(θ)=η(0,θ). By (77), we have
(88)H20(θ)=-g20q(0)-g-20q-(0)+2τ*(-rTmax⁡-kγkγe-2iω*τ*00),(89)H11(θ)=-g11q(0)-g-11q-(0)+2τ*(-rTmax⁡-kReγkReγ00).Substituting (83) and (88) into (86), we obtain
(90)(2ω*τ*I-∫-10‍e2iω*τ*θdη(θ))E1=2τ*(-rTmax⁡-kγkγe-2iω*τ*00),
which leads to(91)(2iω*-Ω-ηa-b0kT*-kV1*e-2iω*τ*2iω*+a+b+d0-kT*e-2iω*τ*0-a(1-η)2iω*+δ000-(1-p)Nδ2iω*+c)×E1=2(-rTmax⁡-kγkγe-2iω*τ*00).It follows that(92)E1(1)=2Δ|-rTmax⁡-kγ-ηa-b0kT*kγe-2iω*τ*2iω*+a+b+d0-kT*e-2iω*τ*0-a(1-η)2iω*+δ000-(1-p)Nδ2iω*+c|,E1(2)=2Δ|2iω*-Ω-rTmax⁡-kγ0kT*-kV1*e-2iω*τ*kγe-2iω*τ*0-kT*e-2iω*τ*0-a(1-η)2iω*+δ000-(1-p)Nδ2iω*+c|,E1(3)=2Δ|2iω*-Ω-ηa-b-rTmax⁡-kγkT*-kV1*e-2iω*τ*2iω*+a+b+dkγe-2iω*τ*-kT*e-2iω*τ*00δ00002iω*+c|,E1(4)=2Δ|2iω*-Ω-ηa-b0-rTmax⁡-kγ-kV1*e-2iω*τ*2iω*+a+b+d0kγe-2iω*τ*0-a(1-η)2iω*+δ000-(1-p)Nδ0|,
where
(93)Δ=|2iω*-Ω-ηa-b0kT*-kV1*e-2iω*τ*2iω*+a+b+d0-kT*e-2iω*τ*0-a(1-η)2iω*+δ000-(1-p)Nδ2iω*+c|.Similarly, substituting (85) and (89) into (87), we can get
(94)(-Ω-ηa-b0kT*-kV1*a+b+d0-kT*0-a(1-η)δ000-(1-p)Nδc)E2=(-rTmax⁡-kReγkReγ00),
and hence(95)E2(1)=2Δ1×|-rTmax⁡-kReγ-ηa-b0kT*kReγa+b+d0-kT*0-a(1-η)δ000-(1-p)Nδc|,E2(2)=2Δ1|-Ω-rTmax⁡-kReγ0kT*-kV1*kReγ0-kT*00δ000-(1-p)Nδc|,E2(3)=2Δ1|-Ω-ηa-b-rTmax⁡-kReγkT*-kV1*a+b+dkReγ-kT*0-a(1-η)00000c|,E2(4)=2Δ1|-Ω-ηa-b0-rTmax⁡-kReγ-kV1*a+b+d0kReγ0-a(1-η)δ000-(1-p)Nδ0|,where
(96)Δ1=|-Ω-ηa-b0kT*-kV1*a+b+d0-kT*0-a(1-η)δ000-(1-p)Nδc|.Thus, we can determineW20(θ) and W11(θ) from (83) and (85). Furthermore, g21 in (75) can be expressed by the parameters and delay. Then, we can compute the following values:
(97)c1(0)=i2ω*τ*(g20g11-2|g11|2-|g02|23)+g212,μ2=-Re{c1(0)}Re{λ′(τ*)},β2=2Re{c1(0)},T2=-Im⁡{c1(0)}+μ2Im⁡{λ′(τ*)}ω*τ*.
By the result of Hassard et al. [29], we have the following.Theorem 7.
In (97), the sign of μ2 determined the direction of Hopf bifurcation: if μ2>0 (μ2<0), then the Hopf bifurcation is supercritical (subcritical) and the bifurcating periodic solution exists for τ>τ* (τ<τ*). β2 determines the stability of the bifurcating periodic solution: the bifurcating periodic solution is stable (unstable) if β2<0 (β2>0), and T2 determines the period of the bifurcating periodic solution: the period increases (decreases) if T2>0 (T2<0).


## 6. Numerical Simulation

In the previous sections, we introduced the analytical tools proposed and used for a qualitative analysis of the system obtaining some results about the dynamics of the system. In this section, we perform a numerical analysis of the model based on the previous results.Our model involves 13 parameters, including the delayτ. In the following, we choose a set of parameters in Table 1. Correspondingly, ℛ0=17.68614487>1 and E*(72.94117647,167.1814866,94.73617573,21315.63954). We can compute that ω*=0.3205862054, and τ*=1.201514430. By Theorem 5, equilibrium E* is locally asymptotically stable when τ<τ* (see Figure 1), and Hopf bifurcation occurs at τ=τ*; a periodic solution exists when τ>τ* (see Figure 2). Furthermore, we compute c1(0)=-13.62874664-2.28357223i. Therefore, Re(c1(0))<0. By Theorem 7, we know that the Hopf bifurcation is supercritical: the bifurcating periodic solutions exist for τ>τ* and they are orbitally asymptotically stable.Table 1
Variables and parameters for viral spread.
Variables and meaning parameters
Values
T
Uninfected CD4+ T-cell population size
1000 mm−3
I
1
Pre-RT
0 mm−3
I
2
Post-RT
0 mm−3
V
1
Infectious virus
10−3 mm−3
s
Source term for uninfected CD4+ T-cell
5 day−1 mm−3
k
Rate at which CD4+ T-cell becomes infected with virus
0.00005 mm3 day−1
μ
1
Death rate of healthy CD4+ T-cell
0.01 day−1
η
Efficacy of RTIs
0.15
a
Transition rate from pre-RT infected CD4+ T cells to productively post-RT
0.4 day−1
b
Reverting rate of infected cells to uninfected class
0.05 day−1
d
Death rate of infected T-cells
0.015 day−1
δ
Death rate of actively infected T-cellsI2
0.6 day−1
N
Number of virions produced by infected CD4+ T-cells
1000 virions cell−1
c
Clearance rate of virus
2.4 day−1
r
Growth rate of T-cells
0.8 day−1
T
max
⁡
Carrying capacity of T-cells
1300 mm3
p
Protease inhibitor efficacy
0.1Whenτ=0.8<τ*, the positive equilibrium E* is stable.
(a)
(b)
(c)
(d)Whenτ=2.4>τ*, the positive equilibrium E* losses its stability and periodic solution occurs.
(a)
(b)
(c)
(d)The ranges of time delayτ are reported in [30, 31], are between 0 and 2 days. By the theory of Hopf bifurcation, we have shown that sustained oscillations are possible in the realistic parameter space. This shows that our model is reasonable.


## 7. Discussion

We have considered a mathematical model for drugs therapy to the infection ofCD4+ T cells in vivo by HIV. The model incorporates the effects of antiretroviral therapy, logistic growth of the CD4+ T cell, and intracellular delay. We have carried out a rigorous mathematical analysis of global dynamics of the model and have shown that the time delay can destabilize the positive equilibrium and lead to periodic solutions through Hopf bifurcation.If we cannot consider the effect of “intracellular” delay, the viral oscillation will not occur [10]. Intracellular delay can induce rich dynamics in the viral system. Moreover, in system (5), we used a logistic term to model the generation and death of target cells. In fact, we can find the logistic term to model the generation by using simulation. And Li and Su have studied that both the “intracellular” delay and target cell can proliferate on virus dynamics [32]. All in all, based on the analytic and simulation results, we can conclude that both the “intracellular” delay and logistic term may give rise to the viral oscillation in the host. Hence, the oscillation behaviors of virus population can be understood in these ways. We will discuss the effect of the “intracellular” delay and logistic term in theory in the future.It is well to know that current treatment regimens cannot eradicate the virus. And the single drug may be highly effective. From the expression of the basic reproductive numberℛ0=T0/T*=(Tmax⁡[r-μ1+(r-μ1)2+4rs/Tmax⁡]/2r)((1-p)(1-η)kNa/c(d+a+b)), we can find that ℛ0 is a decreasing function for p. The value of ℛ0 is smaller for a larger p. That is to say, PIs are positive for the treatment of HIV. Hence, our results show that we need a combination therapy to obtain the better results of drug therapy.Finally, if we assume a constant death ratem for infected but not yet virus-producing cells, the probability of surviving from time t-τ to time t is just e-mτ [14]. Thus the refined model can be rewritten as (4). Hence, we have the following system:
(98)dTdt=s+rT(1-TTmax⁡)-μ1T-kTV1+(ηa+b)I1,dI1dt=ke-mτT(t-τ)V1(t-τ)-(d+a+b)I1,dI2dt=(1-η)aI1-δI2,dV1dt=(1-p)NδI2-cV1,dV2dt=pNδI2-cV2.
It is easy to obtain that the characteristic equation about the positive equilibrium of model (98) is delay dependent coefficients. We can deduce that the stability switches around the positive equilibrium may occur. We leave it in the future.



---

*Source: 294052-2014-01-16.xml*