About this Journal Submit a Manuscript Table of Contents
Abstract and Applied Analysis
Volume 2014 (2014), Article ID 294052, 14 pages
Research Article

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

1College of Mathematics and Information Science, Xinyang Normal University, Xinyang, Henan 464000, China
2College of Forest, Beijing Forest University, Beijing 100083, China

Received 2 October 2013; Accepted 3 November 2013; Published 16 January 2014

Academic Editor: Weiming Wang

Copyright © 2014 Xueyong Zhou and Xiangyun Shi. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


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 by Hence, 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 to that 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 to Eliminating , 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 as where 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 have Thus, 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 </