Global stability of an HIV infection model with saturated CTL immune response and intracellular delay

1.   Introduction

Mathematical modeling and analysis of HIV pathogenesis are crucial in understanding the transmission mechanism of HIV. Many authors are interested in studying the global stability of HIV infection models (see, for example, ^[1,2,3,4]). In addition, the analysis of mathematical models contributes to developing new antiviral drugs and designing effective therapies for HIV infection (see, for example, ^[5,6,7]). Since cytotoxic T lymphocytes (CTLs) play a significant role in the within-host anti-HIV defense by attacking infected cells in HIV infection, in 1996, Nowak and Bangham ^[8] proposed the following HIV model describing the interaction between host cells, free viruses and CTLs:

Here, $x(t)$, $y(t)$, $v(t)$ and $z(t)$ are the concentrations of uninfected cells, infected cells, free viruses and CTLs at time $t$, respectively. The parameter $\lambda$ denotes the production rate of uninfected cells. Uninfected cells, infected cells, free viruses and CTLs die at rate $dx(t)$, $ay(t)$, $uv(t)$ and $bz(t)$, respectively. The term $\beta x(t)v(t)$ is the rate at which uninfected cells are infected by free viruses, and the term $py(t)z(t)$ represents the rate for infected cells to be killed by CTLs. Free viruses are released from infected cells at rate $ky(t)$ and CTL immune response is activated by infected cells at rate $cy(t)z(t)$.

We note that system (1.1) assumes that CTL immune response is activated at bilinear rate. However, in ^[9], De Boer suggested that bilinear rate cannot model several immune responses that are together controlling a chronic infection, and proposed immune response functions with saturation effect. In ^[10], Jiang and Wang considered saturated immune response function $cy(t)z(t)/(h+z(t))$ to replace the bilinear rate, here $h$ is a saturation constant.

It is assumed in system (1.1) that as soon as the free viruses enter a target cell, the target cell is immediately infected and new free viruses are produced simultaneously. In ^[11], a mathematical model with the intracellular phase of the viral life-cycle was first proposed by Herz et al. There is a fixed time delay $\tau$ between infection of a cell and production of new free viruses. Some researchers have incorporated the intracellular delay into HIV infection models and investigated the effect of the intracellular delay on HIV infection dynamics (see, for example, ^{[12,13,14,15]}).

On the other hand, system (1.1) assumes that the infection rate of virus-to-cell is bilinear. In fact, however, in ^[16], Ebert et al. observed a nonlinear relationship between parasite dose and infection rate in experiments. Furthermore, in ^[17], Regoes et al. showed that the infection rate is a sigmoidal function. In ^[18], Song and Neumann considered the infection rate with saturation effect, $\beta xv^{p}/(1+\alpha v^{q})$, where $p$, $q$ and $\alpha$ are positive constants, and investigated the global stability of the viral model with saturated infection rate $\beta xv/(1+\alpha v)$.

Inspired by the biological reasons mentioned above, in this paper, we study the effect of saturated infection rate, intracellular delay and saturated CTL immune response on the dynamics of HIV infection. To this end, we consider the following delay differential equations:

Here, $\tau$ denotes the time between viral entry into a cell and production of free viruses; $e^{-m\tau}$ is the probability of cellular survival from time $t-\tau$ to time $t$. All parameters of system (1.2) are positive constants.

We assume that the initial condition of system (1.2) satisfies:

where $(\phi_1(\theta), \phi_2(\theta))\in \textbf{C}([-\tau, 0], \mathbb{R}_{+0}^{2})$, $\mathbb{R}_{+0}^{2} = \{(x_1, x_2):x_i\geq0, i = 1, 2\}.$

According to the basic theory of functional differential equations ^[19], we can prove that system (1.2) has a unique solution $(x(t), y(t), v(t), z(t))$ that satisfies the initial condition (1.3). Furthermore, it is easy to show that under the initial condition (1.3), all solutions of system (1.2) are defined on $[0, +\infty)$ and are positive for all $t\geq0$.

The main goal of this paper is to make a complete mathematical analysis of system (1.2) and investigate its dynamical behaviors. This paper is organized as follows. In the next section, we calculate the reproduction numbers of system (1.2) and discuss the existence of feasible equilibria. In Section 3, we study the local asymptotic stability of each of feasible equilibria. The global asymptotic stability of each of feasible equilibria is investigated in Section 4. In Section 5, we perform sensitivity analyses to illustrate the effects of parameters on the immunity-activated reproduction number $\mathscr{R}_{1}$ and the viral load.

2.   Equilibria and reproduction numbers

It is clear that system (1.2) always admits an infection-free equilibrium $E_0(\lambda/d, 0, 0, 0)$.

Now, we calculate the immunity-inactivated reproduction number of system (1.2). According to the method of next generation matrix proposed by van den Driessche and Watmough ^[20], we have

Sequentially, we obtain that

Hence the next generation matrix is

We therefore derive the immunity-inactivated reproduction number

representing the number of second infected cells produced by a infective cell in a whole susceptible population ^[20]. It is easy to prove that if $\mathscr{R}_0 > 1$, system (1.2) has a unique immunity-inactivated equilibrium $E_1(x_1, y_1, v_1, 0)$, here

Denote

where $\mathscr{R}_1$ is called immunity-activated reproduction number of system (1.2). Besides, we can show that if $\mathscr{R}_1 > 1$, system (1.2) has an immunity-activated equilibrium $E_2 (x_2, y_2, v_2, z_2)$, where

and

here,

3.   Local asymptotic stability

In this section, we study the local dynamics of system (1.2).

Theorem 3.1. If $\mathscr{R}_0 < 1$, the infection-free equilibrium $E_0 (\lambda/d, 0, 0, 0)$ of system (1.2) is locally asymptotically stable; if $\mathscr{R}_0 > 1$, $E_0$ is unstable.

Proof. By calculation, we have the following characteristic equation of system (1.2) at $E_0$:

Clearly, Eq. (3.1) has negative real roots $\xi = -b$ and $\xi = -d$, and other roots are governed by the following equation:

Substituting $\mathscr{R}_0$ into Eq. (3.2), we can rewrite (3.2) as

Next, we verify that all roots of Eq. (3.3) have negative real parts. By contradiction, let $\xi_1 =$Re$\xi_1+\text{i}$Im$\xi_1$ with Re$\xi_1\geq0$ be a root of Eq. (3.3). It follows that

Further, we have

which contradicts (3.3). Hence, we conclude that if $\mathscr{R}_0 < 1$, all roots of Eq. (3.1) have negative real parts. Accordingly, $E_0$ is locally asymptotically stable. On the other hand, if $\mathscr{R}_0 > 1$, we denote

It is easy to see that $G(0) = au(1-\mathscr{R}_0) < 0$ and $G(\xi)\rightarrow +\infty$ as $\xi\rightarrow +\infty$. Noting that $G(\xi)$ is a continuous function in respect to $\xi$, it follows that if $\mathscr{R}_0 > 1$, Eq. (3.1) has a positive real root. Thus, $E_0$ is unstable.

Theorem 3.2. If $\mathscr{R}_1 < 1 < \mathscr{R}_0$, the immunity-inactivated equilibrium $E_1 (x_1, y_1, v_1, 0)$ of system (1.2) is locally asymptotically stable.

Proof. By calculation, we have the following characteristic equation of system (1.2) at $E_1$:

Since

Eq. (3.5) has a negative real root $\xi = \frac{cy_1}{h}-b$, and the remaining roots are determined by the following equation:

We claim that all roots of Eq. (3.7) have negative real parts. Assume the contrary, Eq. (3.7) has a root $\xi_2 =$Re$\xi_2+\text{i}$Im$\xi_2$ with Re$\xi_2\geq0$. It is easy to show that

Furthermore, we have

which contradicts (3.7). Consequently, if $\mathscr{R}_1 < 1 < \mathscr{R}_0$, all roots of Eq. (3.5) have negative real parts. Accordingly, $E_1$ is locally asymptotically stable.

Theorem 3.3. If $\mathscr{R}_1 > 1$, the immunity-activated equilibrium $E_2 (x_2, y_2, v_2, z_2)$ of system (1.2) is locally asymptotically stable.

Proof. The characteristic equation of system (1.2) at $E_2$ is:

In the following, we show that all roots of Eq. (3.9) have negative real parts. If Eq. (3.9) has a root $\xi_3$ with nonnegative real part, we obtain that

From the expression of $E_2$, we derive that $u(a+pz_2) = \frac{\beta kx_2e^{-m\tau}}{1+\alpha v_2}$. Thus, we have

This leads to a contradiction. Therefore, all roots of Eq. (3.9) have negative real parts if $\mathscr{R}_1 > 1$. Accordingly, $E_2$ is locally asymptotically stable.

4.   Global asymptotic stability

In this section, we establish the global dynamics of system (1.2).

Now, we discuss the boundedness of solutions, and have the following result.

Lemma 4.1. All solutions of system (1.2) with initial condition (1.3) are ultimately bounded for all $t\geq0$.

Proof. Let $(x(t), y(t), v(t), z(t))$ be any solution of system (1.2) satisfying the initial condition (1.3). Denote $B(t) = x(t-\tau)+e^{m\tau}y(t)$.

Differentiating $B(t)$ with respect to $t$, we obtain that

which implies that $\limsup_{t\to\infty}B(t)\leq\frac{\lambda}{\min\{a, d\}}.$ Then $x(t)$ and $y(t)$ are ultimately bounded. From the third and the fourth equations of system (1.2), we derive that

and thus

Hence $v(t)$ and $z(t)$ are ultimately bounded. Therefore, the positively invariant set for system (1.2) is given by:

Next, we define a function $g(x) = x-1-\ln x$, where $g(1) = 0$ and $g(x)$ attains its minimum at $x = 1$.

Theorem 4.2. If $\mathscr{R}_0 < 1$, the infection-free equilibrium $E_0 (\lambda/d, 0, 0, 0)$ of system (1.2) is globally asymptotically stable.

Proof. Let $(x(t), y(t), v(t), z(t))$ be any positive solution of system (1.2) satisfying the initial condition (1.3). Define

where $x_0 = \lambda/d$. Calculating the time derivative of $V_1(t)$ along positive solutions of system (1.2), we have

Since $\mathscr{R}_0 < 1$, we have $\dot{V_1}(t)\leq0$, and $\dot{V_1}(t) = 0$ if and only if $x = x_0$, $y = v = z = 0$. It is easy to show that $M_{0} = \{E_{0}\} \subset {\bf{\Omega}}$ is the largest invariant subset of $\{(x(t), y(t), v(t), z(t)) : \dot{V_{1}}(t) = 0\}$. In addition, from Theorem 3.1, we know that if $\mathscr{R}_{0} < 1$, $E_{0}$ is locally asymptotically stable. Thus, by LaSalle's invariance principle ^[19], we conclude that $E_0$ is globally asymptotically stable.

Theorem 4.3. If $\mathscr{R}_1 < 1 < \mathscr{R}_0$, the immunity-inactivated equilibrium $E_1(x_1, y_1, v_1, 0)$ of system (1.2) is globally asymptotically stable.

Proof. Let $(x(t), y(t), v(t), z(t))$ be any positive solution of system (1.2) satisfying the initial condition (1.3). Define

Differentiating $V_2(t)$ along positive solutions of system (1.2), we obtain that

Combining $\lambda = dx_1+\frac{\beta x_1v_1}{1+\alpha v_1}$, $\frac{\beta e^{-m\tau}x_1v_1}{1+\alpha v_1} = ay_1$ and $ky_1 = uv_1$ with (4.2), we have

From (3.6), we have $\frac{cy_1}{h+z(t)}-b < \frac{cy_1}{h}-b < 0$. Thus $\dot{V_2}(t)\leq0$, and $\dot{V_2}(t) = 0$ if and only if $x = x_1$, $y = y_1$, $v = v_1$, $z = 0$. It is easy to see that the largest invariant subset of $\{(x(t), y(t), v(t), z(t)) : \dot{V_{2}}(t) = 0\}$ is the singleton $M_{1} = \{E_{1}\} \subset {\bf{\Omega}}$. Furthermore, from Theorem 3.2, we know that if $\mathscr{R}_1 < 1 < \mathscr{R}_0$, $E_{1}$ is locally asymptotically stable. According to LaSalle's invariance principle ^[19], $E_1$ is globally asymptotically stable.

Theorem 4.4. If $\mathscr{R}_1 > 1$, the immunity-activated equilibrium $E_2(x_2, y_2, v_2, z_2)$ of system (1.2) is globally asymptotically stable.

Proof. Let $(x(t), y(t), v(t), z(t))$ be any positive solution of system (1.2) satisfying the initial condition (1.3). Define

Calculating the derivative of $V_3(t)$ along positive solutions of system (1.2), we obtain that

Substituting $\lambda = dx_2+\frac{\beta x_2v_2}{1+\alpha v_2}$, $\frac{\beta e^{-m\tau}x_2v_2}{1+\alpha v_2} = ay_2+py_{2}z_{2}$, $ky_2 = uv_2$ and $\frac{cy_{2}z_{2}}{h+z_{2}} = bz_{2}$ into (4.4), we get

Obviously, $\dot{V_3}(t)\leq0$ and $\dot{V_3}(t) = 0$ if and only if $x = x_2$, $y = y_2$, $v = v_2$, $z = z_2$. $M_{2} = \{E_{2}\} \subset {\bf{\Omega}}$ is the largest invariant subset of $\{(x(t), y(t), v(t), z(t)) : \dot{V_{3}}(t) = 0\}$. From Theorem 3.3, if $\mathscr{R}_1 > 1$, $E_{2}$ is locally asymptotically stable. Hence, using LaSalle's invariance principle ^[19], we obtain that $E_2$ is globally asymptotically stable.

5.   Sensitivity analyses

In this section, we perform sensitivity analyses to show the effects of parameters $\beta$, $k$, $u$ and $c$ on the immunity-activated reproduction number $\mathscr{R}_{1}$ and the viral load.

We choose the parameter values as follows ^[21,22,23]:

Meanwhile, we assume $\alpha = 0.15 \ \ \text{ml virion}^{-1}$ and $h = 0.01\ \ \text{cells ml}^{-1}$.

First, using the methods of Latin Hypercube Sampling and Partial Rank Correlation Coefficients (PRCCs) developed in ^[24], we investigate the sensitivity of the immunity-activated reproduction number $\mathscr{R}_{1}$ on the parameters $\beta$, $k$, $u$ and $c$, which have important effects on HIV infection (see Figure 1). From Figure 1, we see that $\beta$, $k$ and $c$ are positively correlated with $\mathscr{R}_{1}$ while $u$ are negatively correlated with $\mathscr{R}_{1}$. In addition, we obtain that decreasing the activation rate of CTL immune response and the infection rate of virus-to-cell infection can more effectively reduce $\mathscr{R}_{1}$.

Second, we study the sensitivity of the viral load with respect to the parameters $\beta$, $k$, $u$ and $c$ (see Figure 2). In Figure 2, it is easy to see that $\beta$, $k$ and $c$ are positively correlated with the viral load while $u$ are negatively correlated with the viral load. Besides, we get that decreasing the release rate of free viruses or increasing the death rate of free viruses can more effectively reduce the viral load, which is helpful in controling the viral infection.

Acknowledgments

The authors wish to thank the Editor and the reviewers for their valuable comments that greatly improved the presentation of this work.

This work was supported by the National Natural Science Foundation of China (Nos. 11871316, 11801340), and the Natural Science Foundation of Shanxi Province (Nos. 201801D121006, 201801D221007).

Conflict of interest

The authors have no conflicts of interest in this paper.