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.

1. Introduction

Recently there has been a substantial effort in the mathematical modelling of virus dynamics [18]. 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]:

In system (1), the following variables are includes: uninfected cells at time (unit is cells ), infected cells at time (unit is cells ), and virus at time (unit is virions ). Parameters , and are the death rates of the uninfected cells, the infected cells, and the virus particles, respectively. is the contact rate between uninfected cells and the virus particles. is the average number of virus particles produced by an infected 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 ) and post-RT (denoted by ). They assumed that a virus enters a resting 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 in pre-RT class reverts back to uninfected class and the remaining proceeds to post-RT class and becomes productively infected due to presence of RT inhibitors. The model of Srivastava et al. is as follows where is the efficacy of reverse transcriptase inhibitors (RTIs), is the transition rate from pre-RT (i.e., ) infected T cells class to productively post-RT (i.e., ) which is a productively infected class, and 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 where variables and denote infectious and non-infectious virus at time , respectively. And is the total virus concentration at time . Parameter denotes the effectiveness of PIs with 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 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 [1620]. 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 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 T cell, and intracellular delay. Hence, we can obtain the following model: In model (4), , , , , and represent the density of susceptible T cells, infected T cells before reverse transcription (i.e., those infected cells which are in pre-RT class), infected T cells in which reverse transcription is completed (post-RT class), infectious virus, and noninfectious virus at time , respectively. The meaning of the parameters are as follows: is the source term for uninfected T cell, is the rate at which T cell becomes infected with virus, is the death rate of healthy T cell, is the efficacy of RTIs, is the transition rate from pre-RT infected T cells to productively post-RT, is the reverting rate of infected cells to uninfected class, is the death rate of infected cells, is the death rate of actively infected cells , is the number of virions produced by infected T cells, is the clearance rate of virus, is the maximum proliferation rate, is the cell population density at which proliferation shuts off, is the efficacy of protease inhibitor, and is the “intracellular” delay.

Note that the non-infectious HIV virus does not appear in the first four equations of system (4). Thus, we can consider the following subsystem of system (4):

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 with norm where . As usual, the initial condition of (5) is given as where the initial function belongs to the Banach space of continuous functions mapping the initial into . For biological reasons, the initial functions are assumed as 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 to system (5) with initial conditions (8).

Firstly, we present the positivity of the solutions. System (5) can be put into the matrix form

where and is given by

Let be the nonnegative octant in ;  , (where is a function of the variable ) is locally Lipschitz and satisfies the condition where , and .

Due to lemma in [23] any solution of (9) with , say , is such that for all .

System (5) has an uninfected (boundary) equilibrium and an infected (positive) steady state. The uninfected equilibrium is , where

The infected equilibrium is , where

The basic reproductive number is given as . The basic reproductive number 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 ensures the existence of the infected equilibrium .

3. Stability of Uninfected Equilibrium

In this section, we will discuss the stability of the uninfected equilibrium .

Let be any arbitrary equilibrium. To study the stability of the steady state , let us define Then, the linearized system of (5) around the equilibrium is given by where and are matrices given byHence, the characteristic equation of system (5) at is given by where is a identity matrix that is,

Theorem 1. (1) If , is locally asymptotically stable for any time delay . (2) If , is unstable for any time delay . (3) If , it is a critical case.

Proof. For uninfected equilibrium , (18) reduces to where
It is clear that (19) has the characteristic root .
Next, we will consider the transcendental polynomial
For , we have that Obviously, , , and since . We also get This shows that all the roots of (22) have negative real parts for by using Routh-Hurwitz theorem.
In the following, we investigate the existence of purely imaginary roots , , of (21). If and with is a solution of (21), then separating the real and imaginary parts gives Squaring and adding both equations of (24) yields where
Letting yields
If , then . 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 . Therefore, the uninfected equilibrium is locally asymptotically stable for any provided that .
If , the transcendental polynomial (21) becomes It is clear that is a simple root of (28). We further show that any root of (28) must have negative real part except .
In fact, if (28) has imaginary root for some , , and , from (28) we have which, together with , implies that 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 , which implies that the trivial solution of (5) is stable for any time delay .
If , let Note that since and . It follows from the continuity of the function on that equation has at least one positive root. Hence, characteristic equation (19) has at least one positive. Thus, is unstable. Therefore, our results in this theorem are proved.

4. Dynamical Behavior of Endemic Equilibrium

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 . Some conditions for Hopf bifurcation around equilibrium to occur are obtained by using the time delay as a bifurcation parameter.

For endemic equilibrium , (18) reduces tothat is, where

Obviously, . In addition, in view of Routh-Hurwitz criteria, we can easily know that all roots of (33) with have negative real parts if the following condition holds:

Let us consider and assume , where . Substituting and rewriting (33) in terms of its real and imaginary parts, we obtain

Let be such that and ; then (36a) and (36b) reduce toEliminating , we have

Suppose that is the last positive simple root of (38). We will now show that, with this value of , there is a such that and . Given , (37a) and (37b) can be written aswhere where .

Equations determine a unique . With this value of ,Hence,Equations (43a) and (43b) determine uniquely in and hence uniquely in . 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 is the last such root. Then, is a simple root of (33) and is differentiable with respect to in a neighbourhood of .
Next, to establish Hopf bifurcation at , we need to verify the transversality condition Differentiating equations (36a) and (36b) with respect to , setting and , solving for and , and using (37a) and (37b), we obtain Here, as is a simple root of (33). Let ; then (38) reduces to , where Hence,
If is the first positive simple root of (38), then

Hence, using (45) and (48) we deduce that

Theorem 3. Suppose that (38) has at least one simple positive root and is the last such root. Then, there is a Hopf bifurcation for the system (5) as passes upwards through leading to a periodic solution that bifurcates from .

Next, we will give the sensible conditions that the Hopf bifurcation occurs around equilibrium . Firstly, we need the following important lemma.

Define , , , , and , then (38) becomes

Lemma 4 (see [28]). If , then the quartic equation has a strictly positive triple root if and only if(1); (2) or ;(3) satisfies the equation ;(4);(5).

We also need the following mild condition.

Condition 1. Either(i); (ii) and ;(iii)or if and also either or , then if is a strictly positive root of the quadratic equation, ; either or .Equation (38) has at least one positive real root for if . 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) and the unique endemic equilibrium exists; and(ii)Condition 1 holds and so .
Then, there is a Hopf bifurcation for the system (5) as passes upwards through leading to a periodic solution that bifurcates from .

Remark 6. If (38) has a positive root , from (37a) and (37b) we can obtain

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 equilibrium at the critical value of . 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 . 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 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 for , and then is corresponding purely imaginary roots of the characteristic equation at the positive equilibrium . In the this section, for convenience, we use and instead of and , respectively.

Let , , , , , , and ; system (5) is transformed into an functional differential equation (FDE) in as where and , are given, respectively, by By the Riesz representation theorem, there exists a function of bounded variation for , such that for .

In fact, we can choose where is the Dirac delta function. For , define Then, system (54) is equivalent to where for .

For , define and a bilinear inner product where . Then, and are adjoint operators. By the discussion in Section 4, we know that are eigenvalues of . Thus, they are also eigenvalues of . We first need to compute the eigenvector of and corresponding to and , respectively.

Suppose that is the eigenvector of corresponding to ; then . It follows from the definition of and (55), (57), and (58) that

Thus, we can easily obtain , where

Similarly, let be the eigenvector of corresponding to . By the definition of and (55)–(57), we can compute

In order to assure , we need to determine the value of . From (62), we haveThus, we can choose as

In the remainder of this section, we use the same notations as in [29]; we first compute the coordinates to describe the center manifold at . Let be the solution of (60) when . Define On the center manifold , we have where and are local coordinates for center manifold in the direction of and . Note that is real if is real. We only consider real solutions. For solution of (60), since , we have We rewrite this equation as where It follows from (68) and (70) that It follows together with (56) that

Comparing the coefficients with (73), we have

Since there are and in , we still need to compute them. From (60) and (68), we have where Substituting the corresponding series into (77) and comparing the coefficients, we obtain From (77), we know that for , Comparing the coefficients with (78) gives From (79), (81), and the definition of , it follows that Notice that ; hence, where is a constant vector. Similarly, from (79) and (82), we obtain where is also a constant vector.

In what follows, we will seek appropriate and . From the definition of and (79), we obtain where . By (77), we have

Substituting (83) and (88) into (86), we obtain which leads to

It follows that where

Similarly, substituting (85) and (89) into (87), we can get and hencewhere

Thus, we can determine and from (83) and (85). Furthermore, in (75) can be expressed by the parameters and delay. Then, we can compute the following values: By the result of Hassard et al. [29], we have the following.

Theorem 7. In (97), the sign of determined the direction of Hopf bifurcation: if (), then the Hopf bifurcation is supercritical (subcritical) and the bifurcating periodic solution exists for (). determines the stability of the bifurcating periodic solution: the bifurcating periodic solution is stable (unstable) if (), and determines the period of the bifurcating periodic solution: the period increases (decreases) if ().

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, and . We can compute that , and . By Theorem 5, equilibrium is locally asymptotically stable when (see Figure 1), and Hopf bifurcation occurs at ; a periodic solution exists when (see Figure 2). Furthermore, we compute . Therefore, . By Theorem 7, we know that the Hopf bifurcation is supercritical: the bifurcating periodic solutions exist for and they are orbitally asymptotically stable.

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 of T cells in vivo by HIV. The model incorporates the effects of antiretroviral therapy, logistic growth of the 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 , we can find that is a decreasing function for . The value of is smaller for a larger . 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 rate for infected but not yet virus-producing cells, the probability of surviving from time to time is just [14]. Thus the refined model can be rewritten as (4). Hence, we have the following system: 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.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This work is supported by the National Natural Science Foundation of China (nos. 11171284 and 11371306), Basic and Frontier Technology Research Program of Henan Province (nos. 132300410025 and 132300410364), and the Key Project for the Education Department of Henan Province (no. 13A110771).