A DELAY DYNAMIC MODEL FOR HIV INFECTED IMMUNE RESPONSE

Abstract

Human Immune Deficiency Virus (or simply HIV) induces a persistent infection that leads to AIDS causing death in almost every infected individual. As HIV affects the immune system directly by attacking the CD4+ T cells, to exterminate the infection, the natural immune system produces virus-specific cytotoxic T lymphocytes(CTLs) that kills the infected CD4+ T cells. The reduced CD4+ T cell count produce reduced amount of cytokines to stimulate the production of CTLs to fight the invaders that weakens the body immunity succeeding to AIDS. In this paper, we introduce a mathematical model with discrete time-delay to represent this cell dynamics between CD4+ T cells and the CTLs under HIV infection. A modified functional form has been considered to describe the infection mechanism. Characteristics of the system are studied through mathematical analysis. Numerical simulations are carried out to illustrate the analytical findings.

1. Introduction

CD4+ T cells or Helper T cells play the key role in human immune system. When human body recognizes foreign antigen, it is delivered to the lymph system, where the virious and infected cells are ingested by macrophage. This macrophage, which tags itself for destruction, displays the antigens on its own exterior. The CD4+ T cells identify this protein tag on the macrophage and produce cytokines to stimulate the production of virus-specific cytotoxic T lymphocytes (CTLs), also known as killer T cell, that denature the foreign body and exterminate the infection [7,8,19].

The interaction between HIV and immune system is different from any other virus and it is also more complex. While immune responses have the potential to fight the virus, HIV disrupts the natural immune process by directly infecting the helper T-cells. Initially immune response does get rid of a great deal of HIV, but some of it manages to survive and infect these important cells. Once the infected helper T-cells are activated, they work to create new virious instead of doing the job they are supposed to do in natural immune system. The virus uses CD4+ T cells to replicate. The CTLs then destroy the viral load and infected cells in the process. Because of this the CD4 count fall rapidly and eventually the level of virus in human body back down. Normal CD4 counts are between 800 and 1,200 cells/mm3. As the number of CD4 cells begins to fall below 200 cells per cubic millimeter of blood, the infected person will be diagnosed as having AIDS or Acquired Immune Deficiency Syndrome. As the immune system is badly damaged, the victim become vulnerable to opportunistic infections. Without treatment, people who are diagnosed with AIDS typically survive about 3 years. Once someone has a dangerous opportunistic infection, life-expectancy falls to about 1 year.

According to UNAIDS [1] and amfAR, an estimated 35.3 (32.2 - 38.8) million people were living with HIV in 2012 and 3.3 million of them are under the age of 15, globally. There were 2.3 (1.9 - 2.7) million new HIV infections globally including 260,000 under the age of 15. Every day nearly 6,300 people contract HIV i.e. nearly 262 every hour. In 2012 the number of AIDS deaths was 1.6 (1.4 - 1.9) million. 210,000 of them were under the age of 15. Since the beginning of the epidemic, more than 75 million people have contracted HIV and nearly 36 million have died of HIV-related causes. About 25 million (i.e. 70 percent of all people living with HIV worldwide) live in sub-Saharan Africa-including 88 percent of the worlds HIV-positive children. In 2012 it was estimated that 1.6 million people in this region became newly infected, and about 1.2 million adults and children died of AIDS, accounting for 75 percent of the worlds AIDS deaths. In India, the number of people living with HIV was 2,100,000 in 2012; and during the year, deaths due to AIDS was 140,000. Till now there is no complete cure from AIDS. Though HAART (i.e. highly active retroviral therapy) is very effective in blocking HIV spread within the body, it is not a cure as the viral loads readily rebound when treatment is interrupted [37]. Furthermore, ultrasensitive detection assays have shown that a low-level viraemia persists even after years of therapy [12]. This therapy can involve a complicated medication regime including many drugs and has unpleasant side effects [5,7,11,16]. The only way to prevent the disease is to avoid the contact with the virus.

Over the years, mathematical models using differential equations have been used to gain an understanding of HIV dynamics and the models have evolved over time including more parameters [27,29,34,35]. These models help us identifying important parameters and factors which have dominant effect on the development, transmission and spread of the disease. These give a better understanding for development of treatment strategies and we can draw biologically relevant interpretations. The initial models were introduced by May and Anderson [2,24]. Bonhoeffer et al. [7,8] started with a basic model and went on to include the development of drug resistant virus. Bachar and Dorfmayr [4] shows that treatment without reduction of risky behaviours may even increase the proportion of infected population. Several refinements have been added into modelling frameworks of HIV/AIDS and specific issues have been addressed by researchers [13,38]. Work has also been done on characterizing HIV dynamics incorporating virus resistance [26]. Gumel et al. [18] proposed an HIV vaccine model that considers possible vaccine induced bypass of primary infection and reversal from AIDS to chronic stage of infection together with staged progression and transmission by AIDS patients. Cai et al. [10] have developed a stage-structured HIV/AIDS epidemic model with treatment. Vergu et al.[36] showed the impact of viral diversity on the immune response and disease dynamics. Nelson et al. [25] included less than perfect drug effects and a delay in the initiation of virus production. Stafford et al.[33] suggested that cytotoxic T lymphocytes (CTLs) productively destruct infected cells and is responsible for lowering the viral load.

However, none of these papers have included immune response as a specific component. Wodarz and Nowak [39] assumed that treatment negatively affects the immune response cell in their four dimensional epidemiological model with the viral load and immune response as two different population. They conclude that, if therapy be interrupted to rebuild the immune response, the immune system would be beneficial to the long term clinical outcome of the patient. Culshaw et al. [15] consider a optimal treatment model in which the interaction between HIV and the specific immune response has been measured by the levels of natural killer cells.

Time-delays play important roles in epidemiological models, and time-delays can arise in almost every situation in epidemiology (Hethcote et al. [21]. The first model that include ‘intracellular’ delay was developed by Herz et al.[20]. They reported that including a delay, estimated value of viral clearance load was changed, but the productively infected T cell loss was unchanged. Culshaw and Ruan [14] considered the time-delay between infection of a CD4+ T-cell and the emission of viral particles. Perelson et al.[28] have considered two types of time delays: (i) pharmacological delay that occurs between the ingestion of drug and its appearance within cells and (ii) intracellular delay between initial infection of a cell by HIV and the release of new virion. The latter type of time-delay is also considered by Samanta [30,31,32]. Cai et al. [10] have considered HIV/AIDS model with two infective stages before full blown AIDS by introducing a time-delay of treatment effect, i.e. the delay between the time from the start of treatment in the symptomatic stage until the treatment effect becomes visible. Time-delay is also used to model the gestation lag, the incubation time for a infectious vector, etc. Delay differential equations exhibit much more complicated dynamics than ordinary differential equations since a time-delay could cause a stable equilibrium to become unstable and cause the population to fluctuate.

The main aim of this paper is to describe internal HIV dynamics between CD4+ T cells and the immune response cells (i.e. CTLs) of an untreated individual. We have taken three populations in our model–the healthy CD4+ T cells, the infected CD4+ T cells and CTLs. A modified functional form of infection rate have been used to get a better insight. We perform a stability analysis and supplement our theoretical results by numerical simulations. The time-lag between the infection of a CD4+ T cell and the attack of such an infected cell by CTLs has been considered. A brief discussion concludes the paper.

 

2. The basic model

Here we present an ODE model of the dynamics of an HIV-infected immune system. In this model, we consider x(t) and y(t) as the populations of uninfected and infected CD4+ T cells at time t, respectively. Usually, a mass action interaction is considered as the infection mechanism in epidemic models. But according to Hwang and Kuang [22] when the rate of infection is slow, the mass action term bxy should be replaced by . As the HIV infection progression remains very slow for a long period, we consider this modification of infection mechanism. As free virus is thought to be short lived relative to infected cells [3], we consider the viral load proportional to the infected cells. We take z(t) as the population of immune response cell (i.e. of CTLs). A single pool of immune response, proportional to the infection level, has been considered and it is dependent upon both the CD4+ T cells as well as levels of CTLs themselves. It is also assumed that there is a time-lag τ between production of the CTLs to destroy the infected CD4+ T cell and receiving the signal that a cell has been infected. They remain in this stage of development for τ units of time, decaying exponentially at the rate −μτ.

The system is defined as follows:

Here r is the source term and a is the natural death rate for healthy CD4+T cells. b is the rate at which they are being infected by HIV. c is the rate at which an infected CD4+T cell becomes infectious. The death rate of an infected CD4+T cell other than by CTLs is represented by d. p is the rate at which the infected cells are killed by CTLs. q is a generation constant for CTL pool and m represents the natural death rate for CTLs.

The rest of the paper is organized as follows. In next section positivity and boundedness of solutions of the model (2.1) is discussed. Dynamical behaviours i.e. existence of equilibrium points and stability there are presented in section 4. The stability analysis of the model in presence of time-delay investigated in section 5. In section 6, the analytical findings are verified through computer simulations. Section 7 contains the general discussions of the paper and biological implications of our mathematical findings.

 

3. Non-negativity and boundedness

Model (2.1) is a system of delay differential equations. Hence initial functions need to be specified. Let X = C([−τ, 0];ℝ3) be the Banach space of continuous mapping from [−τ, 0] → ℝ equipped with the sup-norm. By the fundamental theory of functional differential equation (FDE) we know that, there is a unique solution (x(t), y(t), z(t)) of system (2.1) with initial conditions

For biological reasons, the initial conditions are assumed to be non-negative, i.e.

The following theorem gives the criterion for boundedness and positivity of the system (2.1).

Theorem 3.1. All solutions of the system (2.1) satisfying conditions (3.1) and (3.2) are non-negative and bounded for all t ≥ 0 at which the solution exists.

Proof. Since the right hand side of system (2.1) is completely continuous and locally Lipschitzian on C, the solution (x(t), y(t), z(t)) of (2.1) with initial conditions (3.1) and (3.2) exists and is unique on [0, ζ), where 0 < ζ ≤ +∞. Now, we show that x(t) > 0 for all t ∈ [0, ζ), where 0 < ζ ≤ +∞. Otherwise, there exists a t1 ∈ [0, ζ) such that x(t1) = 0, ẋ(t1) < 0 and x(t) > 0 for all t ∈ [0, t1). From the first equation of (2.1), we have:

which is a contradiction with ẋ(t1) < 0. So x(t) > 0 for all t ≥ 0.

Next, we show that y(t) ≥ 0 for all t ∈ [0, ζ), where 0 < ζ ≤ +∞. Otherwise, there exists a t2 ∈ [0, ζ) such that y(t2) = 0, ẏ(t2) < 0 and y(t) ≥ 0 for all t ∈ [0, t2]. From the second equation of (2.1), we have:

which is a contradiction with ẏ(t2) < 0. So y(t) ≥ 0 for all t ≥ 0.

Finally, we must have z(t) ≥ 0 for all t ∈ [0, ζ), where 0 < ζ ≤ +∞. Otherwise, there exists a t3 ∈ [0, ζ) such that z(t3) = 0, ż(t3) < 0 and z(t) ≥ 0 for all t ∈ [0, t3]. From the third equation of (2.1), we have:

which is a contradiction with ż(t3) < 0. So z(t) ≥ 0 for all t ≥ 0. Therefore,

Let (x(t), y(t), z(t)) ∈ be any solution of the system (2.1). Since

we obtain

Then for any given ϵ > 0, ∃ tϵ > 0 such that x(t) < + ϵ, ∀t > tϵ.

Let us define

By non-negativity of the solutions we have

Therefore

Applying a theorem on differential inequality by Birkhoff and Rota [9] we obtain

and for t → ∞

This implies that W(t) is bounded and hence x(t), y(t), z(t) is also bounded. This completes the proof. □

 

4. Dynamics in absence of time delay (τ = 0)

4.1. Boundary equilibria and stability. Before studying the model (2.1), we consider the situation with no time-delay (i.e.τ = 0). Then the system (2.1) reduced to:

The following lemma gives all the boundary equilibrium points of the system (4.1).

Lemma 4.1. The system (4.1) always has the axial equilibrium point E1(, 0, 0). The boundary equilibrium point E2(, ŷ, 0) exists if and only if c > d. When this condition is satisfied, , ŷ are given by

The variational matrix V (E1) at E1 is given by

The eigenvalues of V (E1) are −a,−(d − c),−m. If c < d (i.e. if the rate at which an infected CD4+ T cell becomes infectious is less than the natural death rate of CTLs), then an asymptotically stable steady state exists with no infected CD4+ T cell.

The eigenvalues of V (E1) are −a,−(d − c),−m. If c < d (i.e. if the rate at which an infected CD4+ T cell becomes infectious is less than the natural death rate of CTLs), then an asymptotically stable steady state exists with no infected CD4+ T cell.

The variational matrix V (E2) at the equilibrium point E2(, ŷ, 0) is given by

The characteristic equation of V (E2) is

and the corresponding eigen values are

where

Then we have the following theorem.

Theorem 4.2. If E2(, ŷ, 0) exists and md(c(a + b) − bd)2 − c2qr2(c − d) > 0 then E2 is locally asymptotically stable in xy plane.

Proof. From lemma 4.1, if E2 exists then c > d. Then, all parameters being positive, P1 > 0 and P2 > 0. Therefore λ1, λ2 are either negative or with negative real part. Also, for the given condition λ3 < 0. Hence the theorem. □

4.2. The interior equilibrium point: its existence and stability.

Lemma 4.3. The interior equilibrium point E∗(x∗, y∗, z∗) of the system (4.1) exists if c > d and the cubic

has a positive real root greater than . Then x∗ is a positive real root of F1(x) = 0, and y∗, z∗ are given by

The variational matrix of the system (4.1) at E∗ is given by

where

The characteristic equation of V (E∗) is

where

Now we have the following theorem guaranteeing the local stability at E∗.

Theorem 4.4. If E∗(x∗, y∗, z∗) exists and A∗ = aq2x∗4+2aqmx∗2+(a+b)m2−bmqx∗2 > 0, then E∗(x∗, y∗, z∗) is locally asymptotically stable.

Proof. Since all parameters are non-negative and x∗, y∗, z∗ > 0 we have,

Hence

Then the theorem follows from the Routh-Herwitz criterion. □

 

5. Dynamics in presence of time delay (τ ≠ 0)

In presence of time-delay, the system (2.1) has the same set of boundary equilibrium points as the system (4.1). The interior equilibrium point of system (2.1) is different from E∗(x∗, y∗, z∗) of system (4.1). Here is a positive real root of the cubic

and

Obviously, c > d and are the conditions for existence of .

In the following, we study the stability behavior of .

We use the transformations : .

Then the linear system is given by

where

and all other mij, nij = 0.

The characteristic equation is given by

where

and

A necessary condition for a stability change of is that the characteristic equation (5.2) has purely imaginary solutions. Here we notice that the coefficients in P(λ, τ) and Q(λ, τ) are delay-dependent as the equilibrium components are also delay-dependent. Characteristic equations with delay-independent coefficients are comparatively simpler to deal with. The theory in such cases is well developed [17,23]. In our case, the main complication arises when we proceed to investigate the existence of a purely imaginary root λ = iω of (5.2). Here we follow the approach developed by Beretta and Kuang [6]. Let τmax be the maximum value of τ for which exists. For τ ∈ [0, τmax], we assume the following:

(i) P(0, τ) + Q(0, τ) = a3(τ) + b3(τ) ≠ 0,

(ii) P(iω, τ) + Q(iω, τ) = a3(τ) + b3(τ) − ω2(a1(τ) + b1(τ)) − i{ω3 − ωb2(τ) − ωa2(τ)} ≠ 0.

Further, we notice that

Now, it is easy to see that F(ω, τ) = |P(iω, τ)|2− |Q(iω, τ)|2 is a polynomial of degree six. Therefore,

(iv) F(ω, τ) has a finite number of zeros.

Finally, by the implicit function theorem, we have

(v) each positive simple root ω(τ) of F(ω, τ) = 0 is continuous and differentiable in τ, whenever it exists.

To obtain the stability criterion of , we set λ = iω. Substituting it in (5.2) we have the real and imaginary parts as,

where we notice that + (b3−b1ω2)2 = |Q(iω, τ)|2 ≠ 0. (Because Q(iω, τ) = 0 would imply P(iω, τ) = 0, a contradiction to (ii) above.)

Further, (5.3) gives

where .

This F(ω, τ) determines ω in terms of τ. For each τ, (5.4) has at most a finite number of real roots, which ensures that there are only a finite number of ‘gates’ for the roots to cross the imaginary axis.

Let I = {τ : τ > 0 and ω(τ) is a positive root of (5.4)}. Then, if τ ∉ I, there is no positive solution of (5.4), and consequently we have the following theorem:

Theorem 5.1. If τ ∉ I, no stability switches occur.

Now, for any τ ∈ I, (where ω(τ) is a positive simple root of (5.4)), we can define θ(τ) ∈ (0, 2π) as the solution of (5.2):

where ϕ, ψ are continuous and differentiable functions of τ such that ϕ2 + ψ2 = |P(iω, τ)|4 and |Q(iω, τ)|2 = |P(iω, τ)|2for τ ∈ I. Substituting ω = ω(τ) in (5.5), θ(τ) ∈ (0, 2π) can be determined as follows:

Here we notice that for τ ∈ I, θ(τ) defined above is continuous at τ. Furthermore if θ(τ) ∈ (0, 2π), τ ∈ I, then θ(τ) is also differentiable at τ (Beretta and Kuang 2002). Now, the relation between the arguments “θ(τ)” in (5.5) and “ω(τ)τ” in (5.3) for τ ∈ I must be

Hence, we can define the maps τn : I → ℝ+0 given by

where ω(τ) is a positive simple root of F(ω, τ) = 0. Let us introduce the functions I → ℝ

that are continuous and differentiable at τ. We notice that the values of τ(∈ I) at which stability switches may occur, are the solutions of Sn(τ) = 0 for some n ∈ ℕ0 provided the corresponding transversality condition is satisfied. To find out the transversality condition, we differentiate the characteristic equation (5.2). Then, after some algebraic manipulations, we obtain

where

Therefore,

Let us define

Then we have the following theorem:

Theorem 5.2. Let ω(τ) be a positive root of (5.4) defined for τ ∈ I, and at some τ∗ ∈ I, Sn(τ∗) = 0 for some n ∈ ℕ0. Then a pair of simple conjugate pure imaginary roots λ+(τ∗) = iω(τ∗), λ−(τ∗) = −iω(τ∗) of (5.2) exists at τ = τ∗ which crosses the imaginary axis from left to right if δ(τ∗) > 0 and crosses the imaginary axis from right to left if δ(τ∗) < 0.

Now, owing to the stability criterion of the interior equilibrium in absence of delay (given in Theorem 4.2) and the Hopf bifurcation theorem, we have the following theorem for the existence of Hopf bifurcation near .

Theorem 5.3. Let A∗ > 0. Further, let ω(τ) be a positive root of (5.4) defined for τ ∈ I, and at some τ∗ ∈ I, Sn(τ∗) = 0 for some n ∈ ℕ0. Then the system (2.1) exhibits a Hopf bifurcation near , provided that δ(τ∗) ≠ 0. As τ increases from zero upwards the first such Hopf bifurcation will be from stable to unstable, the second from unstable back to stable, the third from stable to unstable and so on.

 

6. Numerical simulations

In this section, we discuss the numerical simulations of some solutions of the systems (2.1) and (4.1) using MATLAB. If we take the parameters of the system (4.1) as r = 1.1, a = 0.4, b = 0.9, c = 0.8, d = 0.9, p = 0.4, q = 0.8, m = 0.2 and μ = 0.1 then d − c(= 0.1) being positive, all eigenvalues corresponding to the equilibrium point E1 are negative and hence E1(2.75, 0, 0) is asymptotically stable (see Fig.1). On the other hand, if we choose the parameters of system (4.1) as r = 1.1, a = 0.4, b = 0.9, c = 1.1, d = 0.8, p = 0.4, q = 0.9 and m = 1.5, we get the computed values of and ŷ as 1.7042 and 0.6391, respectively. We also observe here that md(c(a + b) − bd)2 − c2qr2(c − d) = 0.2096 > 0. So E2 is locally asymptotically stable. Fig.2 agrees with this result.

Figure 1.Here r = 1.1, a = 0.4, b = 0.9, c = 0.8, d = 0.9, p = 0.4, q = 0.8, m = 0.2 and μ = 0.1, τ = 0. It shows that, the infected CD4+ T cells (y) as well as the CTLs (z) go to extinction.

Figure 2.Here r = 1.1, a = 0.4, b = 0.9, c = 1.1, d = 0.8, p = 0.4, q = 0.9, m = 0.2 and μ = 0.1, τ = 0.. It shows that the immune response(CTLs) vanishes even the infection persists with time.

To verify the conditions of lemma 4.2, we take the parameters as r = 4, a = 0.6, b = 5, c = 4.5, d = 2.0, p = 3.1, q = 1.5 and m = 0.9. Hence we notice that c > d and x∗ = 5.8229 is a positive root of F1(x) = 0, which is greater than xm = 0.69. All the conditions of lemma 4.2 being satisfied, the interior equilibrium point of system (4.1) exists and is given by E∗(5.8229, 0.1030, 0.7812). Also for these choices of parameters, A∗ = 1382 > 0, and consequently by theorem 4.2 E∗ is locally asymptotically stable. The corresponding phase portrait is shown in Fig.3(a). Fig.3(b) shows the stable behaviour of x, y, z with time.

Figure 3.Here r = 4, a = 0.6, b = 5, c = 4.5, d = 2.0, p = 3.1, q = 1.5, m = 0.9. (a) Phase portrait of system (4.1) showing that E∗ is locally asymptotically stable. (b) For x(0) = 6, y(0) = 0.5, z(0) = 0.5, it is seen that x, y, z approach to their equilibrium values in finite time.

In presence of delay, we have considered the same values of parameters (as in Fig.3) except τ (≠ 0). We observe that for τ = 0.06 < τ∗ = 0.075, = (5.8171, 0.1038, 0.7812) is stable (see Fig.4(a) and Fig.4(b)). For τ = 0.09 > τ∗, is unstable and there is a periodic orbit near = (5.8142, 0.1041, 0.7811) (see Fig.5(a)). The oscillations of x, y, z with time is shown in Fig.5(b).

Figure 4.Here r = 4, a = 0.6, b = 5, c = 4.5, d = 2.0, p = 3.1, q = 1.5, m = 0.9, μ = 0.1 and τ = 0.06. (a) Phase portrait of system (2.1) showing that is locally asymptotically stable. (b) For x(0) = 6, y(0) = 0.5, z(0) = 0.5, it is seen that x, y, z approach to their equilibrium values in finite time.

Figure 5.Here r = 4, a = 0.6, b = 5, c = 4.5, d = 2.0, p = 3.1, q = 1.5, m = 0.9, μ = 0.1 but τ = 0.09. (a) Phase portrait of system (2.1) showing that is unstable, and there is a limit cycle which grows out of . (b) Oscillations of x, y, z with time when x(0) = 6, y(0) = 0.5, z(0) = 0.5.

 

7. Concluding remarks

In this paper we have studied a HIV infection model with immune response cells as a specific compartment. The time-lag between a cell being infected and corresponding immune response cell being activated, has been considered. The assumptions made in formulation of the model has been described in details in section 2. Positivity and boundedness, established in section 3, ensures that the model is biologically well behaved. Throughout the discussion we have observed that the model possesses the axial equilibrium E1 and this equilibrium is locally asymptotically stable if and only if the death rate (killed by CTLs) of infected helper T cells is grater than the infection conversion rate (i.e. d > c). On the other hand, if c > d and md(c(a + b) − bd)2 − c2qr2(c − d) > 0 then the system stabilizes to the immune response free equilibrium. It is undesirable that the infected CD4+T cells are in positive level but CTLs goes to extinction. Hence the parameters must be controlled in such a manner that the CTL-free equilibrium becomes unstable.

The dynamics in absence of time-delay has been discussed in section 4. Lemma 4.2 and theorem 4.2 describe the existence and stability criterion of the system with τ = 0. Using Routh-Herwitz criterion, it has been shown that for A∗ > 0 the system with no delay approaches the co-existence equilibrium point trough stable path in finite time. To analyze the system with delay (τ ≠ 0), we have followed the approach developed by Beretta and Kuang [6]. It is shown that the delay play a key role in stability of the model. There is a critical value τ∗ of τ, so that a stable equilibrium becomes unstable as τ crosses τ∗. The analytical findings have been verified by computer simulations.