Abstract and Applied Analysis

Volume 2014 (2014), Article ID 838396, 18 pages

http://dx.doi.org/10.1155/2014/838396

Research Article

Hopf Bifurcation and Stability of Periodic Solutions for Delay Differential Model of HIV Infection of CD4⁺ T-cells

P. Balasubramaniam,¹ M. Prakash,¹ Fathalla A. Rihan,2,3 and S. Lakshmanan²

¹Department of Mathematics, Gandhigram Rural Institute-Deemed University, Gandhigram, Tamil Nadu 624 302, India

²Department of Mathematical Sciences, College of Science, UAE University, P.O. Box 15551, Al-Ain, UAE

³Department of Mathematics, Faculty of Science, Helwan University, Cairo 11795, Egypt

Received 17 February 2014; Revised 13 June 2014; Accepted 19 June 2014; Published 31 August 2014

Academic Editor: Cemil Tunç

Copyright © 2014 P. Balasubramaniam et al. 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.

Abstract.

This paper deals with stability and Hopf bifurcation analyses of a mathematical model of HIV infection of T-cells. The model is based on a system of delay differential equations with logistic growth term and antiretroviral treatment with a discrete time delay, which plays a main role in changing the stability of each steady state. By fixing the time delay as a bifurcation parameter, we get a limit cycle bifurcation about the infected steady state. We study the effect of the time delay on the stability of the endemically infected equilibrium. We derive explicit formulae to determine the stability and direction of the limit cycles by using center manifold theory and normal form method. Numerical simulations are presented to illustrate the results.

1. Introduction

Since 1980, the human immunodeficiency virus (HIV) or the associated syndrome of opportunistic infections that causes acquired immunodeficiency syndrome (AIDS) has been considered as one of the most serious global public health menaces. When HIV enters the body, its main target is the CD4 lymphocytes, also called CD4 T-cells (including T-cells). When a CD4 cell is infected with HIV, the virus goes through multiple steps to reproduce itself and create many more virus particles. The AIDS term, which is known as the late stage of HIV, covers the range of infections and illnesses which can result from a weakened immune system caused by HIV. Based on the clinical studies, it is known that, for a normal person, the T-cells count is around and for HIV infected patient it gradually decreases to or below, which leads to AIDS. However, this may take several years for the number of CD4 T-cells to reduce to a level where the immune system is weakened [1–6].

Mathematical models, usingdelay differential equations (DDEs), have provided insights in understanding the dynamics of HIV infection. Discrete or continuous time delays have been introduced to the models to describe the time between infection of a T-cell and the emission of viral particles on a cellular level [7–13]. In general, DDEs exhibit much more complicated dynamics than ODEs since the time delay could cause a stable equilibrium to become unstable and cause the populations to fluctuate [14–16]. In studying the viral clearance rates, Perelson et al. [17] 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 virion. In this paper, we incorporate an intracellular delay to the model to describe the time between infection of a T-cell and the emission of viral particles on a cellular level [18]. We study the impact of the presence of such time delay on the dynamics of the model.

The outline of the present paper is as follows. In Section 2, we describe the model. In Section 3, we study the qualitative behavior of the model via stability of the steady states and Hopf bifurcation when time delay is considered as a bifurcation parameter. In Section 4, we provide an explicit formula to determine the direction of bifurcating periodic solution by applying center manifold theory and normal form method. We provide some numerical simulations to demonstrate the effectiveness of the analysis in Section 5 and we conclude in Section 6.

2. Description of the Model

Let us start the analysis with some basic models of the dynamics of target (uninfected) cells and infected T-cells by HIV. As a first approximation, the dynamics between HIV and the macrophage population was described by the simplest model of infection dynamics presented in [19–21]. Denoting uninfected cells by and infected cells by and assuming that viruses are transmitted mainly by cell to cell contact, the model is given by

The target (uninfected) T-cells are produced at a rate , die at a rate , and become infected by virus at a rate . The infected host cells die at a rate . The basic reproductive ratio of the virus is then given by . If there is no infection or if , there is only trivial equilibrium () with no virus-producing cells. Whereas if , the virus can establish an infection and the system converges to the equilibrium with both uninfected cells and infected cells, .

However, in most viral infections, the CTL response plays a crucial part in antiviral defence by attacking viral infected cells [22, 23]. As the the cytotoxic T-lymphocyte (CTL) immune response is necessary to eliminate or control the viral infection, we incorporated the antiviral CTL immune response into the basic model (1). Therefore, if we add CTL response, which is denoted by , into model (1) (see [19]), then the extended model is

Thus, CTLs proliferate in response to antigen at a rate , die at a rate , and lyse infected cells at a rate . We assume that the CTL pool consists of two populations: the precursors and the effectors . In other words, we assume that there are primary and secondary responses to viral infections. Then, the model (2) becomes

The infected cells are killed by CTL effector cells at a rate . Upon contact with antigen, CTLp proliferate at a rate and differentiate into effector cells CTLe at a rate . CTL precursors die at a rate , and effectors die at a rate ; see Figure 1.

Figure 1: A simplified model of virus-CTL interaction. The virus dynamics is described by the basic model of Nowak and Bangham [19]. The uninfected target cells are produced at a rate and die at a rate . They become infected by the virus at a rate . The infected cells produce new virus particle and die at a rate . When recognize antigen on the surface of infected cells, they become activated and expand at a rate , decay at a rate , and differentaite into efector cells at a rate . The effector cells lyse the infected cells at a rate .

Since the proliferation of T-cells is density dependent, that is, the rate of proliferation decreases as T-cells increase and reach the carrying capacity, we then extend the above basic viral infection model to include the density dependent growth of the T-cell population (see [24–26]). It is also known that HIV infection leads to low levels of T-cells via three main mechanisms: direct viral killing of infected cells, increased rates of apoptosis in infected cells, and killing of infected T-cells by cytotoxic T-lymphocytes [26]. Hence, it is reasonable to include apoptosis of infected cells. An average of viral particles is produced by infected cells per day. The treatment with single antiviral drug is considered to be failed, so that the combination of antiviral drugs is needed for the better treatment [25]. Therefore, in the below revised model, we combine the antiretroviral drugs, namely, reverse transcriptase inhibitor (RTI) and protease inhibitor (PI) to make the model realistic (see [27–29]). RTIs can block the infection of target T-cells by infectious virus, and PIs cause infected cells to produce noninfectious virus particles. The modified model takes the form

The first equation of model (4) represents the rate of change in the count of healthy T-cells that produced at rate and become infected at rate , with the mortality . We assume that the uninfected T-cells proliferate logistically, thus the growth rate is multiplied by the term and this term approaches zero when the total number of T-cells approaches the carrying capacity . The effects of combination of RTI and PI antiviral drugs are represented by the term , where , , represents the effects of RTI and , , represents the effects of PI. The second equation of model (4) denotes the rate of change in the count of infected T-cells. The infected T-cells decay at a rate and denotes apoptosis rate of infected cell; infected cells are killed by CTL effectors at a rate . The third equation of the model denotes the rate of change in the CTLp population; proliferation rate of the CTLp is given by and is proportional to the infected cells ; CTLp die at a rate and differentiate into CTL effectors at a rate . The last equation of the model represents the concentration of CTL effectors, which die at a rate . In reality, the specific immune system is not immediately effective following invasion by a novel pathogen. There may be an explicit time delay between infection and immune initiation and there may be a gradual build-up in immune efficacy during which the immune response develops, before reaching maximal specificity to the pathogen ([8, 30, 31]). In order to make model (4) more realistic, time delay in the immune response should be included in the following model:

The range of parameter values of the model are given in Table 1.

Table 1: Parameter definitions and estimations used in the underlying model.

Parameter
Notes
Estimated Value
Range
Source

Source of uninfected CD4⁺ T-cells
10
0–10
[26]

Rate of infection
0.1
0.00001–0.5
[26]

Total carrying capacity
1500
1500
[26]

Logistic growth term
0.03
0.03–3
[26]

Mortality rate of CD4⁺ T-cells
0.06
0.007–0.1
[26]

Antiretroviral (RTI) therapy
0.9
0-1
see text

Infected cells died out naturally
0.3
0.2–1.4
[26]

Apoptosis rate of infected cells
0.2
0.2
[26]

Clearance rate of infected cells
1
0.001–1
[26]

Protease inhibitor therapy
0.9
0, 1
see text

Rate of differentiation of CTLs
0.02
Assumed
—

Death rate of CTL precursors
0.02
0.005–0.15
[26]

Proliferation of CTLs responsiveness
0.1
0.001–1
[26]

Mortality rate or CTL effectors
0.1
0.005–0.15
[26]

We start our analysis by presenting some notations that will be used in the sequel. Let be the Banach space of continuous functions mapping the interval into , where ; the initial conditions are given by

where are smooth functions, for all .   From the fundamental theory of functional differential equations (see [32, 33]), it is easy to see that the solutions of system (5) with the initial conditions as stated above exist for all and are unique. It can be shown that these solutions exist for all and stay nonnegative. In fact, if , then for all . The same argument is true for the ,  , and components. Hence, the interior is invariant for system (5).

3. Steady States

We can obtain the steady state values by setting . The steady state value of the infection-free steady sate is given by , while the infected steady state is given by

and is given by the following quadratic equation:

where , , .

3.1. Stability and Hopf Bifurcation Analysis of Infected Steady State

In order to study full dynamics of model (4) by using time delay as a bifurcation parameter, we need to linearize the model around the steady state and determine the characteristic equation of the Jacobian matrix. The roots of the characteristic equation determine the asymptotic stability and existence of Hopf bifurcation for the model. The characteristic equation of the linearized system is given by

which is equivalent to the equation

where and

Let us consider the following equation:

For the nondelayed model (say ), from (10), we have

where

Lemma 1. For , the unique nontrivial equilibrium is locally asymptotically stable if the real parts of all the roots of (13) are negative.

Proof. The proof of the above lemma is based on holding the following conditions: ,  ,  , and , as proposed by Routh-Hurwitz criterion. We conclude that equilibrium is locally asymptotically stable if and only if all the roots of the characteristic equation (13) have negative real parts which depends on the numerical values of parameters that are shown in the numerical exploration.

3.2. Existence of Hopf Bifurcation

We here study the impact of the time-delay parameter on the stability of HIV infection of T-cells. We deduce criteria that ensure the asymptotic stability of infected steady state , for all . We arrive at the following theorem.

Theorem 2. Necessary and sufficient conditions for the infected equilibrium to be asymptotically stable for all delay are as follows

(i)the real parts of all the roots of are negative;(ii)for all and , , where .

Proof. Assume that Lemma 1 is true. Now, for , we have

Substituting into (5) and separating the real and imaginary parts of the equations yields

After some mathematical manipulations, we obtain the following equations

Let

From (16), we have

where

The conditions and of Theorem 2 hold if and only if (19) has no real positive root.

Let ; then (19) takes the form

If , then (19) has at least one positive root. In the case when (19) has four positive roots, we have

From (16), we have

where and ; we choose .

To establish Hopf bifurcation at , we need to show that

By differentiating (10) with respect to , we can get

It follows that

Then, by combining (10), we get

Substituting in (27) (where and ) yields

where

Thus,

Notice that

By summarizing the above analysis, we arrive at the following theorem.

Theorem 3. The infected equilibrium of the system (5) is asymptotically stable for and it undergoes Hopf bifurcation at .

4. Direction and Stability of Bifurcating Periodic Solutions

In the previous section, we obtained conditions for Hopf bifurcation to occur when , . It is also important to derive explicit formulae from which we can determine the direction, stability, and period of periodic solutions bifurcating around the infected equilibrium at the critical value . We use the cafeteria of normal forms and center manifold proposed by Hassard [34]. We assume that the model (5) undergoes Hopf bifurcation at the infected equilibrium when , , and then are the corresponding purely imaginary roots of the characteristic equation at the infected equilibrium . Assume also that

then using Taylors expansion for system (3) at the equilibrium point yields

Here,

For convenience, let and for . Denote ; has -order continuous derivative. For initial conditions , (33) can be rewritten as

where ,  , and are given, respectively, by

is one parameter family of bounded linear operators in and

From the discussion in the above section, we know that if , then model (5) undergoes a Hopf bifurcation at the infected equilibrium , and the associated characteristic equation of model (5) has a pair of purely imaginary roots . By Reisz representation, there exists a function of bounded variation for such that

In fact, we can choose

where is Dirac delta function. Next, for , define

Since , (35) can be written as

where ,  . For , the adjoint operator of can be defined as

For and , in order to normalize the eigenvalues of operator and adjoint operator , the following bilinear form is defined by

where and is complex conjugate of . It can verify that and are adjoint operators with respect to this bilinear form.

We assume that are eigenvalues of and the other eigenvalues have strictly negative real parts. Thus, they are also eigenvalues of . Now we compute the eigenvector of corresponding to the eigenvalue and the eigenvector of corresponding to the eigenvalue . Suppose that is eigenvector of associated with ; then, . It follows from the definition of and (36), (38), and (40) that

Solving (45), we can easily obtain , where

Similarly, suppose that the eigenvector of corresponding to is , . By the definition of and (36), (38), and (40), one gets

Solving (47), we easily obtain , where

In order to assure that , we need to determine the value of . From (44), one gets

Let

On the center manifold , we have

where

and are local coordinates of the center manifold in the direction of and , respectively. Note that is real if is real. So we only consider real solutions. From (50), we obtain

For the solution of (35), from (41) and (44), since , we have

Rewrite (54) as

where

Substituting (42) and (54) into (50) yields

which can be written as

where

On the center manifold , we have

Substituting (52) and (55) into (61), one obtains

Substituting (52) and (60) into (59) yields

Comparing the coefficients of (62) and (63), one gets

Since , then we have

Thus, we obtain

It is obvious that

So

also

and hence

It follows from (54) that

where

Since , we have

Comparing the coefficients of the above equation with those in (61), we have

We need to compute and for . Equations (62) and (63) imply that

Comparing the coefficients of the above equation with (60), we have

It follows from (40) and (64) that

By solving the above equation for and for , one obtains

where and can be determined by setting in .

In fact, we have

comparing the coefficients of the above equations with those in (61), it follows that

By the definition of and (40) and (64), we get

One can notice that

Thus, we obtain

where , ; the above equation can be written as

From (78), (84), we can calculate , and we can derive the following parameters:

We arrive at the following theorem.

Theorem 4. The periodic solution is supercritical (subcritical) if   ; the bifurcating periodic solutions are orbitally asymptotically stable with asymptotical phase (unstable) if   ; the period of the bifurcating periodic solution increases (decreases) if   .

5. Numerical Simulations

In this section, we provide some simulations of model (4) to exhibit the impact of discrete time delay in the model. We consider the parameters values: , , , , , , , , , , , , , and . According to the given parameters’ values, the threshold critical value from the formula (21) exists. The steady state exists and is asymptotically stable (see Figure 1). We may notice that the solution converges to the equilibrium with damping oscillations as the value of increases. Once the delay crosses the critical value , then the model shows the existence of Hopf bifurcation which is depicted in the Figure 2. In Figure 3, we consider the efficacy of antiretroviral value is 0.9, which may be responsible for the loss of stability. The asymptotic behavior to the infection-free steady state, when we consider antiviral treatment (with ,  , and time delay ), is shown in Figure 4. According to Theorem 4, the parameters , , , , and are estimated. Based on these values one can conclude that bifurcating periodic solutions are unstable and decreases in the period of bifurcating periodic solutions. The existence of periodic solution is subcritical. For numerical treatment of DDEs and related issues; we refer the readers to [35, 36].

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 2: Each panel (from (a) to (h)) shows the time evolution and trajectory of model (4) when (critical value) and the effect of therapies is considered to be and . It shows that the endemic steady state of model is asymptotically stable.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 3: It shows the numerical simulations of model (4), when the time delay of immune activation exceeds the critical value, . The endemic steady state of the model undergoes Hopf bifurcation; stability switch and periodic solutions appear.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 4: It shows the numerical simulations of model (4), when the efficacy rate of antiretroviral treatments is considered to be low; that is, and . It shows that the equilibrium of the model undergoes Hopf bifurcation with oscillatory behavior in solutions even though the delay value is less than the critical value .

Several packages and types of software are available for the numerical integration and/or the study of bifurcations in delay differential equations (see, e.g., [37, 38]. In this paper we utilize MIDDE code [39]) which is suitable to simulate stiff and nonstiff delay differential equations and Volterra delay integrodifferential equations, using monoimplicit RK methods.

6. Concluding Remarks

In this manuscript, we provided a conceptual T-cell infection model which includes the logistic growth term along with two different types of antiretroviral drug therapies. The model includes a discrete time delay in the immune activation response, which plays an important role in the dynamics of the model. The infection-free and endemic steady states of the model are determined (Figure 5). The stability of steady states is analyzed. We deduced a formula that determines the critical value (branch value) . Necessary and sufficient conditions for the equilibrium to be asymptotically stable for all positive delay values are proved. We have seen that if the time delay exceeds the critical value , model (4) undergoes a Hopf bifurcation. The direction and stability of bifurcating periodic solutions are deduced in explicit formulae, using center manifold and normal forms. We also presented some numerical simulations to the underlying model to investigate the obtained results and theory. We have seen also that the antiretroviral treatments help to increase the level of uninfected T-cells. The theoretical results that were confirmed by the numerical simulations show that the delayed CTL response can lead to complex bifurcations, and, in particular, the coexistence of multiple stable periodic solutions. When the time delay exceeds the critical (threshold) value, we may get subcritical behaviour that leads to a loss of uninfected T-cells.

(a)

(b)

(c)

(d)

Figure 5: It shows the numerical simulations model (4) when the efficacy rate of antiretroviral treatment is at expected level, and , and the delay value exceeds the critical value . The solution always lies within the feasible region and the infection-free steady state is asymptotically stable.

Conflict of Interests

The authors declare that they have no competing interests for this paper.

Acknowledgments

This work is supported by United Arab Emirates University, UAE, (Grant no. NRF-7-20886) and the National Board for Higher Mathematics, Mumbai (Grant no. 2/48/2012/NBHM/R&D-II/11020). The authors would like to thank Professor Cemil Tunc and referees for their valuable comments.

References

D. Liu and B. Wang, “A novel time delayed HIV/AIDS model with vaccination & antiretroviral therapy and its stability analysis,” Applied Mathematical Modelling, vol. 37, no. 7, pp. 4608–4625, 2013.

I. Ncube, “Absolute stability and Hopf bifurcation in a Plasmodium falciparum malaria model incorporating discrete immune response delay,” Mathematical Biosciences, vol. 243, no. 1, pp. 131–135, 2013.

X. Song, S. Wang, and J. Dong, “Stability properties and Hopf bifurcation of a delayed viral infection model with lytic immune response,” Journal of Mathematical Analysis and Applications, vol. 373, no. 2, pp. 345–355, 2011.

X. Song, S. Wang, and X. Zhou, “Stability and Hopf bifurcation for a viral infection model with delayed non-lytic immune response,” Journal of Applied Mathematics and Computing, vol. 33, no. 1-2, pp. 251–265, 2010.

T. Wang, Z. Hu, and F. Liao, “Stability and Hopf bifurcation for a virus infection model with delayed humoral immunity response,” Journal of Mathematical Analysis and Applications, vol. 411, no. 1, pp. 63–74, 2014.

T. Wang, Z. Hu, F. Liao, and W. Ma, “Global stability analysis for delayed virus infection model with general incidence rate and humoral immunity,” Mathematics and Computers in Simulation, vol. 89, pp. 13–22, 2013.

Z. Bai and Y. Zhou, “Dynamics of a viral infection model with delayed CTL response and immune circadian rhythm,” Chaos, Solitons & Fractals, vol. 45, no. 9-10, pp. 1133–1139, 2012.

E. Beretta, M. Carletti, D. E. Kirschner, and S. Marino, “Stability analysis of a mathematical model of the immune response with delays,” in Mathematics for Life Science and Medicine, pp. 177–206, Springer, Berlin, Germany, 2007.

Z. Hu, J. Zhang, H. Wang, W. Ma, and F. Liao, “Dynamics analysis of a delayed viral infection model with logistic growth and immune impairment,” Applied Mathematical Modelling, vol. 38, no. 2, pp. 524–534, 2014.

G. Huang, H. Yokoi, Y. Takeuchi, and T. Sasaki, “Impact of intracellular delay, immune activation delay and nonlinear incidence on viral dynamics,” Japan Journal of Industrial and Applied Mathematics, vol. 28, no. 3, pp. 383–411, 2011.

F. A. Rihan and D. H. Abdel Rahman, “Delay differential model for tumor-immune dynamics with HIV infection of CD4⁺ T-cells,” International Journal of Computer Mathematics, vol. 90, no. 3, pp. 594–614, 2013.

J. Tam, “Delay effect in a model for virus replication,” Journal of Mathemathics Applied in Medicine and Biology, vol. 16, no. 1, pp. 29–37, 1999.

Z. Wang and R. Xu, “Stability and Hopf bifurcation in a viral infection model with nonlinear incidence rate and delayed immune response,” Communications in Nonlinear Science and Numerical Simulation, vol. 17, no. 2, pp. 964–978, 2012.

J. E. Mittler, B. Sulzer, A. U. Neumann, and A. S. Perelson, “Influence of delayed viral production on viral dynamics in HIV-1 infected patients,” Mathematical Biosciences, vol. 152, no. 2, pp. 143–163, 1998.

F. A. Rihan, D. H. Abdel Rahman, and S. Lakshmanan, “A time delay model of tumour-immune system interactions: global dynamics, parameter estimation, sensitivity analysis,” Applied Mathematics and Computation, vol. 232, pp. 606–623, 2014.

L. Zhang, C. Zhang, and D. Zhao, “Hopf bifurcation analysis of integro-differential equation with unbounded delay,” Applied Mathematics and Computation, vol. 217, no. 10, pp. 4972–4979, 2011.

A. S. Perelson, A. U. Neumann, M. Markowitz, J. M. Leonard, and D. D. Ho, “HIV-1 dynamics in vivo: virion clearance rate, infected cell life-span, and viral generation time,” Science, vol. 271, no. 5255, pp. 1582–1586, 1996.

A. V. M. Herz, S. Bonhoeffer, R. M. Anderson, R. M. May, and M. A. Nowak, “Viral dynamics in vivo: limitations on estimates of intracellular delay and virus decay,” Proceedings of the National Academy of Sciences of the United States of America, vol. 93, no. 14, pp. 7247–7251, 1996.

M. A. Nowak and C. R. M. Bangham, “Population dynamics of immune responses to persistent viruses,” Science, vol. 272, no. 5258, pp. 74–79, 1996.

R. M. Anderson and R. M. May, “Population biology of infectious diseases: part I,” Nature, vol. 280, no. 5721, pp. 361–367, 1979.

R. J. de Boer and A. S. Perelson, “Target cell limited and immune control models of HIV infection: a comparison,” Journal of Theoretical Biology, vol. 190, no. 3, pp. 201–214, 1998.

D. Wodarz and M. A. Nowak, “Specific therapy regimes could lead to long-term immunological control of HIV,” Proceedings of the National Academy of Sciences of the United States of America, vol. 96, no. 25, pp. 14464–14469, 1999.

D. Wodarz, K. M. Page, R. A. Arnaout, A. R. Thomsen, J. D. Lifson, and M. A. Nowak, “A new theory of cytotoxic T-lymphocyte memory: implications for HIV treatment,” Philosophical Transactions of the Royal Society B: Biological Sciences, vol. 355, no. 1395, pp. 329–343, 2000.

C. Lv, L. Huang, and Z. Yuan, “Global stability for an {HIV}-1 infection model with Beddington-DeAngelis incidence rate and {CTL} immune response,” Communications in Nonlinear Science and Numerical Simulation, vol. 19, no. 1, pp. 121–127, 2014.

A. S. Perelson and P. W. Nelson, “Mathematical analysis of HIV-1 dynamics in vivo,” SIAM Review, vol. 41, no. 1, pp. 3–44, 1999.

Y. Wang, Y. Zhou, F. Brauer, and J. M. Heffernan, “Viral dynamics model with {CTL} immune response incorporating antiretroviral therapy,” Journal of Mathematical Biology, vol. 67, no. 4, pp. 901–934, 2013.

M. Pitchaimani, C. Monica, and M. Divya, “Stability analysis for hiv infection delay model with protease inhibitor,” Biosystems, vol. 114, no. 2, pp. 118–124, 2013.

H. Shu and L. Wang, “Role of cd4⁺ t-cell proliferation in HIV infection u nder antiretroviral therapy,” Journal of Mathematical Analysis and Applications, vol. 394, no. 2, pp. 529–544, 2012.

S. Wang and Y. Zhou, “Global dynamics of an in-host {HIV} 1-infection model with the long-lived infected cells and four intracellular delays,” International Journal of Biomathematics, vol. 5, no. 6, Article ID 1250058, 2012.

A. Fenton, J. Lello, and M. B. Bonsall, “Pathogen responses to host immunity: the impact of time delays and memory on the evolution of virulence,” Proceedings of the Royal Society B: Biological Sciences, vol. 273, no. 1597, pp. 2083–2090, 2006.

M. A. Nowak, R. M. May, and K. Sigmund, “Immune responses against multiple epitopes,” Journal of Theoretical Biology, vol. 175, no. 3, pp. 325–353, 1995.

J. Hale, Theory of Functional Differential Equations, Springer, New York, NY, USA, 1997.

Y. Kuang, Delay Differential Equations with Applications in Population Dynamics, vol. 191 of Mathematics in Science and Engineering, Academic Press, San Diego, Calif, USA, 1993.

B. D. Hassard, Theory and Applications of Hopf Bifurcation, vol. 41, CUP Archive, 1981.

C. T. H. Baker, G. Bocharov, A. Filiz et al., Numerical modelling by delay and volterra functional differential equations, Topics in Computer Mathematics and its Applications, LEA Athens, Hellas, 1999

F. A. Rihan, Numerical treatment of delay differential equations in bioscience [Ph.D. thesis], The University of Manchester, Manchester, UK, 2000.

L. F. Shampine and S. Thompson, “Solving DDEs in Matlab,” Applied Numerical Mathematics, vol. 37, no. 4, pp. 441–458, 2001.

C. Paul, “A user-guide to archi: an explicit runge-kutta code for so lving delay and neutral differential equations and parameter estimation problems,” MCCM Techical Report 283, University of Manchester, 1997.

F. A. Rihan, E. H. Doha, M. I. Hassan, and N. M. Kamel, “Numerical treatments for Volterra delay integro-differential equations,” Computational Methods in Applied Mathematics, vol. 9, no. 3, pp. 292–308, 2009.