Global Stability and Hopf Bifurcation for a Virus Dynamics Model with General Incidence Rate and Delayed CTL Immune Response

Abdoul Samba Ndongo

Department of Mathematics, Faculty of Sciences, University of Nouakchott Al Aasriya, Nouakchott, Mauritania.

DOI: 10.4236/am.2021.1211068                

Abstract

In this work, we investigate an HIV-1 infection model with a general incidence rate and delayed CTL immune response. The model admits three possible equilibria, an infection-free equilibrium E^*₀, CTL-inactivated infection equilibrium E^*₁ and CTL-activated infection equilibrium E^*₂. We prove that in the absence of CTL immune delay, the model has exactly the basic behaviour model, for all positive intracellular delays, the global dynamics are determined by two threshold parameters R₀ and R₁, if R₀ ≤ 1, E^*₀ is globally asymptotically stable, if R₁ ≤ 1 < R₀, E^*₁ is globally asymptotically stable and if R₁ >1, E^*₂ is globally asymptotically stable. But if the CTL immune response delay is different from zero, then the behaviour of the model at E^*₂ changes completely, although R₁ > 1, a Hopf bifurcation at E^*₂ is established. In the end, we present some numerical simulations.

Keywords

Virus Dynamics, Intracellular and CTL Immune Response Delays, Lyapunov Function, Global Asymptotic Stability, Hopf Bifurcation

1. Introduction

Immunity can be broadly categorized into adaptive immunity and innate immunity. Adaptive immunity is mediated by clonally B-cells and T-cells are also called lymphocytes, humoral and cellular immunity is characterized by specificity and memory. Cytotoxic T Lymphocytes (CTLs) play an important role in antiviral defense by attacking infected cells. Many mathematical models have been developed to describe HIV-1 (human immunodeficiency virus type 1) (see for example [1] - [15] ). First of all, we introduced the standard viral infection model with CTL immune response considered by Nowak and Bangham [16] as follows:

$\begin{array}{l}\frac{\text{d}x\left(t\right)}{\text{d}t}=s-dx\left(t\right)-\beta x\left(t\right)v\left(t\right),\\ \frac{\text{d}y\left(t\right)}{\text{d}t}=\beta x\left(t\right)v\left(t\right)-\delta y\left(t\right)-py\left(t\right)z\left(t\right),\\ \frac{\text{d}v\left(t\right)}{\text{d}t}=ky\left(t\right)-uv\left(t\right),\\ \frac{\text{d}z\left(t\right)}{\text{d}t}=cy\left(t\right)z\left(t\right)-bz\left(t\right),\end{array}$ (1.1)

where $x\left(t\right)$, $y\left(t\right)$, $v\left(t\right)$ and $z\left(t\right)$ represent the densities of uninfected target cells, infected cells, virus and CTL cells at time t, respectively. Uninfected cells are produced at rate s, die at rate d, and become infected cells at rate $\beta xv$. Infected cells are produced from uninfected cells at rate $\beta x\left(t\right)v\left(t\right)$ and die at rate $\delta$. The parameter p accounts for the strength of the lytic component. Free virions are produced from infected cells at rate $ky$ and are removed at rate u. The parameter b is the death rate for CTLs, and $cyz$ describes the rate of CTL immune response activated by the infected cells. Recently many studies have been done to improve the model (1.1) by introducing delays and changing the incidence rate according to different practical backgrounds. These studies used different delayed models with different forms of incidence rates. In [17], Wang

etal. used the model (1.1) and replaced the incidence function rate by a Beddington-DeAngelis functional response of the form $\frac{\beta xv}{1+{\alpha }_{1}x+{\alpha }_{2}v}$, which was introduced by Beddington [18] and DeAngelis et al. [19]. In [10], Yuan et al. have

presented a model HIV-1 with an incidence rate of the form $\beta x\left(t\right)f\left(v\left(t\right)\right)$. This incidence rate considered in this paper generalized many forms of commonly used incidence rate, including simple mass action, saturation incidence rate, Beddington-DeAngelis functional response form and the Crowly-Martin functional response form introduced by Crowly-Martin (see [20] ). The global dynamic of a virus dynamics model with Crowly-Martin functional response was discussed in [21]. In this paper, we propose a virus dynamics model with three delays and general incidence rate as follows:

$$\begin{array}{l}\frac{\text{d}x\left(t\right)}{\text{d}t}=s-dx\left(t\right)-f\left(x\left(t\right),v\left(t\right)\right)v\left(t\right),\\ \frac{\text{d}y\left(t\right)}{\text{d}t}={\text{e}}^{-m_1{\tau }_1}f\left(x\left(t-{\tau }_1\right),v\left(t-{\tau }_1\right)\right)v\left(t-{\tau }_1\right)-\delta y\left(t\right)-py\left(t\right)z\left(t\right),\\ \frac{\text{d}v\left(t\right)}{\text{d}t}={\text{e}}^{-m_2{\tau }_2}ky\left(t-{\tau }_2\right)-uv\left(t\right),\\ \frac{\text{d}z\left(t\right)}{\text{d}t}=cy\left(t-{\tau }_3\right)z\left(t-{\tau }_3\right)-bz\left(t\right),\end{array}$$ (1.2)

where the parameters have the same meanings as in system (1.1). ${\tau }_1$ represents the time between viral entry into a target cell and the production of new virus particles, ${\tau }_2$ represent the time necessary for the newly produced viruses to become mature and then infection, that is, the maturation time of newly produce viruses and ${\tau }_3$ represents the time that antigenic stimulation needs for generating immunity response [22]. The terms ${\text{e}}^{-m_1{\tau }_1}$, ${\text{e}}^{-m_2{\tau }_2}$ denote the surviving rate of infected cells and viruses during the delay period, respectively. The function f is assumed to be continuously differentiable in the interior of ${\Re }_{+}^2$ and satisfies the following hypotheses:

Ÿ $f\left(0,v\right)=0,\text{for}\text{all}v\ge 0,$ (H1)

Ÿ $\frac{\partial f\left(x,v\right)}{\partial x}>0,\text{for}\text{all}x>0\text{and}v\ge 0,$ (H2)

Ÿ $\frac{\partial f\left(x,v\right)}{\partial v}\le 0,\text{for}\text{all}x>0\text{and}v\ge 0.$ (H3)

Ÿ $\left(1-\frac{f\left(x\mathrm{,}v\right)}{f\left(x\mathrm{,}{v}^{\prime }\right)}\right)\left(\frac{f\left(x\mathrm{,}{v}^{\prime }\right)}{f\left(x\mathrm{,}v\right)}-\frac{v}{v^{\prime }}\right)\le \mathrm{0,}\text{for}\text{all}x\mathrm{,}v\mathrm{,}{v}^{\prime }>0.$ (H4)

Under the hypotheses (H₁) - (H₄) we show that, for $E_0^{\mathrm{<mo>\; *\; </mo>}}$ and $E_1^{\mathrm{<mo>\; *\; </mo>}}$, the dynamics of the model (1.2) is exactly like as for the standard basic model, if ${\tau }_1\mathrm{,}{\tau }_2\mathrm{,}{\tau }_3\ge 0$, as well as for $E_2^{\mathrm{<mo>\; *\; </mo>}}$ if ${\tau }_3=0$ and ${\tau }_1,{\tau }_2>0$. But for $E_2^{\mathrm{<mo>\; *\; </mo>}}$, if in particular, ${\tau }_1={\tau }_2=0$ and ${\tau }_3\ne 0$, we show that the dynamics of the model changes completely. $E_2^{\mathrm{<mo>\; *\; </mo>}}$ remains locally stable if ${\tau }_3$ remains below a certain critical value ${\tau }_3^{*}>0$. Moreover, in ${\tau }_3^{\mathrm{<mo>\; *\; </mo>}}$, there is appearance of a periodic phenomenon, obtained by Hopf bifurcation from $E_2^{\mathrm{<mo>\; *\; </mo>}}$. When ${\tau }_3$ passes ${\tau }_3^{\mathrm{<mo>\; *\; </mo>}}$, $E_2^{\mathrm{<mo>\; *\; </mo>}}$ becomes unstable. Note that the incidence rate considered in this paper generalizes many forms of commonly used incidence rate, including simple mass action and saturation incident rate and Beddington-DeAngelis functional response form and Crowly-Martin functional response form.

The present paper is organized as follows. In Section 2, we derived the basic reproduction ratios for viral infection and CTL immune response $R_0$ and $R_1$, respectively, and the existence of three equilibria is established. By means of suitable Lyapunov functionals and LaSalle’s invariant principle, we proved that, if $R_0\le 1$, the infection-free equilibrium $E_0^{\mathrm{<mo>\; *\; </mo>}}$ is globally asymptotically stable for all ${\tau }_1\mathrm{,}{\tau }_2\mathrm{,}{\tau }_3\ge 0$ (Section 3). In Section 4 we study the global stability of $E_1^{\mathrm{<mo>\; *\; </mo>}}$ and $E_2^{\mathrm{<mo>\; *\; </mo>}}$ and we prove that if $R_1\le 1<R_0$, the CTL-inactivated infection equilibrium $E_1^{\mathrm{<mo>\; *\; </mo>}}$ is globally asymptotically stable for all $${\tau }_1\mathrm{,}{\tau }_2\mathrm{,}{\tau }_3\ge 0$$, and if $R_1>1$ and ${\tau }_3=0$, the CTL-activated infection equilibrium $E_2^{\mathrm{<mo>\; *\; </mo>}}$ is globally asymptotically stable for all ${\tau }_1\mathrm{,}{\tau }_2\ge 0$. In Section 5, we analyze the local stability of $E_2^{\mathrm{<mo>\; *\; </mo>}}$ in the case ${\tau }_3\ge 0$, ${\tau }_1={\tau }_2=0$ and the existence of Hopf bifurcation from $E_2^{\mathrm{<mo>\; *\; </mo>}}$. In Section 6, we present some numerical simulations to illustrate our results. In Section 7, some examples are given.

The initial conditions of (1.2) are given as:

$\begin{array}{l}x\left(\theta \right)={\varphi }_1\left(\theta \right),y\left(\theta \right)={\varphi }_2\left(\theta \right),v\left(\theta \right)={\varphi }_3\left(\theta \right),z\left(\theta \right)={\varphi }_4\left(\theta \right)\\ {\varphi }_i\left(\theta \right)\ge 0,i=1,2,3,4,\theta \in \left(-\tau ,0\right],\tau =\max \left\{{\tau }_1,{\tau }_2,{\tau }_3\right\},\end{array}$ (1.3)

where $\varphi =\left({\varphi }_1\mathrm{,}{\varphi }_2\mathrm{,}{\varphi }_3\mathrm{,}{\varphi }_4\right)\in C^{+}$, here $C^{+}=C\left(\left(-\tau ,0\right],{\Re }_{+}^4\right)$, with $C\left(\left(-\tau \mathrm{,0}\right]\mathrm{,}{\Re }^4\right)$ denotes the Banach space of continuous functions mapping the interval $\left(-\tau \mathrm{,0}\right]$ into ${\Re }^4$.

2. Preliminary Results

In this section we established the positivity and the boundedness of solutions of (1.2) and we define the basic reproduction numbers $R_0$ and $R_1$ and the existence of three possible equilibrium points is studied. The following theorem establishes the non-negativity and boundedness of solutions of (1.2).

Theorem 2.1. Under hypothesis (H₁) and with the initial conditions (1.3), all solutions $\left(x\left(t\right)\mathrm{,}y\left(t\right)\mathrm{,}v\left(t\right)\mathrm{,}z\left(t\right)\right)$ of system (1.2) are non-negative and bounded on $\left[\mathrm{0,}+\infty \right)$.

Proof. Let put system (1.2) in a vector form $\stackrel{\dot{}}{Z}\left(t\right)=G\left(Z_t\right)$ by setting $Z={\left(x\mathrm{,}y\mathrm{,}v\mathrm{,}z\right)}^{\text{T}}$, and for all $t\ge 0$.

$G\left(Z_t\right)=\left(\begin{array}{c}s-f\left(x\left(t\right),v\left(t\right)\right)v\left(t\right)-dx\left(t\right)\\ {\text{e}}^{-m_1{\tau }_1}f\left(x\left(t-{\tau }_1\right),v\left(t-{\tau }_1\right)\right)v\left(t-{\tau }_1\right)-\delta y\left(t\right)-py\left(t\right)z\left(t\right)\\ {\text{e}}^{-m_2{\tau }_2}ky\left(t-{\tau }_2\right)-uv\left(t\right)\\ cy\left(t-{\tau }_3\right)z\left(t-{\tau }_3\right)-bz\left(t\right)\end{array}\right),$

where $Z_t\left(\theta \right)=Z\left(t+\theta \right)$ and $G\mathrm{<mo>\; :\; </mo>}C\left(\left(-\tau \mathrm{,0}\right]\mathrm{,}{\Re }_{+}^4\right)\to {\Re }^4$. It is easy to check that ${G_i\left(Z\right)|}_{Z_i=0}\ge 0$, $i=1,2,3,4$. Due to lemma 2 in [23], any solution of (1.2) with $Z\left(\theta \right)\in C_{+}$, say $Z\left(t\right)=Z\left(t\mathrm{,}Z\left(\theta \right)\right)$, is such that $Z\left(t\right)\in {\Re }_{+}^4$ for all $t\ge 0$.

Next we show that the solutions are also bounded.

From (1.2)₁ we have $\frac{\text{d}x}{\text{d}t}\le s-dx\left(t\right)$. This implies $\lim {\sup }_{t\to \infty }x\left(t\right)\le \frac{s}{d}$, so $x\left(t\right)$ is bounded.

Let:

$F\left(t\right)={\text{e}}^{-m_1{\tau }_1}x\left(t-{\tau }_1\right)+y\left(t\right)+\frac{p}{c}z\left(t+{\tau }_3\right).$

$\frac{\text{d}F\left(t\right)}{\text{d}t}={\text{e}}^{-m_1{\tau }_1}\left(s-dx\left(t-{\tau }_1\right)\right)-\delta y\left(t\right)-\frac{pb}{c}z\left(t+{\tau }_3\right)\le {\text{e}}^{-m_1{\tau }_1}s-\gamma F\left(t\right),$

where $\gamma =\min \left\{d,\delta ,b\right\}$. This implies that $F\left(t\right)$ is bounded and so is $y\left(t\right)$ and $z\left(t\right)$. Thus, there exists a $\chi >0$ such that $y\left(t\right)\le \chi$. It follows from (1.2)₃ that:

$\frac{\text{d}v\left(t\right)}{\text{d}t}\le {\text{e}}^{-m_2{\tau }_2}k\chi -uv\left(t\right),$

and consequently $v\left(t\right)$ is bounded. Finally, all the solutions of system (1.2) are bounded. This completes the proof.

Global behaviour of system (1.2) may depends on the basic reproduction numbers $R_0$ and $R_1$ given by:

$R_0=\frac{k{\text{e}}^{-m_1{\tau }_1}{\text{e}}^{-m_2{\tau }_2}}{\delta u}f\left(x_0^{*},0\right),$ (2.1)

where, $x_0^{*}=\frac{s}{d}$, and:

$R_1=\frac{k{\text{e}}^{-m_1{\tau }_1}{\text{e}}^{-m_2{\tau }_2}}{\delta u}f\left(M,\frac{kb{\text{e}}^{-m_2{\tau }_2}}{uc}\right),$ (2.2)

with $M=x_0^{*}-\frac{\delta b}{dc{\text{e}}^{-m_1{\tau }_1}}$. Here, $R_0$ and $R_1$ are the basic reproduction ratios

for viral infection and CTL immune response of system (1.2), respectively. Based on the hypotheses (H₂) and (H₃) it is clear that $R_1<R_0$.

System (1.2) always has the infection-free equilibrium $E_0^{\mathrm{<mo>\; *\; </mo>}}=\left(x_0^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,0,0,0}\right)$. For other possible equilibrium , we have the following theorem:

Theorem 2.2. Suppose that the conditions (H₁) - (H₃) are satisfied.

1) If $R_0>1$, then system (1.2) has a CTL-inactivated infection equilibrium of the form $E_1^{*}=\left(x_1^{*},y_1^{*},v_1^{*},0\right)$ with $x_1^{\mathrm{<mo>\; *\; </mo>}}\in \left(\mathrm{0,}{x}_0^{\mathrm{<mo>\; *\; </mo>}}\right)$.

2) If $R_1>1$, then system (1.2) has a CTL-activated infection equilibrium of the form $E_2^{\mathrm{<mo>\; *\; </mo>}}=\left(x_2^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{y}_2^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{v}_2^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{z}_2^{\mathrm{<mo>\; *\; </mo>}}\right)$ with $x_2^{\mathrm{<mo>\; *\; </mo>}}\in \left(\mathrm{0,}M\right)$.

Proof. The steady states of system (1.2) satisfy the following equations:

$\begin{array}{l}s-f\left(x,v\right)v-dx=0,\\ {\text{e}}^{-m_1{\tau }_1}f\left(x,v\right)v-\delta y-pyz=0,\\ {\text{e}}^{-m_2{\tau }_2}ky-uv=0,\\ cyz-bz=0.\end{array}$ (2.3)

From the last equation of (2.3), we have:

$\left(cy-b\right)z=0.$ (2.4)

Equations (2.4) has two possible solutions, $z=0$ or $cy-b=0$.

If $z=0$, (2.3)₃ yields $y=\frac{u}{k{\text{e}}^{-m_2{\tau }_2}}v$.

By substituting this into (2.3)₂, we obtain:

$\left({\text{e}}^{-m_1{\tau }_1}f\left(x,v\right)-\frac{\delta u}{k{\text{e}}^{-m_2{\tau }_2}}\right)v=0,$ (2.5)

which gives $v=0$ or $f\left(x,v\right)=\frac{\delta u}{k{\text{e}}^{-m_1{\tau }_1}{\text{e}}^{-m_2{\tau }_2}}$.

If $v=0$, we obtain the infection-free equilibrium $E_0^{\mathrm{<mo>\; *\; </mo>}}=\left(\frac{s}{d}\mathrm{,0,0,0}\right)$.

If $v\ne 0$, (2.3)₁ and (2.3)₂ yields:

$y=\frac{\left(s-dx\right){\text{e}}^{-m_1{\tau }_1}}{\delta }:=y\left(x\right),$ (2.6)

with Equation (2.3)₃, yields:

$v=\frac{k{\text{e}}^{-m_1{\tau }_1}{\text{e}}^{-m_2{\tau }_2}\left(s-dx\right)}{\delta u}:=v\left(x\right).$ (2.7)

Since, $y\ge 0$ and $v\ge 0$ this implies that $x\le \frac{s}{d}$.

Now, from (H₁), (H₂) and (H₃), the following functional:

$K\left(x\right)=f\left(x\mathrm{,}v\left(x\right)\right)-\frac{\delta u}{k{\text{e}}^{-m_1{\tau }_1}{\text{e}}^{-m_2{\tau }_2}}\mathrm{,}$ (2.8)

satisfy:

$K\left(0\right)K\left(\frac{s}{d}\right)=-{\left(\frac{\delta u}{k{\text{e}}^{-m_1{\tau }_1}{\text{e}}^{-m_2{\tau }_2}}\right)}^2\left(R_0-1\right)<0$ for $R_0>1$

and:

$\stackrel{\dot{}}{K}\left(x\right)=\frac{\partial f}{\partial x}-\frac{kd{\text{e}}^{-m_1{\tau }_1}{\text{e}}^{-m_2{\tau }_2}}{\delta u}\frac{\partial f}{\partial v}>0$.

Hence, we obtain the CTL-inactivated infection equilibrium:

$E_1^{*}=\left(x_1^{*},y_1^{*},v_1^{*},0\right)=\left(x_1^{*},y\left(x_1^{*}\right),v\left(x_1^{*}\right),0\right),$

with $x_1^{\mathrm{<mo>\; *\; </mo>}}$ is the unique zero in $\left(\mathrm{0,}{x}_0^{\mathrm{<mo>\; *\; </mo>}}\right)$ of K and y and v are given by (2.6), (2.7) respectively.

If $z\ne 0$, from (2.4) and (2.3)₃ we obtain:

$y=\frac{b}{c}:=y_2^{*}$,

and:

$v=\frac{kb{\text{e}}^{-m_2{\tau }_2}}{uc}:=v_2^{*}$,

and from the first and second equation of (2.3) we obtain:

$z=z\left(x\right):=\frac{{\text{e}}^{-m_1{\tau }_1}dc\left(M-x\right)}{pb}\ge 0,$ (2.9)

which implies that $x\le M$.

Now, from (2.3)₁ the functional:

$L\left(x\right):=s-dx-f\left(x,v_2^{*}\right)v_2^{*}=0,$

satisfy:

$L\left(0\right)L\left(M\right)=s\frac{\delta b}{c{\text{e}}^{-m_1{\tau }_1}}\left(1-R_1\right)<0$ for $R_1>1$,

$\stackrel{\dot{}}{L}\left(x\right)=-\frac{kb{\text{e}}^{-m_2{\tau }_2}}{uc}\frac{\partial f}{\partial x}-d<0.$

Hence, we obtain the CTL-activated infection equilibrium of the form $E_2^{\mathrm{<mo>\; *\; </mo>}}=\left(x_2^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{y}_2^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{v}_2^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}z\left(x_2^{\mathrm{<mo>\; *\; </mo>}}\right)\right)$ with $x_2^{\mathrm{<mo>\; *\; </mo>}}$ is the unique zero of L in $\left(\mathrm{0,}M\right)$ and z is given by (2.9). This completes the proof.

Remark 1. From (2.8) we have $K\left(M\right)=\frac{\delta u}{k{\text{e}}^{-m_1{\tau }_1}{\text{e}}^{-m_2{\tau }_2}}\left(R_1-1\right)\le 0$ if $R_1\le 1$. So, as K is increasing in the interval $\left[\mathrm{0,}\frac{s}{d}\right]$ ; we deduces that $M\le x_1^{\mathrm{<mo>\; *\; </mo>}}$ and consequently $y_1^{\mathrm{<mo>\; *\; </mo>}}-\frac{b}{c}\le 0$.

3. Global Stability of the Infection-Free Equilibrium

In this section, we study the global stability of the infection-free equilibrium $E_0^{\mathrm{<mo>\; *\; </mo>}}$ of System (1.2).

Theorem 3.1. If $R_0\le 1$, the infection-free equilibrium $E_0^{\mathrm{<mo>\; *\; </mo>}}$ of System (1.2) is globally asymptotically stable for all ${\tau }_1,{\tau }_2,{\tau }_3\ge 0$.

Proof. Define a Lyapunov function $U_0$ as follows:

$\begin{array}{c}{U}_0\left(t\right)=x\left(t\right)-x_0^{*}+{\displaystyle {\int }_{x_0^{*}}^x}\frac{f\left(x_0^{*},0\right)}{f\left(\sigma ,0\right)}\text{d}\sigma +{\text{e}}^{m_1{\tau }_1}y\left(t\right)+\frac{\delta {\text{e}}^{m_1{\tau }_1}{\text{e}}^{m_2{\tau }_2}}{k}v\left(t\right)\\ +\frac{p{\text{e}}^{m_1{\tau }_1}}{c}z\left(t\right)+{\displaystyle {\int }_{t-{\tau }_1}^t}f\left(x\left(\sigma \right),v\left(\sigma \right)\right)v\left(\sigma \right)\text{d}\sigma \\ +{\text{e}}^{m_1{\tau }_1}\delta {\displaystyle {\int }_{t-{\tau }_2}^t}y\left(\sigma \right)\text{d}\sigma +{\text{e}}^{m_1{\tau }_1}p{\displaystyle {\int }_{t-{\tau }_3}^t}y\left(\sigma \right)z\left(\sigma \right)\text{d}\sigma \end{array}$

Calculating the time derivative of $U_0\left(t\right)$ along the positive solutions of System (1.2), we obtain:

$\begin{array}{c}\frac{\text{d}{U}_0}{\text{d}t}=\left(1-\frac{f\left(x_0^{*},0\right)}{f\left(x\left(t\right),0\right)}\right)\stackrel{\dot{}}{x}\left(t\right)+{\text{e}}^{m_1{\tau }_1}\stackrel{\dot{}}{y}\left(t\right)+\frac{{\text{e}}^{m_1{\tau }_1}{\text{e}}^{m_2{\tau }_2}\delta }{k}\stackrel{\dot{}}{v}\left(t\right)+\frac{p{\text{e}}^{m_1{\tau }_1}}{c}\stackrel{\dot{}}{z}\left(t\right)\\ +f\left(x\left(t\right),v\left(t\right)\right)v\left(t\right)-f\left(x\left(t-{\tau }_1\right),v\left(t-{\tau }_1\right)\right)v\left(t-{\tau }_1\right)+{\text{e}}^{m_1{\tau }_1}\delta y\left(t\right)\\ -{\text{e}}^{m_1{\tau }_1}\delta y\left(t-{\tau }_2\right)+{\text{e}}^{m_1{\tau }_1}py\left(t\right)z\left(t\right)-{\text{e}}^{m_1{\tau }_1}py\left(t-{\tau }_3\right)z\left(t-{\tau }_3\right)\end{array}$

$\begin{array}{l}=\left(1-\frac{f\left(x_0^{*},0\right)}{f\left(x\left(t\right),0\right)}\right)\left(s-dx\right)+\frac{f\left(x_0^{*},0\right)}{f\left(x\left(t\right),0\right)}f\left(x\left(t\right),v\left(t\right)\right)v\left(t\right)\\ -\frac{u\delta {\text{e}}^{m_1{\tau }_1}{\text{e}}^{m_2{\tau }_2}}{k}v\left(t\right)-\frac{pb{\text{e}}^{m_1{\tau }_1}}{c}z\left(t\right)\\ =d{x}_0^{*}\left(1-\frac{f\left(x_0^{*},0\right)}{f\left(x\left(t\right),0\right)}\right)\left(1-\frac{x}{x_0^{*}}\right)\\ +\frac{u\delta {\text{e}}^{m_1{\tau }_1}{\text{e}}^{m_2{\tau }_2}}{k}v\left(t\right)\left(\frac{f\left(x\left(t\right),v\left(t\right)\right)}{f\left(x\left(t\right),0\right)}{R}_0-1\right)-\frac{pb{\text{e}}^{m_1{\tau }_1}}{c}z\left(t\right)\\ \le d{x}_0^{*}\left(1-\frac{f\left(x_0,0\right)}{f\left(x\left(t\right),0\right)}\right)\left(1-\frac{x}{x_0}\right)+\frac{u\delta }{k}v\left(t\right)\left(R_0-1\right)-\frac{p}{c}z\left(t\right).\end{array}$

From (H₁), we have $\left(1-\frac{f\left(x_0^{*},0\right)}{f\left(x\left(t\right),0\right)}\right)\left(1-\frac{x}{x_0^{\mathrm{<mo>\; *\; </mo>}}}\right)\le 0$. Since $R_0\le 1$, we have $\frac{\text{d}{U}_0}{\text{d}t}\le 0$. Hence, the infection-free equilibrium $E_0^{\mathrm{<mo>\; *\; </mo>}}$ is stable under the condition $R_0\le 1$. Let $\Gamma$ be the largest invariant set in the set:

$E=\left\{\left(x\mathrm{,}y\mathrm{,}v\mathrm{,}z\right)\mathrm{<mo>\; |\; </mo>}\frac{\text{d}{U}_0}{\text{d}t}=0\right\}$

Note that for each $t\ge 0$, $\frac{\text{d}{U}_0}{\text{d}t}=0$ if and only if $x=x_0^{*}$, $z=0$ and $v\left(t\right)\left(\frac{f\left(x\left(t\right),v\left(t\right)\right)}{f\left(x\left(t\right),0\right)}{R}_0-1\right)=0$. Therefore we have two cases:

• $R_0<1$, then $v=0$ ;

• $R_0=1$. From $x=x_0^{*}$ and the first equation of (1.2), we have $f\left(x_0^{*},v\right)v=0$. Since $x_0^{*}>0$ from (H₁) and (H₂) we get $f\left(x_0^{*},v\right)>0$. Hence $v=0$.

By the above discussion, we deduce that $\Gamma =\left\{E_0^{*}\right\}$. It follow from LaSalle invariance principle [24] that the infection-free equilibrium $E_0^{\mathrm{<mo>\; *\; </mo>}}$ is globally asymptotically stable.

4. Global Stability of the Infected Equilibria

In this section, we study the global stability of the CTL-inactivated infection equilibrium $E_1^{\mathrm{<mo>\; *\; </mo>}}$ and the CTL-activated infection equilibrium $E_2^{\mathrm{<mo>\; *\; </mo>}}$ of system (1.2) by the Lyapunov direct method.

We set:

$g\left(x\right)=x-1-\ln x,\text{for}x\in \left(0,\infty \right).$

It is clear that for any $x>0$, $g\left(x\right)\ge 0$ and $g\left(x\right)$ has the global minimum $x=1$, with $g\left(1\right)=0$.

Theorem 4.1. Suppose that the conditions (H₁) - (H₄) are satisfied. Then the equilibrium $E_1^{\mathrm{<mo>\; *\; </mo>}}$ is globally asymptotically stable if $R_1\le 1<R_0$ for all ${\tau }_1,{\tau }_2,{\tau }_3\ge 0$, and the equilibrium $E_2^{\mathrm{<mo>\; *\; </mo>}}$ is globally asymptotically stable if $R_1>1$ for all ${\tau }_1,{\tau }_2\ge 0$ and ${\tau }_3=0$.

Proof. Let $E^{\mathrm{<mo>\; *\; </mo>}}=\left(x^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{y}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{v}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{z}^{\mathrm{<mo>\; *\; </mo>}}\right)$ an infected equilibrium point of (1.2), i.e., $E^{*}=E_1^{*}$ or $E_2^{\mathrm{<mo>\; *\; </mo>}}$. Define a Lyapunov functional:

$U_1=V_1+V_2,$

where:

$\begin{array}{c}{V}_1={\displaystyle {\int }_{x^{*}}^x}\left(1-\frac{f\left(x^{*},v^{*}\right)}{f\left(\sigma ,v^{*}\right)}\right)\text{d}\sigma +{\text{e}}^{m_1{\tau }_1}{y}^{*}g\left(\frac{y}{y^{*}}\right)\\ +\frac{f\left(x^{*},v^{*}\right)v^{*}{}^2}{{\text{e}}^{-m_2{\tau }_2}k{y}^{*}}g\left(\frac{v}{v^{*}}\right)+\frac{p}{{\text{e}}^{-m_1{\tau }_1}c}{z}^{*}g\left(\frac{z}{z^{*}}\right),\end{array}$

$\begin{array}{c}{V}_2=f\left(x^{*},v^{*}\right)v^{*}{\displaystyle {\int }_{t-{\tau }_1}^t}g\left(\frac{f\left(x\left(\sigma \right),v\left(\sigma \right)\right)v\left(\sigma \right)}{f\left(x^{*},v^{*}\right)v^{*}}\right)\text{d}\sigma \\ +f\left(x^{*},v^{*}\right)v^{*}{\displaystyle {\int }_{t-{\tau }_2}^t}g\left(\frac{y\left(\sigma \right)}{y^{*}}\right)\text{d}\sigma +{\text{e}}^{m_1{\tau }_1}p{\displaystyle {\int }_{t-{\tau }_3}^t}y\left(\sigma \right)z\left(\sigma \right)\text{d}\sigma ,\end{array}$

The function, $U\mathrm{<mo>\; :\; </mo>}x\to {\displaystyle {\int }_{x^{\mathrm{<mo>\; *\; </mo>}}}^x}\left(1-\frac{f\left(x^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{v}^{\mathrm{<mo>\; *\; </mo>}}\right)}{f\left(\sigma \mathrm{,}{v}^{\mathrm{<mo>\; *\; </mo>}}\right)}\right)\text{d}\sigma$, verifies:

$\stackrel{\dot{}}{U}\left(x\right)=1-\frac{f\left(x^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{v}^{\mathrm{<mo>\; *\; </mo>}}\right)}{f\left(x\mathrm{,}{v}^{\mathrm{<mo>\; *\; </mo>}}\right)}\mathrm{.}$

From (H₂), we have $\stackrel{\dot{}}{U}\left(x\right)<0$ for $x\in \left(\mathrm{0,}{x}^{\mathrm{<mo>\; *\; </mo>}}\right)$, $\stackrel{\dot{}}{U}\left(x\right)>0$ for $x\in \left(x^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}\infty \right)$ and $\stackrel{\dot{}}{U}\left(x^{\mathrm{<mo>\; *\; </mo>}}\right)=0$. So $U\left(x\right)\ge 0$. Consequently $U_2$ is nonnegative defined with respect to the endemic equilibrium $E^{\mathrm{<mo>\; *\; </mo>}}$, which is a global minimum.

We now prove that the time derivative of $U_2$ is non-positive. Calculating the time derivative of $V_1$ along the positive solutions of (1.2), we obtain:

$\begin{array}{c}{\stackrel{\dot{}}{V}}_1=\left(1-\frac{f\left(x^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{v}^{\mathrm{<mo>\; *\; </mo>}}\right)}{f\left(x\mathrm{,}{v}^{\mathrm{<mo>\; *\; </mo>}}\right)}\right)\stackrel{\dot{}}{x}+{\text{e}}^{m_1{\tau }_1}\left(1-\frac{y^{*}}{y}\right)\stackrel{\dot{}}{y}\\ +\frac{f\left(x^{*},v^{*}\right)v^{*}}{{\text{e}}^{-m_2{\tau }_2}k{y}^{*}}\left(1-\frac{v^{*}}{v}\right)\stackrel{\dot{}}{v}+\frac{p}{{\text{e}}^{-m_1{\tau }_1}c}\left(1-\frac{z^{*}}{z}\right)\stackrel{\dot{}}{z}\\ =\left(1-\frac{f\left(x^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{v}^{\mathrm{<mo>\; *\; </mo>}}\right)}{f\left(x\mathrm{,}{v}^{\mathrm{<mo>\; *\; </mo>}}\right)}\right)\left(s-f\left(x,v\right)v-dx\right)\\ +{\text{e}}^{m_1{\tau }_1}\left(1-\frac{y^{*}}{y}\right)\left({\text{e}}^{-m_1{\tau }_1}f\left(x_{{\tau }_1},v_{{\tau }_1}\right)v_{{\tau }_1}-\delta y-pyz\right)\\ +\frac{f\left(x^{*},v^{*}\right)v^{*}}{{\text{e}}^{-m_2{\tau }_2}k{y}^{*}}\left(1-\frac{v^{*}}{v}\right)\left({\text{e}}^{-m_2{\tau }_2}k{y}_{{\tau }_2}-uv\right)+\frac{p}{{\text{e}}^{-m_1{\tau }_1}c}\left(1-\frac{z^{*}}{z}\right)\left(c{y}_{{\tau }_3}{z}_{{\tau }_3}-bz\right).\end{array}$

At $E^{\mathrm{<mo>\; *\; </mo>}}$, by using $s=d{x}^{\mathrm{<mo>\; *\; </mo>}}+f\left(x^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{v}^{\mathrm{<mo>\; *\; </mo>}}\right)v^{\mathrm{<mo>\; *\; </mo>}}$, $\delta y^{\mathrm{<mo>\; *\; </mo>}}={\text{e}}^{-m_1{\tau }_1}f\left(x^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{v}^{\mathrm{<mo>\; *\; </mo>}}\right)v^{\mathrm{<mo>\; *\; </mo>}}-p{y}^{\mathrm{<mo>\; *\; </mo>}}{z}^{\mathrm{<mo>\; *\; </mo>}}$ and $\frac{u}{k}=\frac{{\text{e}}^{-m_2{\tau }_2}{y}^{*}}{v^{*}}$, we have:

$\begin{array}{c}{\stackrel{\dot{}}{V}}_1=\left(1-\frac{f\left(x^{*},v^{*}\right)}{f\left(x,v^{*}\right)}\right)\left(d{x}^{*}-dx\right)+f\left(x^{*},v^{*}\right)v^{*}\left(1-\frac{f\left(x^{*},v^{*}\right)}{f\left(x,v^{*}\right)}\right)-f\left(x,v\right)v\\ +\frac{f\left(x^{*},v^{*}\right)}{f\left(x,v^{*}\right)}f\left(x,v\right)v+f\left(x_{{\tau }_1},v_{{\tau }_1}\right)v_{{\tau }_1}-\frac{f\left(x^{*},v^{*}\right)v^{*}}{y^{*}}y-\frac{y^{*}}{y}f\left(x_{{\tau }_1},v_{{\tau }_1}\right)v_{{\tau }_1}\\ +f\left(x^{*},v^{*}\right)v^{*}+\frac{f\left(x^{*},v^{*}\right)v^{*}}{y^{*}}{y}_{{\tau }_2}-f\left(x^{*},v^{*}\right)v-\frac{f\left(x^{*},v^{*}\right)v^{*}}{y^{*}}\frac{v^{*}}{v}{y}_{{\tau }_2}\\ +f\left(x^{*},v^{*}\right)v^{*}+{\text{e}}^{m_1{\tau }_1}p\left(y^{*}-\frac{b}{c}\right)\left(z-z^{*}\right)+{\text{e}}^{m_1{\tau }_1}p\left(\frac{y_{{\tau }_3}{z}_{{\tau }_3}}{z}-y\right)\left(z-z^{*}\right).\end{array}$

Calculating the time derivative of $V_2$, we obtain:

$\begin{array}{c}{\stackrel{\dot{}}{V}}_2=f\left(x^{*},v^{*}\right)v^{*}\left(\frac{f\left(x,v\right)v}{f\left(x^{*},v^{*}\right)v^{*}}-\frac{f\left(x_{{\tau }_1},v_{{\tau }_1}\right)v_{{\tau }_1}}{f\left(x^{*},v^{*}\right)v^{*}}+\ln \frac{f\left(x_{{\tau }_1},v_{{\tau }_1}\right)v_{{\tau }_1}}{f\left(x,v\right)v}\right)\\ +f\left(x^{*},v^{*}\right)v^{*}\left(\frac{y}{y^{*}}-\frac{y_{{\tau }_2}}{y^{*}}+\ln \frac{y_{{\tau }_2}}{y}\right)+{\text{e}}^{m_1{\tau }_1}pyz-{\text{e}}^{m_1{\tau }_1}p{y}_{{\tau }_3}{z}_{{\tau }_3}\\ =f\left(x,v\right)v-f\left(x_{{\tau }_1},v_{{\tau }_1}\right)v_{{\tau }_1}+f\left(x^{*},v^{*}\right)v^{*}\ln \frac{f\left(x_{{\tau }_1},v_{{\tau }_1}\right)v_{{\tau }_1}}{f\left(x,v\right)v}\\ +\frac{f\left(x^{*},v^{*}\right)v^{*}}{y^{*}}y-\frac{f\left(x^{*},v^{*}\right)v^{*}}{y^{*}}{y}_{{\tau }_2}\\ +f\left(x^{*},v^{*}\right)v^{*}\ln \frac{y_{{\tau }_2}}{y}+{\text{e}}^{m_1{\tau }_1}pyz-{\text{e}}^{m_1{\tau }_1}p{y}_{{\tau }_3}{z}_{{\tau }_3}.\end{array}$

Combining (4.1) and (4.2), we obtain:

$\begin{array}{c}{\stackrel{\dot{}}{U}}_1=\left(1-\frac{f\left(x^{*},v^{*}\right)}{f\left(x,v^{*}\right)}\right)\left(d{x}^{*}-dx\right)+f\left(x^{*},v^{*}\right)v^{*}\left(1-\frac{f\left(x^{*},v^{*}\right)}{f\left(x,v^{*}\right)}\right)\\ +\frac{f\left(x^{*},v^{*}\right)}{f\left(x,v^{*}\right)}f\left(x,v\right)v-\frac{y^{*}}{y}f\left(x_{{\tau }_1},v_{{\tau }_1}\right)v_{{\tau }_1}+f\left(x^{*},v^{*}\right)v^{*}-f\left(x^{*},v^{*}\right)v\\ -\frac{f\left(x^{*},v^{*}\right)v^{*}}{y^{*}}\frac{v^{*}}{v}{y}_{{\tau }_2}+f\left(x^{*},v^{*}\right)v^{*}+f\left(x^{*},v^{*}\right)v^{*}\ln \frac{f\left(x_{{\tau }_1},v_{{\tau }_1}\right)v_{{\tau }_1}}{f\left(x,v\right)v}\\ +f\left(x^{*},v^{*}\right)v^{*}\ln \frac{y_{{\tau }_2}}{y}+{\text{e}}^{m_1{\tau }_1}p\left(y^{*}-\frac{b}{c}\right)\left(z-z^{*}\right)+{\text{e}}^{m_1{\tau }_1}p\left(y-\frac{y_{{\tau }_3}{z}_{{\tau }_3}}{z}\right)z^{*}\end{array}$

$\begin{array}{l}=d\left(1-\frac{f\left(x^{*},v^{*}\right)}{f\left(x,v^{*}\right)}\right)\left(x^{\mathrm{<mo>\; *\; </mo>}}-x\right)+f\left(x^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{v}^{\mathrm{<mo>\; *\; </mo>}}\right)v^{\mathrm{<mo>\; *\; </mo>}}\left(\frac{f\left(x,v^{*}\right)}{f\left(x,v\right)}-\frac{v}{v^{\mathrm{<mo>\; *\; </mo>}}}\right)\left(1-\frac{f\left(x,v\right)}{f\left(x,v^{*}\right)}\right)\\ +f\left(x^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{v}^{\mathrm{<mo>\; *\; </mo>}}\right)v^{\mathrm{<mo>\; *\; </mo>}}\left(1-\frac{f\left(x^{*},v^{*}\right)}{f\left(x,v^{*}\right)}+\ln \frac{f\left(x^{*},v^{*}\right)}{f\left(x,v^{*}\right)}\right)\\ +f\left(x^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{v}^{\mathrm{<mo>\; *\; </mo>}}\right)v^{\mathrm{<mo>\; *\; </mo>}}\left(1-\frac{f\left(x,v^{*}\right)}{f\left(x,v\right)}+\ln \frac{f\left(x,v^{*}\right)}{f\left(x,v\right)}\right)\\ +f\left(x^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{v}^{\mathrm{<mo>\; *\; </mo>}}\right)v^{*}\left(1-\frac{f\left(x_{{\tau }_1},v_{{\tau }_1}\right)v_{{\tau }_1}{y}^{*}}{f\left(x^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{v}^{\mathrm{<mo>\; *\; </mo>}}\right)v^{*}y}+\ln \frac{f\left(x_{{\tau }_1},v_{{\tau }_1}\right)v_{{\tau }_1}{y}^{*}}{f\left(x^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{v}^{\mathrm{<mo>\; *\; </mo>}}\right)v^{*}y}\right)\\ +f\left(x^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{v}^{\mathrm{<mo>\; *\; </mo>}}\right)v^{*}\left(1-\frac{y_{{\tau }_2}{v}^{*}}{y^{*}v}+\ln \frac{y_{{\tau }_2}{v}^{*}}{y^{*}v}\right)\end{array}$

$\begin{array}{l}+{\text{e}}^{m_1{\tau }_1}p\left(y^{*}-\frac{b}{c}\right)\left(z-z^{*}\right)+{\text{e}}^{m_1{\tau }_1}p\left(y-\frac{y_{{\tau }_3}{z}_{{\tau }_3}}{z}\right)z^{*}\\ =d\left(1-\frac{f\left(x^{*},v^{*}\right)}{f\left(x,v^{*}\right)}\right)\left(x^{*}-x\right)+f\left(x^{*},v^{*}\right)v^{*}\left(\frac{f\left(x,v^{*}\right)}{f\left(x,v\right)}-\frac{v}{v^{*}}\right)\left(1-\frac{f\left(x,v\right)}{f\left(x,v^{*}\right)}\right)\\ -f\left(x^{*},v^{*}\right)v^{*}g\left(\frac{f\left(x^{*},v^{*}\right)}{f\left(x,v^{*}\right)}\right)-f\left(x^{*},v^{*}\right)v^{*}g\left(\frac{f\left(x,v^{*}\right)}{f\left(x,v\right)}\right)\\ -f\left(x^{*},v^{*}\right)v^{*}g\left(\frac{f\left(x_{{\tau }_1},v_{{\tau }_1}\right)v_{{\tau }_1}{y}^{*}}{f\left(x^{*},v^{*}\right)v^{*}y}\right)-f\left(x^{*},v^{*}\right)v^{*}g\left(\frac{y_{{\tau }_2}{v}^{*}}{y^{*}v}\right)d\tau \\ +p\left(y^{*}-\frac{b}{c}\right)\left(z-z^{*}\right)+p\left(y-\frac{y_{{\tau }_3}{z}_{{\tau }_3}}{z}\right)z^{*}.\end{array}$ (4.1)

From (H₂), we have:

$\left(1-\frac{f\left(x^{*},v^{*}\right)}{f\left(x,v^{*}\right)}\right)\left(x^{\mathrm{<mo>\; *\; </mo>}}-x\right)\le \mathrm{0,}$

and from (H₄) we have:

$\left(\frac{f\left(x,v^{*}\right)}{f\left(x,v\right)}-\frac{v}{v^{\mathrm{<mo>\; *\; </mo>}}}\right)\left(1-\frac{f\left(x,v\right)}{f\left(x,v^{*}\right)}\right)\le \mathrm{0,}$

and as g is positive, we have:

${\stackrel{\dot{}}{U}}_1\le {\text{e}}^{m_1{\tau }_1}p\left(y^{\mathrm{<mo>\; *\; </mo>}}-\frac{b}{c}\right)\left(z-z^{\mathrm{<mo>\; *\; </mo>}}\right)+{\text{e}}^{m_1{\tau }_1}p\left(y-\frac{y_{{\tau }_3}{z}_{{\tau }_3}}{z}\right)z^{\mathrm{<mo>\; *\; </mo>}}\mathrm{.}$

Now, if we approach $E^{\mathrm{<mo>\; *\; </mo>}}$ by $E_1^{\mathrm{<mo>\; *\; </mo>}}$ ; in particular $z^{*}=0$ and $y^{*}=y_1^{*}$ we obtain:

${\stackrel{\dot{}}{U}}_1\le {\text{e}}^{m_1{\tau }_1}p\left(y_1^{\mathrm{<mo>\; *\; </mo>}}-\frac{b}{c}\right)z\mathrm{.}$

From Remark 1, we have ${\stackrel{\dot{}}{U}}_1\left(t\right)\le 0$ for all $x\mathrm{,}y\mathrm{,}v\mathrm{,}z\ge 0$. In this case, it is easy to verify that from (4.3), the largest invariant set in $\left\{\left(x,y,v,z\right)\backslash {\stackrel{\dot{}}{U}}_1\left(t\right)=0\right\}$ is the singleton $\left\{E_1^{\mathrm{<mo>\; *\; </mo>}}\right\}$. Using LaSalle invariance principle [24], if $R_1\le 1<R_0$ then the equilibrium $E_1^{\mathrm{<mo>\; *\; </mo>}}$ is globally asymptotically stable.

And, if we approach $E^{\mathrm{<mo>\; *\; </mo>}}$ by $E_2^{\mathrm{<mo>\; *\; </mo>}}$ and we suppose ${\tau }_3=0$, in particular $y^{*}=y_2^{*}=\frac{b}{c}$ we obtain ${\stackrel{\dot{}}{U}}_1\le 0.$

Thus, the equilibrium $E_2^{\mathrm{<mo>\; *\; </mo>}}$ is stable. In this case, note that ${\stackrel{\dot{}}{U}}_1\left(t\right)=0$ if and only if $x=x_2^{*}$, $y=y_2^{*}$, and $v=v_2^{*}$ and using the first and second equation of (1.2), we obtain $z=z_2^{*}$. Therefore, it follows from LaSalle's invariance principle [24] that the CTL-activated infection equilibrium $E_2^{\mathrm{<mo>\; *\; </mo>}}$ is globally asymptotically stable. This completes the proof.

We have so obtained the results of [1] [3] [4] [8] [9] [17] on the passage of the global stability from one equilibrium to another, depending on the parameters $R_0$ and $R_1$. In what follows, we will show that stability of $E_2^{\mathrm{<mo>\; *\; </mo>}}$ depends mainly on the delay ${\tau }_3$.

5. The CTL-Activated Equilibrium and Hopf Bifurcation

In this section we will take ${\tau }_1={\tau }_2=0$ and ${\tau }_3=\tau$. The characteristic equation of system (1.2) at the CTL-activated equilibrium $E_2^{\mathrm{<mo>\; *\; </mo>}}$ is of the form:

${\lambda }^4+p_3{\lambda }^3+p_2{\lambda }^2+p_1\lambda +p_0+\left(q_3{\lambda }^3+q_2{\lambda }^2+q_1\lambda +q_0\right){\text{e}}^{-\lambda \tau }=\mathrm{0,}$ (5.1)

where:

$p_0=bu\left(\delta +p{z}_2^{*}\right)\left(d+\frac{\partial f\left(x_2^{*},v_2^{*}\right)v_2^{*}}{\partial x}\right)-dbk\left(f\left(x_2^{*},v_2^{*}\right)+\frac{\partial f\left(x_2^{*},v_2^{*}\right)v_2^{*}}{\partial v}\right),$

$\begin{array}{l}{p}_1=bu\left(\delta +p{z}_2^{*}\right)+\left(d+\frac{\partial f\left(x_2^{*},v_2^{*}\right)v_2^{*}}{\partial x}\right)\left(u\left(\delta +p{z}_2^{*}\right)+b\left(\delta +p{z}_2^{*}+u\right)\right)\\ -k\left(b+d\right)\left(f\left(x_2^{*},v_2^{*}\right)+\frac{\partial f\left(x_2^{*},v_2^{*}\right)v_2^{\mathrm{<mo>\; *\; </mo>}}}{\partial v}\right)\end{array}$

$\begin{array}{l}{p}_2=u\left(\delta +p{z}_2^{*}\right)+b\left(\delta +p{z}_2^{*}+u\right)+\left(\delta +p{z}_2^{*}+u+b\right)\left(d+\frac{\partial f\left(x_2^{*},v_2^{*}\right)v_2^{*}}{\partial x}\right)\\ -k\left(f\left(x_2^{*},v_2^{*}\right)+\frac{\partial f\left(x_2^{*},v_2^{*}\right)v_2^{*}}{\partial v}\right),\end{array}$

$p_3=\delta +p{z}_2^{*}+u+b+d+\frac{\partial f\left(x_2^{*},v_2^{*}\right)v_2^{*}}{\partial x}$

$q_0=-bu\delta \left(d+\frac{\partial f\left(x_2^{*},v_2^{*}\right)v_2^{*}}{\partial x}\right)+dbk\left(f\left(x_2^{*},v_2^{*}\right)+\frac{\partial f\left(x_2^{*},v_2^{*}\right)v_2^{*}}{\partial v}\right),$

$q_1=bk\left(f\left(x_2^{*},v_2^{*}\right)+\frac{\partial f\left(x_2^{*},v_2^{*}\right)v_2^{*}}{\partial v}\right)-bu\delta -b\left(\delta +u\right)\left(d+\frac{\partial f\left(x_2^{*},v_2^{*}\right)v_2^{*}}{\partial x}\right),$

$q_2=-b\left(\delta +u+d+\frac{\partial f\left(x_2^{*},v_2^{*}\right)v_2^{*}}{\partial x}\right),$

$q_3=-b.$ (5.2)

We have already seen in Theorem 4.1 for ${\tau }_1\ge 0$ and ${\tau }_2\ge 0$ that the equilibrium $E_2^{\mathrm{<mo>\; *\; </mo>}}$ is globally asymptotically stable in the case ${\tau }_3=\tau =0$ and in particular is locally asymptotically stable in this case.

Now, let see to what value of ${\tau }_3>0$, this stability persists.

If $iw\left(w>0\right)$ is a solution of (5.1), separating real and imaginary parts, it follows that:

$\begin{array}{l}{w}^4-p_2{w}^2+p_0=\left(q_3{w}^3-q_{1}w\right)\sin w\tau +\left(q_2{w}^2-q_0\right)\cos w\tau ,\\ -p_3{w}^3+p_{1}w=\left(q_3{w}^3-q_{1}w\right)\cos w\tau -\left(q_2{w}^2-q_0\right)\sin w\tau .\end{array}$ (5.3)

Squaring and adding the two equations of (5.3), we obtain that:

$w^8+h_3{w}^6+h_2{w}^4+h_1{w}^2+h_0=0,$

where:

$\begin{array}{l}{h}_0=p_0^2-q_0^2,\\ h_1=p_1^2-q_1^2+2q_0{q}_2-2p_0{p}_2,\\ h_2=p_2^2-q_2^2+2p_0+2q_1{q}_3-2p_1{p}_3,\\ h_3=p_3^2-q_3^2-2p_2.\end{array}$ (5.4)

Letting $z=w^2$, Equation (5.4) can be written as:

$h\left(z\right):=z^4+h_3{z}^3+h_2{z}^2+h_{1}z+h_0=0.$ (5.5)

Denote:

$P=\frac{8h_2-3h_3^2}{16},Q=\frac{h_3^3-4h_2{h}_3+8h_1}{32},D_0=\frac{Q^2}{4}+\frac{P^3}{27},$

and define:

${z^{\prime }}_1=-\frac{h_3}{4}+\sqrt[3]{-\frac{Q}{2}+\sqrt{D_0}}+\sqrt[3]{-\frac{Q}{2}-\sqrt{D_0}},\text{if}{D}_0>0,$

${z^{\prime }}_2=\max \left\{-\frac{h_3}{4}-2\sqrt[3]{\frac{Q}{2}},-\frac{h_3}{4}+\sqrt[3]{\frac{Q}{2}}\right\},\text{if}{D}_0=0,$

${z^{\prime }}_3=\max \left\{-\frac{h_3}{4}+2Re\left\{\theta \right\}\mathrm{,}-\frac{h_3}{4}+2Re\left\{\theta \epsilon \right\}\mathrm{,}-\frac{h_3}{4}+2Re\left\{\theta \bar{\epsilon }\right\}\right\}\mathrm{,}\text{if}{D}_0=\mathrm{0,}$

where $\theta$ is one of cube roots of the complex number $\frac{-Q}{2}+\sqrt{D_0}$ and $\epsilon =-\frac{1}{2}+\frac{\sqrt{3}}{2}i$.

According to [25] we have the following Lemma.

Lemma 1. [25]

For the polynomial Equation (5.6), the following states are true.

1) If $h_0\ge 0$ and one of the following conditions holds,

a) $D_0>0$ and ${z^{\prime }}_1<0$ ;

b) $D_0=0$ and ${z^{\prime }}_2<0$ ;

c) $D_0<0$ and ${z^{\prime }}_3<0$.

Then Equation (5.6) has no positive root.

2) If $h_0<0$ or $h_0\ge 0$ and one of the following conditions holds,

a) $D_0>0$, ${z^{\prime }}_1>0$ and $h\left({z^{\prime }}_1\right)<0$ ;

b) $D_0=0$, ${z^{\prime }}_2>0$ and $h\left({z^{\prime }}_2\right)<0$ ;

c) $D_0<0$, ${z^{\prime }}_3>0$ and $h\left({z^{\prime }}_3\right)<0$.

Then Equation (5.6) has at least one positive root.

Noting that:

$h_3={\left(\delta +p{z}_2^{*}\right)}^2+u^2+{\left(d+\frac{\partial f\left(x_2^{*},v_2^{*}\right)v_2^{*}}{\partial x}\right)}^2+2k\left(f\left(x_2^{*},v_2^{*}\right)+\frac{\partial f\left(x_2^{*},v_2^{*}\right)v_2^{*}}{\partial v}\right)\ge 0.$

So Equation (5.6) has at most three positive roots. Suppose that Equation (5.6) has $m$ $\left(1\le m\le 3\right)$ positive real roots $z_n$ $\left(1\le n\le m\right)$, respectively. Then Equation (5.4) has $m$ positive roots $w_n=\sqrt{z_n}$ $\left(1\le n\le m\right)$. From (5.3), we have:

$\cos \left(w_n\tau \right)=Q_n=\frac{\left(w_n^4-p_2{w}_n^2+p_0\right)\left(q_2{w}_n^2-q_0\right)-\left(p_3{w}_n^3-p_1{w}_n\right)\left(q_3{w}_n^3-q_1{w}_n\right)}{{\left(q_3{w}_n^3-q_1{w}_n\right)}^2+{\left(q_2{w}_n^2-q_0\right)}^2},$

$\sin \left(w_n\tau \right)=P_n=\frac{\left(w_n^4-p_2{w}_n^2+p_0\right)\left(q_3{w}_n^3-q_1{w}_n\right)+\left(p_3{w}_n^3-p_1{w}_n\right)\left(q_2{w}_n^2-q_0\right)}{{\left(q_3{w}_n^3-q_1{w}_n\right)}^2+{\left(q_2{w}_n^2-q_0\right)}^2},$

Let:

${\tau }_n^{\left(j\right)}=\{\begin{array}{ll}\frac{1}{w_n}\left[\arccos \left(Q_n\right)+2j\pi \right]\hfill & \text{if}{P}_n\ge 0\hfill \\ \frac{1}{w_n}\left[2\pi -\arccos \left(Q_n\right)+2j\pi \right]\hfill & \text{if}{P}_n\le 0\hfill \end{array}$ (5.6)

where $\left(1\le n\le m\right)$ and $j=0,1,2,\cdots$. Then $\pm i{w}_n$ is a pair of purely imaginary roots of Equation (5.1) with $\tau ={\tau }_n^{\left(j\right)}$.

Define:

${\tau }_0={\tau }_{n0}^{\left(0\right)}=\underset{\left(1\le n\le m\right)}{\min }\left\{{\tau }_n^{\left(0\right)}\right\},w_0=w_{n0}.$ (5.7)

Let $\lambda \left(\tau \right)=\alpha \left(\tau \right)+iw\left(\tau \right)$ be a root of Equation (5.1) satisfying $\alpha \left({\tau }_n^j\right)=0,w\left({\tau }_n^j\right)=w_n$.

Differentiating the two sides of Equation (5.1) with respect to $\tau$, we obtained

${\left(\frac{\text{d}\lambda }{\text{d}\tau }\right)}^{-1}=-\frac{4{\lambda }^3+3p_3{\lambda }^2+2p_2\lambda +p_1}{\lambda \left({\lambda }^4+p_3{\lambda }^3+p_2{\lambda }^2+p_1\lambda +p_0\right)}+\frac{3q_3{\lambda }^2+2q_2\lambda +q_1}{\lambda \left(q_3{\lambda }^3+q_2{\lambda }^2+q_1\lambda +q_0\right)}-\frac{\tau }{\lambda },$

and we have:

$\begin{array}{l}sign{\left\{\frac{d\left(Re\lambda \right)}{d\tau }\right\}}_{\tau \mathrm{<mo>\; =\; </mo>}{\tau }_n^{\left(j\right)}}=sign{\left\{Re{\left(\frac{\text{d}\lambda }{\text{d}\tau }\right)}^{-1}\right\}}_{\tau ={\tau }_n^{\left(j\right)}}\\ =sign[\frac{4w_n^6+3\left(p_3^2-2p_2\right)w_n^4+2\left(p_2^2+2p_0-2p_1{p}_3\right)w_n^2+p_1^2-2p_0{p}_2}{{\left(p_3{w}_n^2-p_1\right)}^2{w}_n^2+{\left(w_n^4-p_2{w}_n^2+p_0\right)}^2}\\ +\frac{-3q_3^2{w}_n^4+2\left(2q_1{q}_3-q_2^2\right)w_n^2+2q_2{q}_0-q_1^2}{{\left(q_3{w}_n^2-q_1\right)}^2{w}_n^2+{\left(q_0-q_2{w}_n^2\right)}^2}]\mathrm{.}\end{array}$

From (5.4), we get:

${\left(p_3{w}_n^2-p_1\right)}^2{w}_n^2+{\left(w_n^4-p_2{w}_n^2+p_0\right)}^2={\left(q_3{w}_n^2-q_1\right)}^2{w}_n^2+{\left(q_0-q_2{w}_n^2\right)}^2\mathrm{.}$

It therefore follows that:

$sign{\left\{\frac{\text{d}\left(Re\lambda \right)}{\text{d}\tau }\right\}}_{\tau ={\tau }_n^{\left(j\right)}}=sign\left[\frac{\frac{\text{d}h\left(z_n\right)}{\text{d}z}}{{\left(q_3{w}_n^2-q_1\right)}^2{w}_n^2+{\left(q_0-q_2{w}_n^2\right)}^2}\right]\mathrm{.}$

Since $z_n>0$, we conclude that $Re\left[{\frac{\text{d}{\lambda }_n\left(\tau \right)}{\text{d}\tau }|}_{\tau ={\tau }_n^{\left(j\right)}}\right]$ and $\frac{\text{d}h\left(z_n\right)}{\text{d}z}$ have the same sign.

After this discussion, we have the following theorem.

Theorem 5.1. Let ${\tau }_n^{\left(j\right)}$ and $w_0$, ${\tau }_0$ be defined by (5.10) and (5.11), respectively. If $R_1>1$ we have the following results:

1) If $h_0\ge 0$ and one of the following conditions holds,

a) $D_0>0$ and ${z^{\prime }}_1<0$ ;

b) $D_0=0$ and ${z^{\prime }}_2<0$ ;

c) $D_0<0$ and ${z^{\prime }}_3<0$.

Then the CTL-activated infection equilibrium $E_2^{\mathrm{<mo>\; *\; </mo>}}$ of system (1.2) is locally asymptotically stable for all $\tau \ge 0$.

2) If $h_0<0$ or $h_0\ge 0$ and one of the following conditions holds:

a) $D_0>0$, ${z^{\prime }}_1>0$ and $h\left({z^{\prime }}_1\right)<0$ ;

c) $D_0=0$, ${z^{\prime }}_2>0$ and $h\left({z^{\prime }}_2\right)<0$ ;

c) $D_0<0$, ${z^{\prime }}_3>0$ and $h\left({z^{\prime }}_3\right)<0$.

Then the CTL-activated infection equilibrium $E_2^{\mathrm{<mo>\; *\; </mo>}}$ of system (1.2) is locally asymptotically stable for $\tau \in \left[0,{\tau }_0\right)$.

3) System (1.2) undergoes a Hopf bifurcation at CTL-activated infection equilibrium when $\tau ={\tau }_n^{\left(j\right)}$ if the condition as stated in (2) are satisfied and:

$\frac{\text{d}h\left(z_n\right)}{\text{d}z}\ne 0$.

6. Numerical Simulations

In this section, we give some numerical simulations supporting the theoretical analysis given in Sections 3, 4 and 5. We assume $f\left(x\mathrm{,}v\right)=\frac{\beta x}{1+\alpha v}$.

Parameters in Figure 1 are $s=9$, $\beta =0.005$, $\alpha =50$, $d=0.1$, $\delta =0.36$, $k=140$, $u=200$, $c=0.1$, $b=0.2$, $p=0.1$, $m_1=0.01$, $m_2=0.02$, ${\tau }_1=1$, ${\tau }_2=1.2$, ${\tau }_3=2$ (see [26] [27] [28] [29] [30] ).

By making the calculations we have $R_0=0.9053<1$, and that system (1.2) has a infection-free equilibrium $E_0^{*}=\left(90,0,0,0\right)$. By Theorem 3.1, we get that the infection-free equilibrium $E_0^{\mathrm{<mo>\; *\; </mo>}}$ is globally asymptotically stable. Numerical simulation illustrates this fact.

Parameters in Figure 2 are $s=9$, $\beta =0.005$, $\alpha =50$, $d=0.01$, $\delta =0.36$, $k=140$, $u=2$, $c=0.1$, $b=0.2$, $p=0.1$, $m_1=0.01$, $m_2=0.02$, ${\tau }_1=1$, ${\tau }_2=1.2$, ${\tau }_3=2$. By making the calculations we have $R_0=905.26>1$ and $R_1=0.11<1$, in this case, the system (1.2) has a CTL-inactivated infection equilibrium $E_1^{*}=\left(891.09,0.24,16.72,0\right)$ who is globally asymptotically stable according to the Theorem 4.1, Therefore Numerical simulation prove this result.

Parameters in Figure 3 are $s=9$, $\beta =0.005$, $\alpha =1$, $d=0.01$, $\delta =0.36$, $p=0.1$, $k=2000$, $u=2$, $c=0.1$, $b=0.2$, $m_1=0.01$, $m_2=0.02$, ${\tau }_1=1$, ${\tau }_2=1.2$, ${\tau }_3=0$, $m_1=0.01$, $m_2=3$. By making the calculations we have $R_1=5.58>1$, and that system (1.2) has a CTL-activated infection equilibrium $E_2^{*}=\left(600.1,2,54.64,11.2\right)$. By Theorem 4.1, we obtain that the CTL-activated infection equilibrium $E_2^{*}$ of system (1.2) is globally asymptotically stable. Numerical simulation illustrates our result.

Figure 1. Global stability on the infection-free equilibrium $E_0^{\mathrm{<mo>\; *\; </mo>}}$ of system (1.2).

Figure 2. Global stability of the CTL-inactivated infection equilibrium $E_1^{\mathrm{<mo>\; *\; </mo>}}$ of system (1.2).

Figure 3. Global stability of the CTL-activated infection equilibrium $E_2^{\mathrm{<mo>\; *\; </mo>}}$ of system (1.2) in the case ${\tau }_3=0$.

Parameters in Figure 4 are $s=4.5$, $\beta =1.23$, $\alpha =0.1$, $d=2.43$, $\delta =0.01$, $p=3$, $k=2$, $u=3.2$, $c=9.33$, $b=5$, ${\tau }_1={\tau }_2=0$, ${\tau }_3=1.3>{\tau }_0=0.1756$. By making the calculations we have $R_1=137.58>1$. By Theorem 5.1, we see that the CTL-activated infection equilibrium $E_2^{\mathrm{<mo>\; *\; </mo>}}$ is locally asymptotically stable if $0<{\tau }_3<{\tau }_0$ and unstable if $\tau >{\tau }_0$. Further, system 1.2 undergoes a Hopf bifurcation at the CTL-activated infection equilibrium when $\tau ={\tau }_0$. Numerical simulation illustrates our result.

7. Examples

In this section, we give some particular examples.

Example 1. In (1.2), if $f\left(x,v\right)=\beta x$ and ${\tau }_2={\tau }_3=0$, we obtain the model studied by Zhu et al. [1].

Example 2. Let $f\left(x,v\right)=\frac{\beta x}{1+\alpha v}$ and ${\tau }_3=m_1=m_2=0$, in this case we obtain the model presented by Li et al. [4]. So the work presented in [4] is a particular case of (1.2) because the function $\frac{\beta x}{1+\alpha v}$ satisfies the hypotheses (H₁) - (H₄).

Figure 4. Stability and Hopf bifurcation of the CTL-activated infection equilibrium $E_2^{\mathrm{<mo>\; *\; </mo>}}$ of system (1.2) in the case, ${\tau }_1={\tau }_2=0$ and ${\tau }_3\ne 0$.

Example 3. In (1.2), if $f\left(x,v\right)=\frac{\beta x}{1+{\alpha }_{1}x+{\alpha }_{2}v}$ and ${\tau }_1={\tau }_2={\tau }_3=0$, we obtain a virus dynamics model with Beddington-Deangelis incidence rate and CTL immune response. This model is presented in [17]. The function $\frac{\beta x}{1+{\alpha }_{1}x+{\alpha }_{2}v}$

satisfies the hypotheses (H₁) - (H₄), so the model presented in [17] is a particular case of (1.2).

Example 4. Let $f\left(x,v\right)=\frac{\beta x}{\left(1+{\alpha }_{1}x\right)\left(1+{\alpha }_{2}v\right)}$ and ${\tau }_1={\tau }_2={\tau }_3=0$. The

hypotheses (H₁) - (H₄) are satisfied. In this case we obtain a virus dynamics model with Crowly-Martin functional response. The global properties of this model were studied in [21].

Example 5. A last example, in (1.2), if $f\left(x,v\right)=\frac{\beta x}{1+\alpha v}$ and ${\tau }_1={\tau }_2=0$, we obtain the results presented in [13].

8. Conclusion

In this paper, we have considered a virus dynamics model with a general incidence rate and three delays ${\tau }_1,{\tau }_2$ and ${\tau }_3$. This general incidence represents a variety of possible incidence functions that could be used in virus dynamics models. We establish that the global dynamics are determined by two threshold parameters the basic reproduction ratios for viral infection and CTL immune response $R_0$ and $R_1$, respectively. We have proved that the infection-free equilibrium $E_0^{\mathrm{<mo>\; *\; </mo>}}$ is globally asymptotically stable if the basic reproduction ratios viral infection $R_0\le 1$ for all ${\tau }_1\ge 0$, ${\tau }_2\ge 0$ and ${\tau }_3\ge 0$. The hypotheses on the general incidence function are used to assure the existence of the CTL-inactivated infection equilibrium $E_1^{\mathrm{<mo>\; *\; </mo>}}$ and the CTL-activated infection equilibrium $E_2^{\mathrm{<mo>\; *\; </mo>}}$. We prove that if $R_1\le 1<R_0$, the CTL-inactivated infection equilibrium $E_1^{\mathrm{<mo>\; *\; </mo>}}$ is globally asymptotically stable for all ${\tau }_1\mathrm{,}{\tau }_2\mathrm{,}{\tau }_3\ge 0$. $R_1>1$, assure the existence of $E_2^{\mathrm{<mo>\; *\; </mo>}}$ and when ${\tau }_3=0$, $E_2^{\mathrm{<mo>\; *\; </mo>}}$ is globally asymptotically stable for all ${\tau }_1\mathrm{,}{\tau }_2\ge 0$, and when ${\tau }_3>0$ and ${\tau }_1={\tau }_2=0$, $E_2^{\mathrm{<mo>\; *\; </mo>}}$ is only locally asymptotically stable for ${\tau }_3$ less than a critical value ${\tau }_3^{\mathrm{<mo>\; *\; </mo>}}$. Moreover, in ${\tau }_3^{\mathrm{<mo>\; *\; </mo>}}$, there is the appearance of a periodic phenomenon from $E_2^{\mathrm{<mo>\; *\; </mo>}}$, obtained by Hopf bifurcation and when ${\tau }_3$ passes ${\tau }_3^{\mathrm{<mo>\; *\; </mo>}}$, $E_2^{\mathrm{<mo>\; *\; </mo>}}$ becomes unstable. In the end, we have presented some numerical simulations.

Acknowledgements

We would like to thank the anonymous referees for their very helpful suggestions and comments.

Conflicts of Interest

The author declares that there is no conflict of interests regarding the publication of this paper.

References

[1]	Zhu, H. and Zou, X. (2009) Dynamics of a HIV-1 Infection Model with Cell-Mediated Immune Response and Intracellular Delay. Discrete & Continuous Dynamical Systems-B, 12, 511-524. https://doi.org/10.3934/dcdsb.2009.12.511
[2]	Song, X. and Neumann, A.U. (2007) Global Stability and Periodic Solution of the Viral Dynamics. Journal of Mathematical Analysis and Applications, 329, 281-297. https://doi.org/10.1016/j.jmaa.2006.06.064
[3]	Li, D. and Ma, W. (2007) Asymptotic Properties of a HIV-1 Infection Model with Time Delay. Journal of Mathematical Analysis and Applications, 335, 683-691. https://doi.org/10.1016/j.jmaa.2007.02.006
[4]	Li, Y., Xu, R., Li, Z. and Mao, S. (2011) Global Dynamics of a Delayed HIV-1 Infection Model with CTL Immune Response. Discrete Dynamics in Nature and Society, 2011, Article ID: 673843. https://doi.org/10.1155/2011/673843
[5]	Wang, Z. and Xu, R. (2012) Stability and Hopf Bifurcation in a Viral Infection Model with Nonlinear Incidence Rate and Delayed Immune Response. Communications in Nonlinear Science and Numerical Simulation, 17, 964-978. https://doi.org/10.1016/j.cnsns.2011.06.024
[6]	Cai, L.M., Guo, B.Z. and Li, X.Z. (2012) Global Stability for a Delayed HIV-1 Infection Model with Nonlinear Incidence of Infection. Applied Mathematics and Computation, 219, 617-623. https://doi.org/10.1016/j.amc.2012.06.051
[7]	Perelson, A.S., Neumann, A.U., Markowitz, M., Leonard, J.M. and Ho, D.D. (1996) HIV-1 Dynamics in Vivo: Virion Clearance Rate, Infected Cell Life-Span, and Viral Generation Time. Science, 271, 1582-1586. https://doi.org/10.1126/science.271.5255.1582
[8]	Xu, R. (2011) Global Stability of an HIV-1 Infection Model with Saturation Infection and Intracellular Delay. Journal of Mathematical Analysis and Applications, 375, 75-81. https://doi.org/10.1016/j.jmaa.2010.08.055
[9]	Nakata, Y. (2011) Global Dynamics of a Viral Infection Model with a Latent Period and Beddington-DeAngelis Response. Nonlinear Analysis: Theory, Methods and Applications, 74, 2929-2940. https://doi.org/10.1016/j.na.2010.12.030
[10]	Yuan, Z., Ma, Z. and Tang, X. (2012) Global Stability of a Delayed HIV Infection Model with Nonlinear Incidence Rate. Nonlinear Dynamics, 68, 207-214. https://doi.org/10.1007/s11071-011-0219-8
[11]	Nowak, M.A., Bonhoeffer, S., Shaw, G.M. and May, R.M. (1997) Anti-Viral Drug Treatment: Dynamics of Resistance in Free Virus and Infected Cell Populations. Journal of Theoretical Biology, 184, 203-217. https://doi.org/10.1006/jtbi.1996.0307
[12]	Perelson, A.S. and Nelson, P.W. (1999) Mathematical Analysis of HIV-1 Dynamics in Vivo. SIAM Review, 41, 3-44. https://doi.org/10.1137/S0036144598335107
[13]	Tian, X. and Xu, R. (2014) Global Stability and Hopf Bifurcation of an HIV-1 Infection Model with Saturation Incidence and Delayed CTL Immune Response. Applied Mathematics and Computation, 237, 146-154. https://doi.org/10.1016/j.amc.2014.03.091
[14]	Ndongo, A.S. and Talibi Alaoui, H. (2014) Global Stability of a HIV-1 Infection Model with General Incidence Rate and Distributed Delays. International Scholarly Research Notices, 2014, Article ID: 260379. https://doi.org/10.1155/2014/260379
[15]	Wang, W., Liu, M. and Zhao, J. (2013) Analysis of a Delayed Sir Model with Exponential Birth and Saturated Incidence Rate. Applied Mathematics, 4, 60-67. https://doi.org/10.4236/am.2013.410A2006
[16]	Nowak, M.A. and Bangham, C.R. (1996) Population Dynamics of Immune Responses to Persistent Viruses. Science, 272, 74-79. https://doi.org/10.1126/science.272.5258.74
[17]	Wang, X., Tao, Y. and Song, X. (2011) Global Stability of a Virus Dynamics Model with Beddington-DeAngelis Incidence Rate and CTL Immune Response. Nonlinear Dynamics, 66, 825-830. https://doi.org/10.1007/s11071-011-9954-0
[18]	Beddington, J. (1975) Mutual Interference between Parasites or Predators and Its Effect on Searching Efficiency. Journal of Animal Ecology, 44, 331-340. https://doi.org/10.2307/3866
[19]	DeAngelis, D.L., Goldstein, R. and Oneill, R. (1975) A Model for Tropic Interaction. Ecology, 56, 881-892. https://doi.org/10.2307/1936298
[20]	Crowley, P.H. and Martin, E.K. (1989) Functional Responses and Interference within and between Year Classes of a Dragonfly Population. Journal of the North American Benthological Society, 8, 211-221. https://doi.org/10.2307/1467324
[21]	Xu, S. (2012) Global Stability of the Virus Dynamics Model with Crowley-Martin Functional Response. Electronic Journal of Qualitative Theory of Differential Equations, No. 9, 1-10. https://doi.org/10.14232/ejqtde.2012.1.9
[22]	Wang, K., Wang, W., Pang, H. and Liu, X. (2007) Complex Dynamic Behavior in a Viral Model with Delayed Immune Response. Physica D: Nonlinear Phenomena, 226, 197-208. https://doi.org/10.1016/j.physd.2006.12.001
[23]	Yang, X., Chen, L. and Chen, J. (1996) Permanence and Positive Periodic Solution for the Single-Species Nonautonomous Delay Diffusive Models. Computers & Mathematics with Applications, 32, 109-116. https://doi.org/10.1016/0898-1221(96)00129-0
[24]	Lasalle, J. (1976) The Stability of Dynamical Systems Regional Conference Series in Applied Mathematics. SIAM, Philadelphia.
[25]	Yan, X.P. and Li, W.T. (2006) Stability and Bifurcation in a Simplified Four-Neuron BAM Neural Network with Multiple Delays. Discrete Dynamics in Nature and Society, 2006, 032529. https://doi.org/10.1155/DDNS/2006/32529
[26]	Culshaw, R. and Ruan, S. (2000) A Delay-Differential Equation Model of HIV Infection of CD4+ T-Cells. Mathematical Biosciences, 165, 27-39. https://doi.org/10.1016/S0025-5564(00)00006-7
[27]	Kirschner, D. and Webb, G. (1996) A Model of Treatment Strategy in the Chemotherapy of AIDS. Bulletin of Mathematical Biology, 58, 367-390. https://doi.org/10.1007/BF02458312
[28]	Leenheer, P.D. and Smith, H. (2003) Virus Dynamics: A Global Analysis. SIAM Journal on Applied Mathematics, 63, 1313-1327. https://doi.org/10.1137/S0036139902406905
[29]	Dixit, N. and Perelson, A. (2004) Complex Patterns of Viral Load Decay under Antiretroviral Therapy: Influence of Pharmacokinetics and Intracellular Delay. Journal of Theoretical Biology, 226, 95-109. https://doi.org/10.1016/j.jtbi.2003.09.002
[30]	Wang, Y., Zhou, Y., Wu, J. and Heffernan, J. (2009) Oscillatory Viral Dynamics in a Delayed HIV Pathogenesis Model. Mathematical Biosciences, 219, 104-112. https://doi.org/10.1016/j.mbs.2009.03.003