# A Differential Equation Model of HIV Infection

Authors: Junyuan Yang; Xiaoyan Wang; Fengqin Zhang
Journal: Discrete Dynamics in Nature and Society
 (2008)

---

## Abstract

An epidemic model of HIV infection ofCD4+ T-cells with cure rate and delay is studied. We include a baseline ODE version of the model, and a differential-delay model with a discrete time delay. The ODE model shows that the dynamics is completely determined
by the basic reproduction number R0<1. If
R0<1, the disease-free equilibrium is asymptotically stable and the disease dies out. If R0>1, a unique endemic equilibrium exists and is globally stable in the interior of the feasible region. In the DDE model, the delay stands for the incubation time. We prove the effect of that delay on the stability of the equilibria. We show that the introduction of a time delay in the virus-to-healthy cells transmission term can destabilize the system, and periodic solutions can arise through Hopf bifurcation.


---

## Body

## 1. Introduction

In the last decade, many
mathematical  models  have been developed to describe the immunological response
to infection with human immunodeficiency virus (HIV) (e.g., [1–11], etc.).
These models have been used to explain different phenomena. The models proposed
have principally been linear and nonlinear ordinary differential equation
models, both with and without delay terms. These models focus on the interactions
of susceptible cells, infected cells, viruses, and immune cells. Simple HIV
models have played a significant role in the development of better
understanding of the disease and the various drug therapy strategies used
against it.The simplest HIV dynamic model is(1.1)dVdt=P−cV,where P is an unknown function
representing the rate of virus production, c is a constant called the clearance
rate constant, and V is the virus concentration. The population dynamics of CD4+ T-cells in humans is not well understood.
Nevertheless, a reasonable model for this population of cells is(1.2)T˙=s−dT+aT(1−TTmax),where s represents the rate at
which new T-cells are created from sources within the body, such as the thymus,
d is the death rate per T-cell. T-cells can also be created by proliferation of
existing T-cells. Here, we represent the proliferation by a logistic function
in which “a” is the maximum proliferation rate and Tmax is the T-cell
population density at which proliferation shuts off. The human immune system
can mount a highly specific response against virtually any foreign substance,
even those never seen before in the course of evolution.Like most viruses, HIV is a very simple creature.
Viruses do not have the ability to reproduce independently. Therefore, they
must rely on a host to aid reproduction. Most viruses carry copies of their DNA
and insert this into the host cell`s DNA. Then, when the host cell is
stimulated, it reproduces copies
of the virus. When HIV infects the body, its target isCD4+ T-cells. A protein on the surface of the virus
has a high affinity for the CD4+ protein on the surface of the T-cell. Binding
takes place, and the contents of the HIV are injected into the host T-cell. HIV
differs from most viruses in that it is a retrovirus: it carries a copy of its
RNA which must first be transcribed into DNA. One of the mysteries to the
medical community is why this class of virus has evolved to include this extra
step. After the DNA of the virus has been duplicated by the host cell, it is
reassembled, and new virus particles bud from the surface of the host cell.
This budding can take place slowly, sparing the host cell; or rapidly, bursting
and killing the host cell. The course of infection with HIV is not clearcut.
Clinicians are still arguing about what causes the eventual collapse of the
immune system, resulting in death. What is widely agreed upon, however, is that
there are four main stages of disease progression. First is the initial
innoculum when virus is introduced into the body. Second is the initial
transient—a relatively short period of time when both T-cell population and
virus population are in great flux. This is followed by the third stage,
clinical latency—a period of time when there are extremely large numbers of
virus and T-cells undergoing incredible dynamics, the overall result of which
is an appearance of latency (disease steady state). Finally, there is
AIDS—this is characterized by the T-cells dropping to very low numbers (or
zero) and the virus growing without bound, resulting in death. The transitions
between these four stages are not well understood, and presently, there is
controversy concerning whether the virus directly kills all of T-cells in this
final stage or if there is some other mechanism(s) at work.Current combination anti-retroviral therapies are
widely used to treat HIV. The development of the drugs that are effective against
HIV is a shining example of how to understand the basics of the genesis of HIV
infection, which has led to the rapid development of drugs to combat the
disease; and the principles for the treatment of HIV infection were developed
simultaneously as a result of large, randomized, clinically controlled trials,
and because of the increasing understanding of the dynamics of HIV replication.
Chemotherapy affects the virus once it enters the cell. Through chemotherapy, a
part of infected cells can transform to target cells.As with a single drug, the virus concentration in
plasma fell dramatically for one to two weeks. However, under continued
therapy, after this initial 	 first phase of decline, the
virus continued to fall but at a significantly slower rate. This variation may
have been presented in  previous studies. In the
work of Perelson et al.
[12], the results from
[12, Figure  7.1] show a fast phase followed by what
could be a flat second phase. The reason for this variation among individuals
may lie in the important immunologic component of HIV infection. HIV is thought
to be primarily a noncytopathic virus, and infected cells are lost either
through death, mainly immune-mediated killing, or via  cure, that is, loss
of cccDNA. The second-phase decay has been associated with the
increased rate
of loss of productively infected cells. Antiviral therapy partially blocks the
production of new virious and there is a rapid decline of plasma HIV RNA, but a
vigorous immune response may be needed to drive second-phase decline, which
involves the loss of cells still producing virus. Thus, some process may be
slowing HIV clearance. We show that the pattern of HIV RNA decay can be more
complex than the typical biphasic pattern, with some patients exhibiting
additional phases, raising questions about the need to improve the basic viral
dynamic model. We suggest that including both cytolytic and noncytolytic
mechanisms of infected cell loss will make models more realistic as well as more
accurate.In this paper, we shall investigate an epidemic model
of HIV infection ofCD4+ T-cells with delay. The ODE model considers a
set of cells susceptible to infection, that is, target cells, T, which, through interactions with virus, V, become infected. In addition, infected cells
may also revert to the uninfected state by loss of all cccDNA from their
nucleus at a certain rate per infected cell, which is always omitted in many virus
models, such as Perelson et al. [12]. We extend this model to include a
fixed delay in the system for the infected cells in Section 3. We are
interested in Hopf bifurcation and the presence of sustained oscillations.


## 2. The ODE model

In this section, an epidemic model of HIV infection ofCD4+ T-cells with cure rate and delay is studied:(2.1)dTdt=s−dT+aT(1−TTmax)−βTV+ρI,dIdt=βTV−(δ+ρ)I,dVdt=qI−cV,where T is the number of target cells, I is the number of infected cells, and V is the viral load of the virions. The simplest
and most common method of modelling infection is to augment (2.1) with a
“mass-action” term in which the rate of infection is given by βTV with β being the infection rate constant. This type
of term is sensible since the virus must meet T-cells in order to infect them
and the probability of virus encountering a T-cell at low concentrations (when
V and T motions can be regarded as independent) can be assumed to be
proportional to the product of their density, which is called linear infection rate. Thus, in what
follows, the classical models assume that infected T-cells
are at rate β and the generation of infected T-cells
are at rate βTV. In model (2.1), s represents the rate at which new T-cells are
created from sources, a is the maximum proliferation rate of target
cells, Tmax is the T population density at which
proliferation shuts off, d is death rate of the T-cells, δ is the death rate of the infective cells, q is the
reproductive rate of the infected cells, c is the clearance rate constant of virions, and ρ is the rate of “cure,” that is,
noncytolytic loss of infected cells. Thus, the total rate of disappearance of
infected cells is δ+ρ. The average lifespan of a productively
infected cell is 1/δ. An infected cell produces a total of q/δ virions during its lifetime, where the average rate of the virus released by each cell is q.
Standard and simple arguments show that the solutions of (2.1) exist and stay
positive.However, (2.1) needs to be analyzed  with the following
initial conditions:(2.2)T(0)>0,I(0)>0,V(0)>0.We denote(2.3)R+3={(T,I,V)∈R3,T≥0,I≥0,V≥0}.


### 2.1. Equilibria and the Stability

The non-negative equilibria of (2.1) are (2.4)E0=(T0,0,0),E¯=(T¯,I¯,V¯), where T0=(Tmax/2a)(a−d+(a−d)2+4as/Tmax),T¯=c(δ+ρ)/βq,I¯=(1/δ)[s−dT¯+aT¯(1−T¯/Tmax)], and V¯=(q/c)I¯.LetR0=T0/T¯.
It is well-known the importance of the value,
R0, which is called as the basic reproductive ratio of system (2.1). It represents the average
number of secondary infection caused by a single infected cell
in an entirely susceptible cell population
throughout its infectious period; and it determines the dynamical properties of
(2.1) over along period of time. Based on the result of a differential equation
of HIV infection of CD4+ T-cells with cure rate authored
by Zhou et al. [4], we obtain the following results.Theorem 2.1.
IfR0<1,E0=(T0,0,0) is locally stable; if R0>1,E0=(T0,0,0) is unstable.Theorem 2.2.
There is anM>0 such that, for any positive solution (T(t),I(t),V(t)) of (2.1), T(t)≤M,I(t)≤M, and V(t)≤M, for all large t.Theorem 2.3.
Suppose
that
(i)
R0>1,(ii)
(c+δ+ρ+d−a+2aT¯/Tmax)[(−d+a−2aT¯/Tmax)(c+δ+ρ)+βV¯(c+δ)]<0.
Then, (2.1) is an
orbitally stable periodic orbit.


## 2.1. Equilibria and the Stability

The non-negative equilibria of (2.1) are (2.4)E0=(T0,0,0),E¯=(T¯,I¯,V¯), where T0=(Tmax/2a)(a−d+(a−d)2+4as/Tmax),T¯=c(δ+ρ)/βq,I¯=(1/δ)[s−dT¯+aT¯(1−T¯/Tmax)], and V¯=(q/c)I¯.LetR0=T0/T¯.
It is well-known the importance of the value,
R0, which is called as the basic reproductive ratio of system (2.1). It represents the average
number of secondary infection caused by a single infected cell
in an entirely susceptible cell population
throughout its infectious period; and it determines the dynamical properties of
(2.1) over along period of time. Based on the result of a differential equation
of HIV infection of CD4+ T-cells with cure rate authored
by Zhou et al. [4], we obtain the following results.Theorem 2.1.
IfR0<1,E0=(T0,0,0) is locally stable; if R0>1,E0=(T0,0,0) is unstable.Theorem 2.2.
There is anM>0 such that, for any positive solution (T(t),I(t),V(t)) of (2.1), T(t)≤M,I(t)≤M, and V(t)≤M, for all large t.Theorem 2.3.
Suppose
that
(i)
R0>1,(ii)
(c+δ+ρ+d−a+2aT¯/Tmax)[(−d+a−2aT¯/Tmax)(c+δ+ρ)+βV¯(c+δ)]<0.
Then, (2.1) is an
orbitally stable periodic orbit.


## 3. The Delay Model

In this
section, we introduce a time delay into (2.1) and (2.2) to represent the
incubation time that the vectors need to become infectious. The model for the
CD4+ is exactly as before:(3.1)dTdt=s−dT+aT(1−TTmax)−βT(t−τ)V(t−τ)+ρI,dIdt=βT(t−τ)V(t−τ)−(δ+ρ)I,dVdt=qI−cV.The time delay is introduced in the system describing
the dynamics of the healthy cells. At timet,
only healthy cells that have infected by the virus τ time units ago (i.e., at time t−τ) become infectious, provided that they have
survived the incubation period of τ units, and given that they were alive at time t−τ when they infect the healthy cells. Thus, the
incidence term of healthy cells is changed from βTV to βT(t−τ)V(t−τ).
However, (3.1) also satisfies the initial conditions: T(θ)=T0,I(θ)=I0,V(θ)=V0,θ∈[−τ,0]. All the parameters are the same as in (2.1)
except for the positive constant τ which represents the length of the delay.We find, again, an uninfected steady stateE0=(T0,0,0) and an infected state E¯=(T¯,I¯,V¯), where, T¯,I¯, and V¯ are the same as in Section 2, given by (2.4).
Since the uninfected steady state E0 is unstable when τ=0 and R0<1, incorporation of a delay will not change the
instability. Thus, E0 is unstable if R0>1, which is also the feasibility condition for
the infected steady state E¯.We introduce the reproduction number of differential
delay model (3.1), which is given by a similar expression: R0=T0/T¯=βqT0/c(δ+ρ). Its biological meaning is given as follow, if one virus is introduced in a population of uninfected cells
which infect the total number of secondary infectious during their infectious
period 1/c(δ+ρ).


### 3.1. Local and Global Stability of the Disease-Free Equilibrium

In this section, we turn to study the local and global
stability of the disease-free equilibriumE0 of the differential-delay model (3.1). We
consider the local stability in two cases, namely, when R0<1, and when R0>1.Theorem 3.1.
The
disease-free equilibriumE0 of (3.1) is locally asymptotically stable if R0<1.
The disease-free equilibrium is unstable if R0>1.Proof.
Linearizing (3.1) around E0=(T0,0,0), we obtain one negative characteristic
solution: λ1=a−d−2aT0/Tmax and the following transcendental
characteristic equation for the disease-free equilibrium E0 whose solutions (real and complex) give the
remaining eigenvalues:(3.2)λ2+(δ+ρ+c)λ+c(δ+ρ)−qβT0e−λτ=0.For τ=0, we obtain the same quadratic equation as in
the ODE case. In that case, we know from before that all eigenvalues of the
characteristic (3.2) have negative real part. According to Hurwitz criterion,
when τ=0,
the disease-free equilibrium E0 of (3.2) is locally asymptotically stable if R0<1 and it is unstable if R0>1. To see the claim for the general nonzero delay τ≠0, we first consider the case when R0>1. We expect that in this case, (3.2) has a
positive root and the disease-free equilibrium is unstable. Indeed, to see
this, we rearrange (3.2) in the form(3.3)λ2+(δ+ρ+c)λ=qβT0e−λτ−c(δ+ρ).Suppose that λ is real. Denote the left-hand side of (3.3) as F(λ) and the right-hand side as G(λ). We have that F(0)=0 and limλ→∞F(λ)=∞. In contrast, the function G(λ) is a decreasing function of λ and G(0)=c(δ+ρ)[R0−1]>0. Thus, the two functions must intersect for
some λ>0. Consequently, (3.2) has a positive real
solution and the disease-free equilibrium is unstable.
Now, we turn to the caseR0<1. First, we notice that (3.3) does not have
non-negative real roots since in this case F(λ) is increasing for λ≥0 while G(λ) is still decreasing function of λ but G(0)=c(δ+ρ)[R0−1]<0. Thus, if (3.2) has roots with non-negative
real parts, they must be complex and should have been obtained from a pair of
complex conjugate roots which cross the imaginary axis. Consequently, (3.2)
must have a pair of purely imaginary solutions for some τ>0. Assume that λ=iw, and without loss of generality, we may assume
that w>0 is a root of (3.2). That is, the case if and
only if w satisfies(3.4)−w2+i(δ+ρ+c)w+c(δ+ρ)−qβT0coswτ+iqβT0sinwτ=0.Separating the real and
imaginary parts, we have the following system, satisfied by w:(3.5)−w2+c(δ+ρ)=qβT0coswτ,(δ+ρ+c)w=−qβT0sinwτ. To eliminate the trigonometric functions, we square
both sides of each equation above and we add the squared equations (3.5)  to obtain the following forth-order equation in w:(3.6)w4+[(δ+ρ+c)2+2c(δ+ρ)]w2+c2(δ+ρ)2−q2β2T02=0.To reduce this fourth-order
equation in to a quadratic
equation, we let z=w2 and denote the coefficients as(3.7)a10=(δ+ρ+c)2+2c(δ+ρ),a20=c2(δ+ρ)2−q2β2T02.We can rewrite (3.6) as a
quadratic equation in z:(3.8)z2+a10z+a20=0.
Looking back at the coefficients of this quadratic
equation, we see that we can expand the square ina10 while applying the formula
for the difference of squares to a20, we obtain(3.9)a10=(δ+ρ+c)2+2c(δ+ρ)>0,a20=c2(δ+ρ)2−β2T02=c(δ+ρ)[c(δ+ρ)+βT0][1−R0].Since R0<1, thus, the two roots of (3.8) have positive
product which means that they are complex or they are real but they have the
same sign. In addition, they have negative sum which implies that they are
either real and negative or complex conjugate with negative real parts.
Consequently, (3.8) does not have positive real roots which lead to the
conclusion that there is no w such that iw is a solution of (3.2). Therefore, it follows
from  Rouch's theorem [13] that the real parts of all eigenvalues of the characteristic
equation (3.2) are negative for all values of the delay τ≥0. This implies that E0 is locally asymptotically stable if R0<1. This proves the theorem.


### 3.2. Hopf Bifurcation Analysis

In this section, we determine criteria for Hopf
bifurcation to occur using the time delayτ as the bifurcation parameter. Throughout this
subsection, we will assume that R0>1, that is, the endemic equilibrium E¯ exists. To study the stability of the endemic
equilibrium E¯,
we consider the linearization of (3.1) at the point E¯.
The following transcendental characteristic equation is obtained:(3.10)λ3+a1λ2+a2λ+a3=e−λτ(b1λ2+b2λ+b3),where the coefficients in this
equation are expressed as follows:(3.11)a1=c+δ+d−a+2aT¯Tmax,a2=c(δ+ρ)+(c+δ+ρ)(d−a+2aT¯Tmax−ρ),a3=c(δ+ρ)(d−a+2aT¯Tmax−ρ),b1=−βV¯,b2=−βV¯(c+δ),b3=−βV¯cδ.When τ=0, we obtain the same characteristic equation as
in the ODE case. Consequently, all eigenvalues of the characteristic equation
(3.10) have negative real parts as proved in Theorem 2.1. As a result of Hurwitz
criterion, the endemic equilibrium E¯ of (3.1) is locally asymptotically stable when τ=0. Furthermore, observe again that (3.10) does not
have non-negative real solutions for any τ>0. This implies that a1>0,a2>0,
and a3>0. On the other hand, b1<0,b2<0,
and b3<0. Consequently, the left-hand side in (3.10) is
positive for all τ≥0 while the right-hand side is negative for all τ≥0 and the two cannot be equal for any τ≥0.
We conclude that (3.10) cannot have real non-negative solutions. To rule out
complex conjugate solutions with non-negative real parts, we once again assume
that λ=iw with w>0 is a root of (3.14). This is the case if and
only if w satisfies the following equation:(3.12)−iw3−a1w2+ia2w+a3=−b1w2coswτ+ib2wcoswτ+b3coswτ−ib1w2sinwτ−b2wsinwτ+ib3sinwτ.Separating again the real and imaginary parts, we have
the following system that must be satisfied byw:(3.13)a3−a1w2=(b3−b1w2)coswτ−b2wsinwτ,a2w−w3=b2wcoswτ+(b3−b1w2)sinwτ.We eliminate the trigonometric
functions by squaring both sides of each equation above and adding the
resulting equations. We obtain the following sixth-degree equation for w:(3.14)w6+(a12−2a2−b12)w4+(a22−2a1a3+2b1b3−b32)w2+a32−b32=0.Since this equation contains
only even powers of w,
we can reduce the order by letting once again z=w2.
Then, (3.14) becomes a third-order equation in z:(3.15)z3+m1z2+m2z+m3=0,where we have used the following
notation for the coefficients of (3.15):(3.16)m1=a12−2a2−b12,m2=a22−2a1a3+2b1b3−b32,m3=a32−b32.In order to show that the
endemic equilibrium E¯ is locally stable, we have to show that (3.15)
does not have a positive real solution which might give the square of w,
that is, (3.10) cannot have purely imaginary solutions. The lemma below
establishes conditions leading to that result.Lemma 3.2.
Ifm1≥0,m3≥0,
and m2>0, then (3.15) has no positive real roots.Proof.
We denote the left-hand side of (3.15) as h(z)=z3+m1z2+m2z+m3. We take the derivative of h(z) with respect to z,h′(z)=3z2+2m1z+m2. We notice that for z≥0,
the derivative h(z)>0, and, therefore, the function h(z) is an increasing function of z≥0. Since h(0)=m3>0, it follows that (3.15) has no positive real
roots. This completes the proof of the lemma.Lemma3.2 implies that there is no w such that iw is an eigenvalue of the characteristic (3.10).
Therefore, by Rouche's theorem [13, Theorem 9.17.4], the real parts of all
eigenvalues of (3.10) are negative for all values of the delay τ≥0. Summarizing the above analysis, we have the
following theorem.Theorem 3.3.
Assume that
(i)
R0>1;(ii)
m1≥0,m3≥0, and m2>0.
Then the endemic equilibrium E¯ of (3.1) is absolutely stable, that is, E¯ is asymptotically stable for all values of the
delay τ≥0.Remark 3.4.
Theorem3.3 indicates that if the parameters satisfy
conditions (i) and (ii), then the endemic equilibrium E¯ of (3.1) is asymptotically stable for all
values of the delay, that is, the endemic equilibrium E¯ of (3.1) is asymptotically stable independent
of the delay. However, we should point out that if the conditions in Theorem
3.3, particularly any of the inequalities in (ii), are not satisfied, then the
stability of the endemic equilibrium depends on the delay value and as the
delay varies, the endemic equilibrium can lose stability which can lead to
oscillations.For example, ifm3<0, then we have h(0)<0 and limz→∞h(z)=∞. Thus, (3.15) has at least one positive root,
say z0. Consequently, (3.13) has at least one positive
root, denoted by w0=z0.Now, we turn to the bifurcation analysis. We use the
delayτ as bifurcation parameter. We view the
solutions of (3.10) as functions of the bifurcation parameter τ.
Let λ(τ)=η(τ)+iw(τ) be the eigenvalue of (3.13) such that for some
initial value of the bifurcation parameter τ0,
we have η(τ0)=0, and w(τ0)=w0 (without loss of generality, we may assume w0>0). From (3.13), we have(3.17)τj=1w0arccos((a1b1−b2)w04+(a2b2−a3b1−a1b3)w02+a3b3b22w02+(b3−b1w02)2)+2jπw0,j=0,1,2,….Also, we can verify that the following transversal
condition:(3.18)dReλ(τ)dτ|τ=τ0>0holds. By continuity, the real
part of λ(τ) becomes positive when τ>τ0 and the steady state becomes unstable.
Moreover, a Hopf bifurcation occurs when τ passes through the critical value τ0 (see [14]).To apply  the Hopf bifurcation
theorem as stated in
Marsden and McCracken [15], we state and prove the following
theorem:Theorem 3.5.
Suppose
thatw0 is the largest positive simple root of (3.14).
Then, iw(τ0)=iw0 is a simple root of (3.10) and η(τ)+iw(τ) is differentiable with respect to τ in a neighborhood of τ=τ0.After computation, we get thatiw0 is a simple root of (3.10), which is an
analytic equation, and so, using the analytic version of the implicit function
theorem (Chow and Hale [16]), η(τ)+iw(τ) is defined and analytic in a neighborhood of τ=τ0.Lemma 3.6.
Suppose
thatx1,x2,x3 are the roots of g(x)=x3+m1x2+m2x+m3=0(m2<0), and x3 is the largest positive simple root, then (3.19)dg(x)dx|x=x3>0.This proof is omitted.To establish the Hopf bifurcation atτ=τ0, we need to show that (dReλ(τ)/dτ)|τ=τ0>0. From (3.10)
derivation with respect to τ, we get(3.20)(3λ2+2a1λ+a2)dλdτ=[−τe−λτ(b1λ2+b2λ+b3)+e−λτ(2b1λ+b2)]dλdτ−λe−λτ(b1λ2+b2λ+b3).This gives(3.21)(dλdτ)−1=3λ2+2a1λ+a2+τe−λτ(b1λ2+b2λ+b3)−e−λτ(2b1λ+b2)−λe−λτ(b1λ2+b2λ+b3)=3λ2+2a1λ+a2−λe−λτ(b1λ2+b2λ+b3)+2b1λ+b2λ(b1λ2+b2λ+b3)−τλ=2λ3+a1λ2−a3−λ2(λ3+a1λ2+a2λ+a3)+b1λ2−b3λ2(b1λ2+b2λ+b3)−τλ.Thus,
(3.22)Sign{d(Reλ)dτ}|λ=iw0=Sign{Re(dλdτ)−1}λ=iw0=Sign{Re[2λ3+a1λ2−a3−λ2(λ3+a1λ2+a2λ+a3)]λ=iw0+Re[b1λ2−b3λ2(b1λ2+b2λ+b3)]λ=iw0}=Sign{Re[−2w03i−a1w02−a3w02(−w03i−a1w02+a2w0i+a3)]+Re[−b1w02−b3−w02(−b1w02+b2w0i+b3)]}=Sign{2w06+(a12−2a2)w04−a32w02[(a1w02−a3)2+(w03−a2w0)2]+b32−b22w04w02[(b3−b1w02)2+b22w02]}=Sign{3w04+2(a12−2a2−b12)w02+(a22−2a1a3+2b1b3−b22)(a1w02)2+(w03−a2w0)2}.Since(3.23)h(z)=z3+m1z2+m2z+m3,thus,(3.24)dh(z)dz=3z2+2(a12−2a2−b12)z+(a22−2a1a3+2b1b3−b22).As w0 is the largest positive simple of (3.14), from
Lemma 3.6, we have(3.25)dh(z)dz|z=w02>0.Hence,(3.26)dReλdτ|w=w0,τ=τ0=h(w02)/dz(a1w02−a3)2+(w03−a2w0)2>0.The above analysis can be
summarized in the following theorem.Theorem 3.7.
Suppose
that
(i)
R0>1.
If either(ii)
m3<0   or(iii)
m3≥0 and m2<0is satisfied, and w0 is the largest positive simple root of (3.14),
then the endemic equilibrium E¯ of the delay model (3.1) is asymptotically
stable when τ<τ0 and unstable when τ>τ0, where τ0=(1/w0)arccos(((a1b1−b2)w04+(a2b2−a3b1−a1b3)w02+a3b3)/(b22w02+(b3−b1w02)2)),  when τ=τ0, a Hopf bifurcation occurs; that is, a family
of periodic solutions bifurcates from E¯ as τ passes through the critical value τ0.In this way, using time delay as a bifurcation
parameter, Theorem3.7 indicates that the delay model could exhibit Hopf
bifurcation at a certain value τ0 of the delay if the parameters satisfy
conditions (ii) or (iii). They show that the introduction of a time delay in
the virus-to-uninfected cells transmission term can destabilize the system and
periodic solutions can arise through Hopf bifurcation.


## 3.1. Local and Global Stability of the Disease-Free Equilibrium

In this section, we turn to study the local and global
stability of the disease-free equilibriumE0 of the differential-delay model (3.1). We
consider the local stability in two cases, namely, when R0<1, and when R0>1.Theorem 3.1.
The
disease-free equilibriumE0 of (3.1) is locally asymptotically stable if R0<1.
The disease-free equilibrium is unstable if R0>1.Proof.
Linearizing (3.1) around E0=(T0,0,0), we obtain one negative characteristic
solution: λ1=a−d−2aT0/Tmax and the following transcendental
characteristic equation for the disease-free equilibrium E0 whose solutions (real and complex) give the
remaining eigenvalues:(3.2)λ2+(δ+ρ+c)λ+c(δ+ρ)−qβT0e−λτ=0.For τ=0, we obtain the same quadratic equation as in
the ODE case. In that case, we know from before that all eigenvalues of the
characteristic (3.2) have negative real part. According to Hurwitz criterion,
when τ=0,
the disease-free equilibrium E0 of (3.2) is locally asymptotically stable if R0<1 and it is unstable if R0>1. To see the claim for the general nonzero delay τ≠0, we first consider the case when R0>1. We expect that in this case, (3.2) has a
positive root and the disease-free equilibrium is unstable. Indeed, to see
this, we rearrange (3.2) in the form(3.3)λ2+(δ+ρ+c)λ=qβT0e−λτ−c(δ+ρ).Suppose that λ is real. Denote the left-hand side of (3.3) as F(λ) and the right-hand side as G(λ). We have that F(0)=0 and limλ→∞F(λ)=∞. In contrast, the function G(λ) is a decreasing function of λ and G(0)=c(δ+ρ)[R0−1]>0. Thus, the two functions must intersect for
some λ>0. Consequently, (3.2) has a positive real
solution and the disease-free equilibrium is unstable.
Now, we turn to the caseR0<1. First, we notice that (3.3) does not have
non-negative real roots since in this case F(λ) is increasing for λ≥0 while G(λ) is still decreasing function of λ but G(0)=c(δ+ρ)[R0−1]<0. Thus, if (3.2) has roots with non-negative
real parts, they must be complex and should have been obtained from a pair of
complex conjugate roots which cross the imaginary axis. Consequently, (3.2)
must have a pair of purely imaginary solutions for some τ>0. Assume that λ=iw, and without loss of generality, we may assume
that w>0 is a root of (3.2). That is, the case if and
only if w satisfies(3.4)−w2+i(δ+ρ+c)w+c(δ+ρ)−qβT0coswτ+iqβT0sinwτ=0.Separating the real and
imaginary parts, we have the following system, satisfied by w:(3.5)−w2+c(δ+ρ)=qβT0coswτ,(δ+ρ+c)w=−qβT0sinwτ. To eliminate the trigonometric functions, we square
both sides of each equation above and we add the squared equations (3.5)  to obtain the following forth-order equation in w:(3.6)w4+[(δ+ρ+c)2+2c(δ+ρ)]w2+c2(δ+ρ)2−q2β2T02=0.To reduce this fourth-order
equation in to a quadratic
equation, we let z=w2 and denote the coefficients as(3.7)a10=(δ+ρ+c)2+2c(δ+ρ),a20=c2(δ+ρ)2−q2β2T02.We can rewrite (3.6) as a
quadratic equation in z:(3.8)z2+a10z+a20=0.
Looking back at the coefficients of this quadratic
equation, we see that we can expand the square ina10 while applying the formula
for the difference of squares to a20, we obtain(3.9)a10=(δ+ρ+c)2+2c(δ+ρ)>0,a20=c2(δ+ρ)2−β2T02=c(δ+ρ)[c(δ+ρ)+βT0][1−R0].Since R0<1, thus, the two roots of (3.8) have positive
product which means that they are complex or they are real but they have the
same sign. In addition, they have negative sum which implies that they are
either real and negative or complex conjugate with negative real parts.
Consequently, (3.8) does not have positive real roots which lead to the
conclusion that there is no w such that iw is a solution of (3.2). Therefore, it follows
from  Rouch's theorem [13] that the real parts of all eigenvalues of the characteristic
equation (3.2) are negative for all values of the delay τ≥0. This implies that E0 is locally asymptotically stable if R0<1. This proves the theorem.


## 3.2. Hopf Bifurcation Analysis

In this section, we determine criteria for Hopf
bifurcation to occur using the time delayτ as the bifurcation parameter. Throughout this
subsection, we will assume that R0>1, that is, the endemic equilibrium E¯ exists. To study the stability of the endemic
equilibrium E¯,
we consider the linearization of (3.1) at the point E¯.
The following transcendental characteristic equation is obtained:(3.10)λ3+a1λ2+a2λ+a3=e−λτ(b1λ2+b2λ+b3),where the coefficients in this
equation are expressed as follows:(3.11)a1=c+δ+d−a+2aT¯Tmax,a2=c(δ+ρ)+(c+δ+ρ)(d−a+2aT¯Tmax−ρ),a3=c(δ+ρ)(d−a+2aT¯Tmax−ρ),b1=−βV¯,b2=−βV¯(c+δ),b3=−βV¯cδ.When τ=0, we obtain the same characteristic equation as
in the ODE case. Consequently, all eigenvalues of the characteristic equation
(3.10) have negative real parts as proved in Theorem 2.1. As a result of Hurwitz
criterion, the endemic equilibrium E¯ of (3.1) is locally asymptotically stable when τ=0. Furthermore, observe again that (3.10) does not
have non-negative real solutions for any τ>0. This implies that a1>0,a2>0,
and a3>0. On the other hand, b1<0,b2<0,
and b3<0. Consequently, the left-hand side in (3.10) is
positive for all τ≥0 while the right-hand side is negative for all τ≥0 and the two cannot be equal for any τ≥0.
We conclude that (3.10) cannot have real non-negative solutions. To rule out
complex conjugate solutions with non-negative real parts, we once again assume
that λ=iw with w>0 is a root of (3.14). This is the case if and
only if w satisfies the following equation:(3.12)−iw3−a1w2+ia2w+a3=−b1w2coswτ+ib2wcoswτ+b3coswτ−ib1w2sinwτ−b2wsinwτ+ib3sinwτ.Separating again the real and imaginary parts, we have
the following system that must be satisfied byw:(3.13)a3−a1w2=(b3−b1w2)coswτ−b2wsinwτ,a2w−w3=b2wcoswτ+(b3−b1w2)sinwτ.We eliminate the trigonometric
functions by squaring both sides of each equation above and adding the
resulting equations. We obtain the following sixth-degree equation for w:(3.14)w6+(a12−2a2−b12)w4+(a22−2a1a3+2b1b3−b32)w2+a32−b32=0.Since this equation contains
only even powers of w,
we can reduce the order by letting once again z=w2.
Then, (3.14) becomes a third-order equation in z:(3.15)z3+m1z2+m2z+m3=0,where we have used the following
notation for the coefficients of (3.15):(3.16)m1=a12−2a2−b12,m2=a22−2a1a3+2b1b3−b32,m3=a32−b32.In order to show that the
endemic equilibrium E¯ is locally stable, we have to show that (3.15)
does not have a positive real solution which might give the square of w,
that is, (3.10) cannot have purely imaginary solutions. The lemma below
establishes conditions leading to that result.Lemma 3.2.
Ifm1≥0,m3≥0,
and m2>0, then (3.15) has no positive real roots.Proof.
We denote the left-hand side of (3.15) as h(z)=z3+m1z2+m2z+m3. We take the derivative of h(z) with respect to z,h′(z)=3z2+2m1z+m2. We notice that for z≥0,
the derivative h(z)>0, and, therefore, the function h(z) is an increasing function of z≥0. Since h(0)=m3>0, it follows that (3.15) has no positive real
roots. This completes the proof of the lemma.Lemma3.2 implies that there is no w such that iw is an eigenvalue of the characteristic (3.10).
Therefore, by Rouche's theorem [13, Theorem 9.17.4], the real parts of all
eigenvalues of (3.10) are negative for all values of the delay τ≥0. Summarizing the above analysis, we have the
following theorem.Theorem 3.3.
Assume that
(i)
R0>1;(ii)
m1≥0,m3≥0, and m2>0.
Then the endemic equilibrium E¯ of (3.1) is absolutely stable, that is, E¯ is asymptotically stable for all values of the
delay τ≥0.Remark 3.4.
Theorem3.3 indicates that if the parameters satisfy
conditions (i) and (ii), then the endemic equilibrium E¯ of (3.1) is asymptotically stable for all
values of the delay, that is, the endemic equilibrium E¯ of (3.1) is asymptotically stable independent
of the delay. However, we should point out that if the conditions in Theorem
3.3, particularly any of the inequalities in (ii), are not satisfied, then the
stability of the endemic equilibrium depends on the delay value and as the
delay varies, the endemic equilibrium can lose stability which can lead to
oscillations.For example, ifm3<0, then we have h(0)<0 and limz→∞h(z)=∞. Thus, (3.15) has at least one positive root,
say z0. Consequently, (3.13) has at least one positive
root, denoted by w0=z0.Now, we turn to the bifurcation analysis. We use the
delayτ as bifurcation parameter. We view the
solutions of (3.10) as functions of the bifurcation parameter τ.
Let λ(τ)=η(τ)+iw(τ) be the eigenvalue of (3.13) such that for some
initial value of the bifurcation parameter τ0,
we have η(τ0)=0, and w(τ0)=w0 (without loss of generality, we may assume w0>0). From (3.13), we have(3.17)τj=1w0arccos((a1b1−b2)w04+(a2b2−a3b1−a1b3)w02+a3b3b22w02+(b3−b1w02)2)+2jπw0,j=0,1,2,….Also, we can verify that the following transversal
condition:(3.18)dReλ(τ)dτ|τ=τ0>0holds. By continuity, the real
part of λ(τ) becomes positive when τ>τ0 and the steady state becomes unstable.
Moreover, a Hopf bifurcation occurs when τ passes through the critical value τ0 (see [14]).To apply  the Hopf bifurcation
theorem as stated in
Marsden and McCracken [15], we state and prove the following
theorem:Theorem 3.5.
Suppose
thatw0 is the largest positive simple root of (3.14).
Then, iw(τ0)=iw0 is a simple root of (3.10) and η(τ)+iw(τ) is differentiable with respect to τ in a neighborhood of τ=τ0.After computation, we get thatiw0 is a simple root of (3.10), which is an
analytic equation, and so, using the analytic version of the implicit function
theorem (Chow and Hale [16]), η(τ)+iw(τ) is defined and analytic in a neighborhood of τ=τ0.Lemma 3.6.
Suppose
thatx1,x2,x3 are the roots of g(x)=x3+m1x2+m2x+m3=0(m2<0), and x3 is the largest positive simple root, then (3.19)dg(x)dx|x=x3>0.This proof is omitted.To establish the Hopf bifurcation atτ=τ0, we need to show that (dReλ(τ)/dτ)|τ=τ0>0. From (3.10)
derivation with respect to τ, we get(3.20)(3λ2+2a1λ+a2)dλdτ=[−τe−λτ(b1λ2+b2λ+b3)+e−λτ(2b1λ+b2)]dλdτ−λe−λτ(b1λ2+b2λ+b3).This gives(3.21)(dλdτ)−1=3λ2+2a1λ+a2+τe−λτ(b1λ2+b2λ+b3)−e−λτ(2b1λ+b2)−λe−λτ(b1λ2+b2λ+b3)=3λ2+2a1λ+a2−λe−λτ(b1λ2+b2λ+b3)+2b1λ+b2λ(b1λ2+b2λ+b3)−τλ=2λ3+a1λ2−a3−λ2(λ3+a1λ2+a2λ+a3)+b1λ2−b3λ2(b1λ2+b2λ+b3)−τλ.Thus,
(3.22)Sign{d(Reλ)dτ}|λ=iw0=Sign{Re(dλdτ)−1}λ=iw0=Sign{Re[2λ3+a1λ2−a3−λ2(λ3+a1λ2+a2λ+a3)]λ=iw0+Re[b1λ2−b3λ2(b1λ2+b2λ+b3)]λ=iw0}=Sign{Re[−2w03i−a1w02−a3w02(−w03i−a1w02+a2w0i+a3)]+Re[−b1w02−b3−w02(−b1w02+b2w0i+b3)]}=Sign{2w06+(a12−2a2)w04−a32w02[(a1w02−a3)2+(w03−a2w0)2]+b32−b22w04w02[(b3−b1w02)2+b22w02]}=Sign{3w04+2(a12−2a2−b12)w02+(a22−2a1a3+2b1b3−b22)(a1w02)2+(w03−a2w0)2}.Since(3.23)h(z)=z3+m1z2+m2z+m3,thus,(3.24)dh(z)dz=3z2+2(a12−2a2−b12)z+(a22−2a1a3+2b1b3−b22).As w0 is the largest positive simple of (3.14), from
Lemma 3.6, we have(3.25)dh(z)dz|z=w02>0.Hence,(3.26)dReλdτ|w=w0,τ=τ0=h(w02)/dz(a1w02−a3)2+(w03−a2w0)2>0.The above analysis can be
summarized in the following theorem.Theorem 3.7.
Suppose
that
(i)
R0>1.
If either(ii)
m3<0   or(iii)
m3≥0 and m2<0is satisfied, and w0 is the largest positive simple root of (3.14),
then the endemic equilibrium E¯ of the delay model (3.1) is asymptotically
stable when τ<τ0 and unstable when τ>τ0, where τ0=(1/w0)arccos(((a1b1−b2)w04+(a2b2−a3b1−a1b3)w02+a3b3)/(b22w02+(b3−b1w02)2)),  when τ=τ0, a Hopf bifurcation occurs; that is, a family
of periodic solutions bifurcates from E¯ as τ passes through the critical value τ0.In this way, using time delay as a bifurcation
parameter, Theorem3.7 indicates that the delay model could exhibit Hopf
bifurcation at a certain value τ0 of the delay if the parameters satisfy
conditions (ii) or (iii). They show that the introduction of a time delay in
the virus-to-uninfected cells transmission term can destabilize the system and
periodic solutions can arise through Hopf bifurcation.


## 4. Simulation

In the previous
sections, we introduced the analytical tools proposed and used them 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.Clinical data are becoming more available, making it
possible to get actual values (or orders of values) directly for the individual
parameters in the model. By this, it is meant that it is possible to calculate
the actual rates for the different processes described above based on data
collected from clinical experiments. For example, it has been shown that
infectedCD4+ T cells live less than 1-2 days; therefore, we choose the rate of loss of
infected T-cells, δ,
to values between 0.2 and 1.0. When this type of information is not
available, estimation of the parameters can be determined from simulations
through behavior studies. Periodic solution and sensitivity analyzes can be
carried out for each parameter to get a good understanding of the different
behaviors seen for variations of these values. For example, the parameter a in the model (representing the maximum
proliferation rate of target cells) is not verifiable clinically; however,
since it is a bifurcation parameter, we know that for small values, the
infection would die out and that for large values, the infection persists. This
may be an indication to clinicians that finding a drug which lowers this viral
production may aid in suppressing the disease. In general, this process can be
helpful to clinicians, as a range for possible parameter values can be
suggested. A complete list of parameters and their estimated values for this
model is given in Table 1.Table 1
Variables and parameters for viral spread.
Parameters and variablesValuesDependent variablesTUninfectedCD4+  T-cell population size50mm−3IInfectedCD4+ T-cell density80VInitial density of HIV RNA100mm−3Parameters and constantssSource term for uninfectedCD4+ T-cells5(day)−1(mm−3)dNatural death rate ofCD4+ T-cells0.01day−1aGrowth rate ofCD4+ T-cell population0.8day−1TmaxMaximal population level ofCD4+ T-cells1200mm3day−1βRateCD4+ T-cells become infected with virus0.00024mm−3ρRate of cure0.01day−1δBlanket death rate of infectedCD4+ T-cells0.4day−1qReproductively rate of the infectedCD4+ T-cells1000mm3day−1cDeath rate of free virus8day−1Simulation of the model in this situation
shows stable dynamics as presented in Figure1.
Figures 1(a)–1(c) show that uninfected cells, infected cells, and virus converge
to their equilibrium with the parametric values as stated in Table 1. They show
that the equilibrium E¯ is asymptotically stable.(a)-(c) show that uninfected
cells, infected cells and virus converge to their equilibrium with parametric values as stated in
the text withτ=0.4. They show that the equilibrium is asymptotically stable.
(a)(b)(c)Next, we use the same set of parameter values as those
in Table1, but we vary the value of “a”: (a=5). Thus, the conditions of Theorem 3.7 are
satisfied. Figures 2(a)–2(c) are the
oscillations  of uninfected cells, infected cells, and virus.
Figure 2(d)
shows that there is a periodic solution.(a)–(c) are the oscillations of uninfected
cells, infected cells, and virus, (d) shows that there is bifurcation.
(a)(b)(c)(d)We also find that the infection would always keep
stability when the cure rateρ is larger. This can be analyzed from the
expression of R0 and  the conditions of Theorem 3.7 For
example, we know the oscillations of uninfected cells, infected cells, and
virus in Figure 3; and if we select ρ=0.3,ρ=0.3,
and a=5 (the value a is same as in Figure 3) and the other parameter
values are same in Table 1, then the infection would be
differential stabilities (see Figure 3). Thus, we can claim that the
cure rate ρ is a very important parameter. The results
show that if we improve the cure rate, we will control the disease.(a)–(c) show the uninfected cells, infected cells, and virus withρ=0.01 and ρ=0.3. They show that the cure rate is an important
parameter.
(a)(b)(c)


## 5. Conclusion

An epidemic model of HIV infection ofCD4+ T-cells with cure rate and delay is studied.
Mathematical analyzes of the model equations with regard to invariance of
non-negativity, boundedness of solutions, nature of equilibria, as well as
permanence and global stability are analyzed. The basic reproduction number is
obtained and it completely determines the dynamics of the ODE model. If R0<1, the disease-free equilibrium is locally stable
and the disease dies out. If R0>1, a unique endemic equilibrium exists and
be absolutely  stable. We determine criteria for
Hopf bifurcation using the time delay as the bifurcation parameter based on the
differential-delay model. They show that positive equilibrium is locally
asymptotically stable when time delay is suitably small, while a loss of
stability by a Hopf bifurcation can occur as the delay increases. Hopf
bifurcation has helped us in finding the existence of a region of instability
in the neighborhood of a nonzero endemic equilibrium where the population will
survive undergoing regular fluctuations.There is still some work to do for this model. The
first one is knowing under what condition the
disease equilibrium be globally stable. The second
one is that we want to know that the mass action law is standard incidence or
other general interaction, the survival probability is exp(−dτ). The third one is that we add the delay on
the term of qI(t−τ).



---

*Source: 903678-2009-02-16.xml*