A Differential Equation Model of HIV Infection of CD𝟒⁺ T-Cells with Delay

Abstract

An epidemic model of HIV infection of CD4⁺ 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 𝑅₀<1. If 𝑅₀<1, the disease-free equilibrium is asymptotically stable and the disease dies out. If 𝑅₀>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.

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)where 𝑃 is an unknown function representing the rate of virus production, 𝑐 is a constant called the clearance rate constant, and 𝑉 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−𝑇𝑇_max), (1.2)where 𝑠 represents the rate at which new T-cells are created from sources within the body, such as the thymus, 𝑑 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 “𝑎” is the maximum proliferation rate and 𝑇_max 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 is CD4⁺ 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 of CD4⁺ T-cells with delay. The ODE model considers a set of cells susceptible to infection, that is, target cells, 𝑇, which, through interactions with virus, 𝑉, 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 of CD4⁺ T-cells with cure rate and delay is studied:𝑑𝑇𝑇𝑑𝑡=𝑠−𝑑𝑇+𝑎𝑇(1−𝑇_max)−𝛽𝑇𝑉+𝜌𝐼,𝑑𝐼𝑑𝑡=𝛽𝑇𝑉−(𝛿+𝜌)𝐼,𝑑𝑉𝑑𝑡=𝑞𝐼−𝑐𝑉, (2.1)where 𝑇 is the number of target cells, 𝐼 is the number of infected cells, and 𝑉 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 𝛽𝑇𝑉 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 𝑉 and 𝑇 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 𝛽𝑇𝑉. In model (2.1), 𝑠 represents the rate at which new T-cells are created from sources, 𝑎 is the maximum proliferation rate of target cells, 𝑇_max is the 𝑇 population density at which proliferation shuts off, 𝑑 is death rate of the T-cells, 𝛿 is the death rate of the infective cells, 𝑞 is the reproductive rate of the infected cells, 𝑐 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 𝑞/𝛿 virions during its lifetime, where the average rate of the virus released by each cell is 𝑞. 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:𝑇(0)>0,𝐼(0)>0,𝑉(0)>0. (2.2)

We denote𝑅³₊={(𝑇,𝐼,𝑉)∈𝑅³,𝑇≥0,𝐼≥0,𝑉≥0}. (2.3)

2.1. Equilibria and the Stability

The non-negative equilibria of (2.1) are 𝐸⁰=(𝑇₀,0,0),𝐸=(𝑇,𝐼,𝑉), (2.4) where 𝑇₀=(𝑇_max/2𝑎)(𝑎−𝑑+(𝑎−𝑑)²+4𝑎𝑠/𝑇_max),𝑇=𝑐(𝛿+𝜌)/𝛽𝑞,𝐼=(1/𝛿)[𝑠−𝑑𝑇+𝑎𝑇(1−𝑇/𝑇_max)], and 𝑉=(𝑞/𝑐)𝐼.

Let 𝑅₀=𝑇₀/𝑇. It is well-known the importance of the value, 𝑅₀, 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. If 𝑅₀<1,𝐸⁰=(𝑇₀,0,0) is locally stable; if 𝑅₀>1,𝐸⁰=(𝑇₀,0,0) is unstable.

Theorem 2.2. There is an 𝑀>0 such that, for any positive solution (𝑇(𝑡),𝐼(𝑡),𝑉(𝑡)) of (2.1), 𝑇(𝑡)≤𝑀,𝐼(𝑡)≤𝑀, and 𝑉(𝑡)≤𝑀, for all large 𝑡.

Theorem 2.3. Suppose that
(i)𝑅₀>1,(ii)(𝑐+𝛿+𝜌+𝑑−𝑎+2𝑎𝑇/𝑇_max)[(−𝑑+𝑎−2𝑎𝑇/𝑇_max)(𝑐+𝛿+𝜌)+𝛽𝑉(𝑐+𝛿)]<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:𝑑𝑇𝑇𝑑𝑡=𝑠−𝑑𝑇+𝑎𝑇(1−𝑇_max)−𝛽𝑇(𝑡−𝜏)𝑉(𝑡−𝜏)+𝜌𝐼,𝑑𝐼𝑑𝑡=𝛽𝑇(𝑡−𝜏)𝑉(𝑡−𝜏)−(𝛿+𝜌)𝐼,𝑑𝑉𝑑𝑡=𝑞𝐼−𝑐𝑉. (3.1)

The time delay is introduced in the system describing the dynamics of the healthy cells. At time 𝑡, only healthy cells that have infected by the virus 𝜏 time units ago (i.e., at time 𝑡−𝜏) become infectious, provided that they have survived the incubation period of 𝜏 units, and given that they were alive at time 𝑡−𝜏 when they infect the healthy cells. Thus, the incidence term of healthy cells is changed from 𝛽𝑇𝑉 to 𝛽𝑇(𝑡−𝜏)𝑉(𝑡−𝜏). However, (3.1) also satisfies the initial conditions: 𝑇(𝜃)=𝑇₀,𝐼(𝜃)=𝐼₀,𝑉(𝜃)=𝑉₀,𝜃∈[−𝜏,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 state 𝐸₀=(𝑇₀,0,0) and an infected state 𝐸=(𝑇,𝐼,𝑉), where, 𝑇,𝐼, and 𝑉 are the same as in Section 2, given by (2.4). Since the uninfected steady state 𝐸₀ is unstable when 𝜏=0 and 𝑅₀<1, incorporation of a delay will not change the instability. Thus, 𝐸₀ is unstable if 𝑅₀>1, which is also the feasibility condition for the infected steady state 𝐸.

We introduce the reproduction number of differential delay model (3.1), which is given by a similar expression: 𝑅₀=𝑇₀/𝑇=𝛽𝑞𝑇₀/𝑐(𝛿+𝜌). 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/𝑐(𝛿+𝜌).

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 equilibrium 𝐸₀ of the differential-delay model (3.1). We consider the local stability in two cases, namely, when 𝑅₀<1, and when 𝑅₀>1.

Theorem 3.1. The disease-free equilibrium 𝐸₀ of (3.1) is locally asymptotically stable if 𝑅₀<1. The disease-free equilibrium is unstable if 𝑅₀>1.

Proof. Linearizing (3.1) around 𝐸₀=(𝑇₀,0,0), we obtain one negative characteristic solution: 𝜆₁=𝑎−𝑑−2𝑎𝑇₀/𝑇_max and the following transcendental characteristic equation for the disease-free equilibrium 𝐸₀ whose solutions (real and complex) give the remaining eigenvalues:𝜆²+(𝛿+𝜌+𝑐)𝜆+𝑐(𝛿+𝜌)−𝑞𝛽𝑇₀𝑒^−𝜆𝜏=0. (3.2)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 𝐸₀ of (3.2) is locally asymptotically stable if 𝑅₀<1 and it is unstable if 𝑅₀>1. To see the claim for the general nonzero delay 𝜏≠0, we first consider the case when 𝑅₀>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)Suppose that 𝜆 is real. Denote the left-hand side of (3.3) as 𝐹(𝜆) and the right-hand side as 𝐺(𝜆). We have that 𝐹(0)=0 and lim_𝜆→∞𝐹(𝜆)=∞. In contrast, the function 𝐺(𝜆) is a decreasing function of 𝜆 and 𝐺(0)=𝑐(𝛿+𝜌)[𝑅₀−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 case 𝑅₀<1. First, we notice that (3.3) does not have non-negative real roots since in this case 𝐹(𝜆) is increasing for 𝜆≥0 while 𝐺(𝜆) is still decreasing function of 𝜆 but 𝐺(0)=𝑐(𝛿+𝜌)[𝑅₀−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 𝜆=𝑖𝑤, and without loss of generality, we may assume that 𝑤>0 is a root of (3.2). That is, the case if and only if 𝑤 satisfies−𝑤²+𝑖(𝛿+𝜌+𝑐)𝑤+𝑐(𝛿+𝜌)−𝑞𝛽𝑇₀cos𝑤𝜏+𝑖𝑞𝛽𝑇₀sin𝑤𝜏=0. (3.4)Separating the real and imaginary parts, we have the following system, satisfied by 𝑤:−𝑤²+𝑐(𝛿+𝜌)=𝑞𝛽𝑇₀cos𝑤𝜏,(𝛿+𝜌+𝑐)𝑤=−𝑞𝛽𝑇₀sin𝑤𝜏. (3.5) 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 𝑤:𝑤⁴+[(𝛿+𝜌+𝑐)²+2𝑐(𝛿+𝜌)]𝑤²+𝑐²(𝛿+𝜌)²−𝑞²𝛽²𝑇²₀=0. (3.6)To reduce this fourth-order equation in to a quadratic equation, we let 𝑧=𝑤² and denote the coefficients as𝑎₁₀=(𝛿+𝜌+𝑐)²𝑎+2𝑐(𝛿+𝜌),₂₀=𝑐²(𝛿+𝜌)²−𝑞²𝛽²𝑇²₀. (3.7)We can rewrite (3.6) as a quadratic equation in 𝑧:𝑧²+𝑎₁₀𝑧+𝑎₂₀=0. (3.8)
Looking back at the coefficients of this quadratic equation, we see that we can expand the square in 𝑎₁₀ while applying the formula for the difference of squares to 𝑎₂₀, we obtain𝑎₁₀=(𝛿+𝜌+𝑐)²𝑎+2𝑐(𝛿+𝜌)>0,₂₀=𝑐²(𝛿+𝜌)²−𝛽²𝑇²₀=𝑐(𝛿+𝜌)[𝑐(𝛿+𝜌)+𝛽𝑇₀][1−𝑅₀]. (3.9)Since 𝑅₀<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 𝑤 such that 𝑖𝑤 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 𝐸₀ is locally asymptotically stable if 𝑅₀<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 𝑅₀>1, that is, the endemic equilibrium 𝐸 exists. To study the stability of the endemic equilibrium 𝐸, we consider the linearization of (3.1) at the point 𝐸. The following transcendental characteristic equation is obtained:𝜆³+𝑎₁𝜆²+𝑎₂𝜆+𝑎₃=𝑒^−𝜆𝜏(𝑏₁𝜆²+𝑏₂𝜆+𝑏₃), (3.10)where the coefficients in this equation are expressed as follows:𝑎₁=𝑐+𝛿+𝑑−𝑎+2𝑎𝑇𝑇_max,𝑎₂=𝑐(𝛿+𝜌)+(𝑐+𝛿+𝜌)(𝑑−𝑎+2𝑎𝑇𝑇_max𝑎−𝜌),₃=𝑐(𝛿+𝜌)(𝑑−𝑎+2𝑎𝑇𝑇_max𝑏−𝜌),₁=−𝛽𝑏𝑉,₂=−𝛽𝑏𝑉(𝑐+𝛿),₃=−𝛽𝑉𝑐𝛿. (3.11)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 𝐸 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 𝑎₁>0,𝑎₂>0, and 𝑎₃>0. On the other hand, 𝑏₁<0,𝑏₂<0, and 𝑏₃<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 𝜆=𝑖𝑤 with 𝑤>0 is a root of (3.14). This is the case if and only if 𝑤 satisfies the following equation:−𝑖𝑤³−𝑎₁𝑤²+𝑖𝑎₂𝑤+𝑎₃=−𝑏₁𝑤²cos𝑤𝜏+𝑖𝑏₂𝑤cos𝑤𝜏+𝑏₃cos𝑤𝜏−𝑖𝑏₁𝑤²sin𝑤𝜏−𝑏₂𝑤sin𝑤𝜏+𝑖𝑏₃sin𝑤𝜏. (3.12)

Separating again the real and imaginary parts, we have the following system that must be satisfied by 𝑤:𝑎₃−𝑎₁𝑤²=(𝑏₃−𝑏₁𝑤²)cos𝑤𝜏−𝑏₂𝑤sin𝑎𝑤𝜏,₂𝑤−𝑤₃=𝑏₂𝑤cos𝑤𝜏+(𝑏₃−𝑏₁𝑤²)sin𝑤𝜏. (3.13)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 𝑤:𝑤⁶+(𝑎²₁−2𝑎₂−𝑏²₁)𝑤⁴+(𝑎²₂−2𝑎₁𝑎₃+2𝑏₁𝑏₃−𝑏²₃)𝑤²+𝑎²₃−𝑏²₃=0. (3.14)Since this equation contains only even powers of 𝑤, we can reduce the order by letting once again 𝑧=𝑤². Then, (3.14) becomes a third-order equation in 𝑧:𝑧³+𝑚₁𝑧²+𝑚₂𝑧+𝑚₃=0, (3.15)where we have used the following notation for the coefficients of (3.15):𝑚₁=𝑎²₁−2𝑎₂−𝑏²₁,𝑚₂=𝑎²₂−2𝑎₁𝑎₃+2𝑏₁𝑏₃−𝑏²₃,𝑚₃=𝑎²₃−𝑏²₃. (3.16)In order to show that the endemic equilibrium 𝐸 is locally stable, we have to show that (3.15) does not have a positive real solution which might give the square of 𝑤, that is, (3.10) cannot have purely imaginary solutions. The lemma below establishes conditions leading to that result.

Lemma 3.2. If 𝑚₁≥0,𝑚₃≥0, and 𝑚₂>0, then (3.15) has no positive real roots.

Proof. We denote the left-hand side of (3.15) as ℎ(𝑧)=𝑧³+𝑚₁𝑧²+𝑚₂𝑧+𝑚₃. We take the derivative of ℎ(𝑧) with respect to 𝑧,ℎ^′(𝑧)=3𝑧²+2𝑚₁𝑧+𝑚₂. We notice that for 𝑧≥0, the derivative ℎ(𝑧)>0, and, therefore, the function ℎ(𝑧) is an increasing function of 𝑧≥0. Since ℎ(0)=𝑚₃>0, it follows that (3.15) has no positive real roots. This completes the proof of the lemma.

Lemma 3.2 implies that there is no 𝑤 such that 𝑖𝑤 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)𝑅₀>1; (ii)𝑚₁≥0,𝑚₃≥0, and 𝑚₂>0. Then the endemic equilibrium 𝐸 of (3.1) is absolutely stable, that is, 𝐸 is asymptotically stable for all values of the delay 𝜏≥0.

Remark 3.4. Theorem 3.3 indicates that if the parameters satisfy conditions (i) and (ii), then the endemic equilibrium 𝐸 of (3.1) is asymptotically stable for all values of the delay, that is, the endemic equilibrium 𝐸 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, if 𝑚₃<0, then we have ℎ(0)<0 and lim_𝑧→∞ℎ(𝑧)=∞. Thus, (3.15) has at least one positive root, say 𝑧₀. Consequently, (3.13) has at least one positive root, denoted by 𝑤₀=√𝑧₀.

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 𝜆(𝜏)=𝜂(𝜏)+𝑖𝑤(𝜏) be the eigenvalue of (3.13) such that for some initial value of the bifurcation parameter 𝜏₀, we have 𝜂(𝜏₀)=0, and 𝑤(𝜏₀)=𝑤₀ (without loss of generality, we may assume 𝑤₀>0). From (3.13), we have𝜏_𝑗=1𝑤₀arccos((𝑎₁𝑏₁−𝑏₂)𝑤⁴₀+(𝑎₂𝑏₂−𝑎₃𝑏₁−𝑎₁𝑏₃)𝑤²₀+𝑎₃𝑏₃𝑏²₂𝑤²₀+(𝑏₃−𝑏₁𝑤²₀)²)+2𝑗𝜋𝑤₀,𝑗=0,1,2,…. (3.17)

Also, we can verify that the following transversal condition:𝑑Re𝜆(𝜏)|𝑑𝜏^𝜏=𝜏0>0 (3.18)holds. By continuity, the real part of 𝜆(𝜏) becomes positive when 𝜏>𝜏₀ and the steady state becomes unstable. Moreover, a Hopf bifurcation occurs when 𝜏 passes through the critical value 𝜏₀ (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 that 𝑤₀ is the largest positive simple root of (3.14). Then, 𝑖𝑤(𝜏₀)=𝑖𝑤₀ is a simple root of (3.10) and 𝜂(𝜏)+𝑖𝑤(𝜏) is differentiable with respect to 𝜏 in a neighborhood of 𝜏=𝜏₀.

After computation, we get that 𝑖𝑤₀ 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]), 𝜂(𝜏)+𝑖𝑤(𝜏) is defined and analytic in a neighborhood of 𝜏=𝜏₀.

Lemma 3.6. Suppose that 𝑥₁,𝑥₂,𝑥₃ are the roots of 𝑔(𝑥)=𝑥³+𝑚₁𝑥²+𝑚₂𝑥+𝑚₃=0(𝑚₂<0), and 𝑥₃ is the largest positive simple root, then 𝑑𝑔(𝑥)|𝑑𝑥^𝑥=𝑥3>0. (3.19)

This proof is omitted.

To establish the Hopf bifurcation at 𝜏=𝜏₀, we need to show that (𝑑Re𝜆(𝜏)/𝑑𝜏)|_𝜏=𝜏0>0. From (3.10) derivation with respect to 𝜏, we get(3𝜆²+2𝑎₁𝜆+𝑎₂)𝑑𝜆𝑑𝜏=[−𝜏𝑒^−𝜆𝜏(𝑏₁𝜆²+𝑏₂𝜆+𝑏₃)+𝑒^−𝜆𝜏(2𝑏₁𝜆+𝑏₂)]𝑑𝜆𝑑𝜏−𝜆𝑒^−𝜆𝜏(𝑏₁𝜆²+𝑏₂𝜆+𝑏₃). (3.20)This gives(𝑑𝜆)𝑑𝜏⁻¹=3𝜆²+2𝑎₁𝜆+𝑎₂+𝜏𝑒^−𝜆𝜏(𝑏₁𝜆²+𝑏₂𝜆+𝑏₃)−𝑒^−𝜆𝜏(2𝑏₁𝜆+𝑏₂)−𝜆𝑒^−𝜆𝜏(𝑏₁𝜆²+𝑏₂𝜆+𝑏₃)=3𝜆²+2𝑎₁𝜆+𝑎₂−𝜆𝑒^−𝜆𝜏(𝑏₁𝜆²+𝑏₂𝜆+𝑏₃)+2𝑏₁𝜆+𝑏₂𝜆(𝑏₁𝜆²+𝑏₂𝜆+𝑏₃)−𝜏𝜆=2𝜆³+𝑎₁𝜆²−𝑎₃−𝜆²(𝜆³+𝑎₁𝜆²+𝑎₂𝜆+𝑎₃)+𝑏₁𝜆²−𝑏₃𝜆²(𝑏₁𝜆²+𝑏₂𝜆+𝑏₃)−𝜏𝜆. (3.21)Thus, Sign{𝑑(Re𝜆)𝑑𝜏}|_{𝜆=𝑖𝑤0}=Sign{Re(𝑑𝜆)𝑑𝜏⁻¹}_{𝜆=𝑖𝑤0}=Sign{Re[2𝜆³+𝑎₁𝜆²−𝑎₃−𝜆²(𝜆³+𝑎₁𝜆²+𝑎₂𝜆+𝑎₃)]_{𝜆=𝑖𝑤0}+Re[𝑏₁𝜆²−𝑏₃𝜆²(𝑏₁𝜆²+𝑏₂𝜆+𝑏₃)]_{𝜆=𝑖𝑤0}}=Sign{Re[−2𝑤³₀𝑖−𝑎₁𝑤²₀−𝑎₃𝑤²₀(−𝑤³₀𝑖−𝑎₁𝑤²₀+𝑎₂𝑤₀𝑖+𝑎₃)]+Re[−𝑏₁𝑤²₀−𝑏₃−𝑤²₀(−𝑏₁𝑤²₀+𝑏₂𝑤₀𝑖+𝑏₃)=]}Sign{2𝑤⁶₀+(𝑎²₁−2𝑎₂)𝑤⁴₀−𝑎²₃𝑤²₀[(𝑎₁𝑤²₀−𝑎₃)²+(𝑤³₀−𝑎₂𝑤₀)²]+𝑏²₃−𝑏²₂𝑤⁴₀𝑤²₀[(𝑏₃−𝑏₁𝑤²₀)²+𝑏²₂𝑤²₀]}=Sign{3𝑤⁴₀+2(𝑎²₁−2𝑎₂−𝑏²₁)𝑤²₀+(𝑎²₂−2𝑎₁𝑎₃+2𝑏₁𝑏₃−𝑏²₂)(𝑎₁𝑤²₀)²+(𝑤³₀−𝑎₂𝑤₀)²}. (3.22)Sinceℎ(𝑧)=𝑧³+𝑚₁𝑧²+𝑚₂𝑧+𝑚₃, (3.23)thus,𝑑ℎ(𝑧)𝑑𝑧=3𝑧²+2(𝑎²₁−2𝑎₂−𝑏²₁)𝑧+(𝑎²₂−2𝑎₁𝑎₃+2𝑏₁𝑏₃−𝑏²₂). (3.24)As 𝑤₀ is the largest positive simple of (3.14), from Lemma 3.6, we have𝑑ℎ(𝑧)|𝑑𝑧^𝑧=𝑤20>0. (3.25)Hence,𝑑Re𝜆|𝑑𝜏^{𝑤=𝑤0,𝜏=𝜏0}=ℎ(𝑤²₀)/𝑑𝑧(𝑎₁𝑤²₀−𝑎₃)²+(𝑤³₀−𝑎₂𝑤₀)²>0. (3.26)The above analysis can be summarized in the following theorem.

Theorem 3.7. Suppose that
(i)𝑅₀>1. If either(ii)𝑚₃<0 or (iii)𝑚₃≥0 and 𝑚₂<0is satisfied, and 𝑤₀ is the largest positive simple root of (3.14), then the endemic equilibrium 𝐸 of the delay model (3.1) is asymptotically stable when 𝜏<𝜏₀ and unstable when 𝜏>𝜏₀, where 𝜏₀=(1/𝑤₀)arccos(((𝑎₁𝑏₁−𝑏₂)𝑤⁴₀+(𝑎₂𝑏₂−𝑎₃𝑏₁−𝑎₁𝑏₃)𝑤²₀+𝑎₃𝑏₃)/(𝑏²₂𝑤²₀+(𝑏₃−𝑏₁𝑤²₀)²)), when 𝜏=𝜏₀, a Hopf bifurcation occurs; that is, a family of periodic solutions bifurcates from 𝐸 as 𝜏 passes through the critical value 𝜏₀.

In this way, using time delay as a bifurcation parameter, Theorem 3.7 indicates that the delay model could exhibit Hopf bifurcation at a certain value 𝜏₀ 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 infected CD4⁺ 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 𝑎 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.

Simulation of the model in this situation shows stable dynamics as presented in Figure 1. 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 𝐸 is asymptotically stable.

(a)

(b)

(c)

Figure 1 

(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.

Next, we use the same set of parameter values as those in Table 1, but we vary the value of “𝑎”: (𝑎=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)

(b)

(c)

(d)

Figure 2 

(a)–(c) are the oscillations of uninfected cells, infected cells, and virus, (d) shows that there is bifurcation.

We also find that the infection would always keep stability when the cure rate 𝜌 is larger. This can be analyzed from the expression of 𝑅₀ 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 𝑎=5 (the value 𝑎 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)

(b)

(c)

Figure 3 

(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.

5. Conclusion

An epidemic model of HIV infection of CD4⁺ 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 𝑅₀<1, the disease-free equilibrium is locally stable and the disease dies out. If 𝑅₀>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(−𝑑𝜏). The third one is that we add the delay on the term of 𝑞𝐼(𝑡−𝜏).

Acknowledgments

The authors thank the anonymous referee for his (or her) valuable comments and suggestions on the previous version of this paper. This work is supported by the National Sciences Foundation of China (10471040), Sciences Exploited Foundation of Shanxi (20081045), and the Foundation of Yuncheng University (20060218).