Viral dynamics of an HIV stochastic model with cell-to-cell infection, CTL immune response and distributed delays

1.   Introduction

Acquired immune deficiency syndrome (AIDS) is caused by human immunodeficiency virus (HIV), which is a serious threat to human health. HIV infects the human body by infecting healthy target T-cells. Meanwhile, the cellular immune response mediated by cytotoxic T lymphocytes (CTLs) can kill some infected T-cells, thus inhibiting further replication of the virus. Hence, CTLs play a significant role in the suppression of HIV by killing viral infected T-cells ^[1]. More and more scholars pay attention to the research of HIV infection modelling. The mathematical models have been revealed as a powerful tool for understanding the mechanism of HIV infection.

Most of these earlier models have focused on the interaction between virus and target cells based on a hypothesis that infected T-cells produce new virus particles immediately (that is virus-to-cell infection) ^{[2,3,4,5,6,7,8]}. However, recent studies have demonstrated that cell-to-cell transmission is largely unaffected by some obstacles compared with virus-to-cell infection, so that the cell-to-cell transmission is more efficient than virus-to-cell infection ^[9,10,11]. Inspired by the experimental data, Sigal et al. have showed that under the action of antiretroviral drugs, the virus infection caused by virus-to-cell was significantly reduced, while the drug sensitivity of infection involving cell-to-cell transmission was significantly reduced ^[12]. Cell-to-cell infection may adversely affect the immune system, leading to the persistence of the virus, thus becoming an obstacle to the treatment of HIV infection ^[12]. These studies suggest that cell-to-cell transmission contributes to the pathogenesis of viral infection. Consequently, in the process of HIV infection, we should not ignore the mode of cell-to-cell infection, which may play an important role in viral spread in vivo. Therefore, the HIV model with both virus-to-cell and cell-to-cell infection modes is of research value and significance, which we will concentrate on.

The first mathematical model involving cell-to-cell infection was proposed by Culshaw et al. ^[13], and they also considered the intracellular delay caused by cell-to-cell infection. They studied the effect of time delay on the stability of positive equilibrium, and the model exhibit Hopf bifurcation with intracellular delay as the bifurcation parameter. In fact, there are two main types of delays in the process of HIV infection within-host: (ⅰ) the intracellular delay, that is, the time it takes for a virus to infect healthy T-cell to become productively infected T-cell; (ⅱ) the immune delay, that is, the time when viral infection activates the CTL immune response. In this paper, we incorporate two routes of infection modes (both virus-to-cell and cell-to-cell) and three time delays (two intracellular delays caused by two infection modes and one immune response delay). Many scholars have examined the effects of two infection modes and multiple delays on viral dynamics, see for example ^{[14,15,16,17,18,19,20,21,22,23,24]}, and references therein.

The parameters of the previous references are considered to be a fixed constant in the average sense, while randomness is an inevitable factor in real life. In the macroscopical field of infectious diseases, many scholars have introduced the stochastic fluctuation into the process of mathematical modelling and have examined the effect of the stochastic perturbation on model behaviors ^{[25,26,27,28,29]}. In the microscopic field of HIV infection in vivo, it has been proved that HIV transcription is an inherent random process and produces strong fluctuations in virus gene expression ^[30]. Thus, random-generated expression variability is increasingly considered to have important phenotypic consequences in different cellular environments, such as multicellular development, cancer progression, and viral latency ^[31]. Mao et al. have further demonstrated that even a small random disturbance could suppress population explosion through rigorous mathematical analysis ^[32]. Hence, stochastic perturbation can be included in the process of modelling to accurately depict reality.

At present, some kinds of stochastic HIV infection model have been studied ^{[33,34,35,36]}, but these models only considered the factor of virus-to-cell infection mode, and did not involve the factors of cell-to-cell infection mode and time delays. Lately, some researchers ^[37,38] have investigated the asymptotic behaviors of a two-dimensional cell-to-cell HIV model with random noise, while they did not refer to the virus-to-cell mode. In this paper, we extend the deterministic model with virus-to-cell infection, cell-to-cell infection, CTL immune response and distributed delays by including the random fluctuations. As far as we know, few people have studied the random HIV model of virus-to-cell infection, cell-to-cell transmission and time delays.

Some authors have studied stochastic differential equations with time delays for epidemic infectious disease ^[39], Lotka-Volterra system ^[40] and algal bloom ^[41]. However, the introduction of time delay into stochastic viral dynamics model is rare. Here, we apply stochastic delay differential equations to the field of HIV within-host, and the main purpose of this paper are: (ⅰ) study the existence of stationary Markov process of a degenerate stochastic differential equations; (ⅱ) investigate the influences of noise intensity, cell-to-cell infection and time delays on virus dynamics under realistic parameter values.

The organization of this article is as follows. In Section 2, a degenerate stochastic HIV model with cell-to-cell infection and CTL immune response is derived, and the existence and uniqueness of the global positive solution are also shown. In Section 3, by formulating appropriate Lyapunov functions, we obtain the existence of a stationary Markov process. The extinction of the virus is given in Section 4. In Section 5, we take numerical simulations to verify our theoretical analysis results based on realistic parameter values of HIV in published references, and we also investigate the effects of noise intensity, cell-to-cell infection and delays on virus dynamics, respectively. Finally, we conclude our work.

2.   Model and preliminaries

According to the existing literature ^{[14,18,19,20,21,22,23,24]}, a deterministic HIV infection model including cell-to-cell infection, CTL immune response and distribution delays is as follows

where, $T(t)$, $I(t)$, $V(t)$ and $E(t)$ represent the concentrations of healthy T-cells, infected cells, virions and CTLs at time $t$, respectively. Parameter $\lambda$ is the source of CD4$^+$ T-cells from precursors. The mitosis of healthy T-cells is described as the logistic term $rT(t)\big{(}1-\frac{T(t)}{T_{max}}\big{)}$, where $r$ is the intrinsic mitosis rate and $T_{max}$ is the carrying capacity of the healthy T-cells. $\mu_i$ $(i = 1, 2, 3, 4)$ are the death rates of $T(t)$, $I(t)$, $V(t)$ and $E(t)$ populations, respectively. $\beta_1$ is the infection rate of free virus by virus-to-cell infection mode, and $\beta_2$ is the infection rate of productively infected cells by cell-to-cell infection mode. The probability distribution functions $f_1(\tau)$ and $f_2(\tau)$ stand for the time for infected T-cells to become productively infected due to virus-to-cell infection and cell-to-cell infection modes, respectively. $e^{-s_1\tau}$ and $e^{-s_2\tau}$ are the survival rates of cells that are infected by virus and infected cells at time $t$ and become activated infected $\tau$ time. The delay of the mature viral particles is described by the probability distribution $f_3(\tau)$ and $e^{-s_3\tau}$ accounts for the survival probability during the delay period $\tau$. $k$ represents the average production rate of virus from an infected T-cell. $q$ is CTL effectiveness and $p$ is CTL responsiveness.

We assume that $f_i(\tau): [0, \infty)\rightarrow [0, \infty)$ are probability distributions with compact support, $f_i(\tau)\geq 0$ and $\int_0^\infty f_i(\tau)d\tau = 1$, $i = 1, 2, 3$. For the distributed delays, the kernels are usually been chosen as a gamma distribution ^[42,43],

where $n$ is a nonnegative integer. For convenience of this study, we take all the kernels as weak kernels case, that is the gamma distribution with $n = 0$,

For system (2.1) with weak kernels (2.2), these authors ^{[14,19,22,23,24]} have studied its dynamics theoretically. It is shown that system (2.1) with weak kernels (2.2) always has an infection-free equilibrium $E_0$ $(T_0, 0, 0, 0)$, where

The basic reproduction number is,

where, basic reproduction number $\mathcal{R}_{01}$ stands for the infection by virus-to-cell infection mode, and $\mathcal{R}_{02}$ stands for the infection by cell-to-cell infection mode. Summarizing the results of references ^{[14,19,22,23,24]}, the main theoretical results of system (2.1) with weak kernels (2.2) are as follows:

(Ⅰ) If $\mathcal{R}_0 < 1$, the infection-free equilibrium $E_0$ is globally asymptotically stable under the condition $s_1\tau_1 = s_2\tau_2$.

(Ⅱ) If $\mathcal{R}_0 > 1$, the positive equilibrium $E^*$ $(T^*, I^*, V^*, E^*)$ is globally attractive under the conditions $s_1\tau_1 = s_2\tau_2$ and $r\left(1-\frac{T^*}{T_{max}}\right) < \mu_1$.

In this paper, considering random fluctuation of system (2.1), we assume that the stochastic fluctuation is the white noise type, that is

System (2.1) with random fluctuations can be written as the following stochastic integro-delay differential system

where, $B_i(t)$ $(1\leq i\leq4)$ are independent standard Brownian motions with $B_i(0) = 0$, and $\sigma_i^2 > 0$ $(1\leq i\leq4)$ represent the intensities of the white noises. The remaining parameters meanings are the same as in system (2.1).

In the following, we mainly focus on the weak kernels (2.2) case for system (2.4). Let

By calculation, we derive that

Therefore, system (2.4) with weak kernels (2.2) can be rewritten as the following degenerate stochastic differential system

Throughout this article, let $B_i(t)$ $(1\leq i\leq4)$ are Brownian motions defined on the complete probability space $(\Omega, \mathscr{F}, \mathbb{P})$ adopting to the filtration $\{ \mathscr{F}_t\}_{t\geq 0}$. We also let $X(t) = (T(t), I(t), V(t), E(t), Z_1(t)$, $Z_2(t), Z_3(t))$, $X_0 = (T(0), I(0), V(0), E(0), Z_1(0)$, $Z_2(0), Z_3(0))$, and $\mathbb{R}_+^7 = \{X = (X_1, X_2$, $X_3, X_4, X_5, X_6, X_7)\in\mathbb{R}^7: X_j > 0, 1\leq j\leq7\}$. Thus, we use $a\wedge b$ to denote $\min\{a, b\}$, and use $a\vee b$ to represent $\max\{a, b\}$.

For convenience, we introduce the following symbols

The following result shows that system (2.5) has a unique positive global solution.

Theorem 2.1. System (2.5) has a unique and positive solution $X(t)$ with the initial value $X_0\in \mathbb{R}_{+}^7$ for all $t\geq 0$, and the solution will remain in $\mathbb{R}_+^7$ with probability one, namely, $X(t)\in\mathbb{R}_+^7$ for all $t\geq 0$ almost surely (a.s.).

Proof. Following the theory of stochastic differential equation in Mao's book ^[44], it is clear that the coefficients of system (2.5) are locally Lipschitz continuous. Therefore, stochastic system (2.5) exists a unique local solution $X(t)$ on $t \in [0, \rho_e)$, where $\rho_e$ is the explosion time.

Next, we demonstrate the solution is global, that is, we need to show $\rho_e = \infty$ a.s.. By using reduction to absurdity, suppose that there exists a finite time, such that every component of solution $X(t)$ could not explode to infinity. Let $m_0 > 0$ be large enough such that for every component of $X_0$ located in the interval $[\frac{1}{m_0}, m_0]$. For each integer $m\geq m_0$, define the stopping time

where we set $\inf \emptyset = \infty$ ($\emptyset$ is the empty set). Obviously, $\rho_m$ is increasing as $m\rightarrow \infty$. Denote $\rho_\infty = \underset{m\rightarrow \infty}{\lim} \rho_m$, then $\rho_\infty \leq \rho_e$ a.s.. To illustrate $\rho_\infty = \infty$, we validate it in two cases: (ⅰ) If $\rho_\infty = \infty$ is true a.s., then $\rho_e = \infty$ a.s., so $X(t)\in \mathbb{R}_+^7$ for all $t\geq 0$ a.s.. (ⅱ) If $\rho_\infty < \infty$, assuming there exists a pair of constants $\widetilde{t} > 0$ and $\epsilon\in(0, 1)$ such that $\mathbb{P}\{ \rho_\infty \leq \widetilde{t}\} > \epsilon$. Then, there exists an integer $m_1\geq m_0$ such that

is established.

Construct a $C^2$-function $W$: $\mathbb{R}_+^7 \rightarrow \mathbb{R}_+$,

where, $a$ and $b$ are positive constants which will be determined later, and $\eta_i$ $(1\leq i\leq5)$ are defined in Eq (2.6). Applying It$\hat{\text{o}}$'s formula ^[44] to $W$, we have

where $LW$: $\mathbb{R}^7_+ \rightarrow \mathbb{R}_+$ is defined by

Here, $K_1 = \lambda+a\mu_1+b\mu_2$ $+\mu_3+\mu_4 +r_1+s_1+r_2+s_2+r_3+s_3+\frac{1}{2}a\sigma_1^2+\frac{1}{2}b\sigma_2^2+\frac{1}{2}\sigma_3^2+\frac{1}{2}\sigma_4^2$, and $K_1$ is a positive constant. We further choose the constants $0 < a\leq\frac{1}{4}\left(\frac{\mu_2}{\beta_2} \wedge \frac{\eta_1\mu_3}{\beta_1}\right)$, $0 < b\leq\frac{2\eta_2\mu_4}{q}$, such that $a\beta_2-\frac{\mu_2}{4}\leq0$, $a\beta_1-\eta_1\mu_3\leq0$ and $bq -2\eta_2\mu_4\leq0$ simultaneously. Thus

We therefore obtain

Integrating both sides from 0 to $\rho_m \wedge \widetilde{t}$, we obtain

Taking expectations on both sides, we further get

Set $\Omega_m = \{\rho_m\leq \widetilde{t}\}$ for $m\geq m_1$. By Eq (2.7), we derive $\mathbb{P}(\Omega_m)\geq \varepsilon$. Note that for every $\omega\in\Omega_m$, there is some $j$ $(1\leq j\leq 7)$ such that $X_j(\rho_m, \omega)$ equals either $m$ or $\frac{1}{m}$, and hence $W(X(\rho_m, \omega))$ is no less than either

or

where $\eta = \eta_1\wedge(2\eta_2)\wedge\eta_3\wedge\eta_4\wedge\eta_5$. Thereafter, we have

From Eq (2.8), we derive

where ${\bf{I}}_{\Omega_m}$ is the indicator function of $\Omega_m$. Letting $m\rightarrow\infty$, then

This leads to the contradictions, so we must have $\rho_\infty = \infty$ a.s.. This completes the proof.

3.   Stationary Markov process

In this section, we mainly focus on the persistence of each population. For deterministic model, we need to show the global asymptotic stability of the positive equilibrium. For stochastic differential equation system with distributed delays, we need to prove the existence of stationary Markov process.

Firstly, we introduce some knowledge about stochastic differential equations. Consider the $d$-dimensional time-homogeneous stochastic differential equation of It$\hat{\text{o}}$ type

with initial value $X(t_0) = X_0\in\mathbb{R}^d$. By the definition of stochastic differential, this equation is equivalent to the following stochastic integral equation,

Summarizing Theorems 3.4 and 3.7 in Khasminskii's book ^[45], we derive the following lemma.

Lemma 3.1. It is assumed that the vectors $b(X)$, $\sigma_1(X)$, $\ldots$, $\sigma_d(X)$ ($t\geq t_0$, $X\in\mathbb{R}^d$) are continuous functions of $X$.

(A1) There is a constant $B$ with the properties

(A2) There exists a non-negative $C^2$-function $U(X)$ in $\mathbb{R}^d$ such that $LU(X)\leq -1$ outside some compact set.

If these two conditions are satisfied, then there exists a solution of system (2.5) which is a stationary Markov process.

Define the critical condition $\mathcal{R}_0^s$ as follows

Denote

Theorem 3.1. If $\mathcal{R}_0^s > 1$, then the solution $X(t)$ of system (2.5) is a stationary Markov process.

Proof. We employ Lemma 3.1 to prove the existence of stationary Markov process of system (2.5). In Lemma 3.1, condition (A1) is to ensure that the solution of the system is a Markov process, and condition (A2) is to show that the solution of the system is stationary. Obviously, the coefficients of system (2.5) are continuous functions of $X$ and satisfy condition (A1) in Lemma 3.1, which means that the solution of the system is a Markov process. In the following, we illustrate condition (A2) in Lemma 3.1 by constructing appropriate nonnegative functions.

Denote

$f$ is a quadratic function with respect to $T$. Suppose $T_0$ and $-T_1$ are the roots of $f(T) = 0$, then

As $-(T-T_0)^2 = -(T-T_0)(T+T_1)+(T-T_0)(T_0+T_1)\leq 0$, then $(T-T_0)(T+T_1) \geq (T-T_0)(T_0+T_1)$, and

Define a function

then,

Applying It$\hat{\text{o}}$'s formula ^[44] to $W_1$, we obtain

Here, $\alpha_1 = \beta_1+\frac{\beta_2\mu_3}{k}$, and we use the equality $\mu_1-r+\frac{rT_0}{T_{max}} = \frac{\lambda}{T_0}$ at the infection-free equilibrium $E_0$ of system (2.1). Denote

and we have

where $\alpha_2 = \frac{pqr_3\mu_3}{k\mu_4(r_3+s_3)}$.

Construct a $C^2$-function $U_1$: $\mathbb{R}_+^7\rightarrow \mathbb{R}$,

where $a_i$ $(i = 1, 2, 3)$ and $b_j$ $(j = 1, 2)$ are positive constants which will be determined later. By It$\hat{\text{o}}$'s formula, we have

Let

and

We calculate that,

Denote $\alpha_3 = a_3\alpha_1+b_1\alpha_1+\alpha_2$. Consequently,

where, $\mathcal{R}_0^s$ is defined in Eq (3.1), and

Define a $C^2$-function $U$: $\mathbb{R}_+^7\rightarrow \mathbb{R}$, in the following form,

where,

$0 < \theta < \min\left\{ 1, \frac{1}{4\sigma^2}(\mu_2\wedge 4\mu_3\wedge4\mu_4) \right\}$ (see Eq (3.2) for the definition of $\sigma^2$), and $\eta_i$ $(1\leq i\leq 5)$ are defined in Eq (2.6). We choose a suitable constant $M > 0$ to satisfy the following condition

where,

and $K_2$ is a positive constant. It is easy to obtain that

where $D_l = (\frac{1}{l}, l)\times (\frac{1}{l}, l)\times(\frac{1}{l}, l)\times(\frac{1}{l}, l)\times(\frac{1}{l}, l)\times(\frac{1}{l}, l)\times(\frac{1}{l}, l)$ and $l$ is a positive integer. Since $U(X)$ is a continuous function, $U(X)$ must have a minimum point $\overline{X}_0$ in the interior of $\mathbb{R}_+^7$. Hence, we define a nonnegative $C^2$-function $\overline{U}: \mathbb{R}_+^7\rightarrow \mathbb{R}$ as follows

Using It$\hat{\text{o}}$ formula, we obtain

For the convenience of calculation, we simplify the following

where $\eta_6 = \frac{r_3+s_3}{2}$. Hence,

See Eq (3.2) for the definition of $\sigma^2$, and

Consequently, we summarize the above calculations and obtain

See Eq (3.3) for the expression of $K_2$.

Next, we construct a compact subset $D_\varepsilon$ to make $L\overline{U} < 1$ valid. Define a bounded closed set as below

where $\varepsilon$ is a sufficiently small positive constant. In set $\mathbb{R}_+^7\backslash D_\varepsilon$, this sufficiently small positive constant $\varepsilon$ satisfies the following conditions

where,

Then, we separate $\mathbb{R}_+^7\backslash D_\varepsilon$ to fourteen domains,

Clearly, $D_\varepsilon^c = \bigcup\limits_{i = 1}^{14}D_{j}$.

Case 1. When $X\in D_1$,

According to (3.4), it implies that $L\overline{U} \leq -1$ for any $X\in D_1$.

Case 2. When $X\in D_2$,

In view of (3.5), we have $L\overline{U} \leq -1$ for any $X\in D_2$.

Case 3. When $X\in D_3$,

According to (3.6), we deduce that $L\overline{U} \leq -1$ for any $X\in D_3$.

Case 4. When $X\in D_4$,

According to (3.7), it implies that $L\overline{U} \leq -1$ for any $X\in D_4$.

Case 5. When $X\in D_5$,

Based on (3.8), we derive that $L\overline{U} \leq -1$ for any $X\in D_5$.

Case 6. When $X\in D_6$,

For any $X\in D_6$, we obtain that $L\overline{U} \leq -1$ under the condition (3.9).

Case 7. When $X\in D_7$,

By condition (3.10), we conclude that $L\overline{U} \leq -1$ for any $X\in D_7$.

Case 8. When $X\in D_8$,

It follows that $L\overline{U} \leq -1$ for any $X\in D_8$ if condition (3.11) is satisfied.

Case 9. When $X\in D_9$,

By condition (3.12), we derive that $L\overline{U} \leq -1$ for all $X\in D_9$.

Case 10. When $X\in D_{10}$,

From condition (3.13), we get that $L\overline{U} \leq -1$ for any $X\in D_{10}$.

Case 11. When $X\in D_{11}$,

It follows that $L\overline{U} \leq -1$ for any $X\in D_{11}$ if the condition (3.14) is satisfied.

Case 12. When $X\in D_{12}$,

In view of (3.15), we have $L\overline{U} \leq -1$ for any $X\in D_{12}$.

Case 13. When $X\in D_{13}$,

It leads to $L\overline{U} \leq -1$ for any $X\in D_{13}$ if the condition (3.16) is satisfied.

Case 14. When $X\in D_{14}$,

Under the condition (3.17), we conclude that $L\overline{U} \leq -1$ is satisfied for any $X\in D_{14}$.

Consequently, under the conditions (3.4)–(3.17), there exists a sufficiently small $\varepsilon$, such that

According to Lemma 3.1, we obtain that the solution of system (2.5) is a stationary Markov process. This completes the proof.

By the theory of Khasminskii ^[45], we derive that system (2.5) has a stationary Markov process when the critical condition $\mathcal{R}_0^s$ is greater than one. We should mention that

This means that when there is no white noises, the critical condition $\mathcal{R}_0^s$ of stochastic differential equation (2.5) is reduced to the basic reproduction number $\mathcal{R}_0$ of its corresponding deterministic differential equation (2.1). The result shows that the existence of stationary Markov process in our stochastic model is the extension of its corresponding deterministic model to the stability of the positive equilibrium.

4.   Extinction

In the course of viral infection, we are also concerned about the extinction of the virus. In this section, we derive the sufficient conditions to ensure the extinction of HIV virus theoretically.

Denote

where,

$Q$ is a constant such that $\int_0^\infty \mu(x)dx = 1$, and see the expression of $\eta_1$ in Eq (2.6) of Section 1.

Theorem 4.1. Suppose $X(t)$ be the solution of system (2.5) with the initial value $X_0\in\mathbb{R}_+^7$, then the solution $X(t)$ of system (2.5) has the following property

In particular, if $\widehat{\mathcal{R}}_0 < 1$ holds, then

and

It indicates that the virus can be eradicated with probability one a.s..

Proof. From the first equation of system (2.5), we obtain that

Consider the following auxiliary equation with stochastic differential equation

Let $x(t)$ be the solution of system (4.1) with the initial value $x(0) = T(0) > 0$. By Theorem 3.1 in literature ^[36], we obtain that system (4.1) has the ergodic property with ergodic distribution

where $Q$ is a constant such that $\int_0^\infty \mu(x)dx = 1$. Then, we have

By the comparison theorem of stochastic differential equation ^[46], we further obtain that

Define

and see Eq (2.6) for the expressions of $\eta_1$, $\eta_2$, $\eta_5$, $\eta_7$ and $\eta_8$. Applying It$\hat{\text{o}}$'s formula, we obtain

Notice that

Then, we have

Applying the inequality $(a+b+c)^2\leq 3(a^2+b^2+c^2)\; (a, b, c > 0)$, we get

We further have

For inequality (4), integrating both sides from $0$ to $t$, and dividing by $t$ on both sides, we obtain

Taking the superior limit on both sides of inequality (4.3) and combining with inequality (4.2), under the critical condition $\widehat{\mathcal{R}}_0 < 1$, we outline that

which means that

This completes the proof.

5.   Numerical simulations

We have theoretically analyzed the existence of stationary Markov process and the extinction for virus in Sections 3 and 4. In this section, in order to study the viral dynamics of a delayed HIV stochastic model with cell-to-cell infection and CTL immune response, we carry out numerical simulations on two aspects: (ⅰ) the influence of random fluctuations on the virions and the CTLs populations; (ⅱ) the effect of cell-to-cell infection and time delays on the number of target cells, infected T-cells, virions and CTLs.

In the following, we give numerical simulations to show the effect of the random fluctuations and the delays on the long time behavior around the positive equilibrium $E^*$. By employing the Milstein's higher order method in Higham ^[47], the discretization form of model (2.5) is

where the time increment $\Delta t = 0.01$ in our simulations, and $\xi_{1, m}$, $\xi_{2, m}$, $\xi_{3, m}$ and $\xi_{4, m}$, $m = 1, 2, \cdots, n$, are the $m$th realization of the four independent Gaussian random variables with distribution $N(0, 1)$.

For the weak kernels $f_i(\tau) = r_i e^{-r_i\tau}$ $(i = 1, 2, 3)$, we choose $r_1 = r_2 = r_3 = 10$, $s_1 = s_2 = 0.2$, $s_3 = 0.5$. For the deterministic model (2.1), all the other parameter values are from Table 1. By Matlab software, we compute that

and the unique positive equilibrium $E^* = (253.2461, 3.5997,156.5071, 6.8565)$. Following the theoretical results, we know that the positive equilibrium $E^*$ is globally attractive.

Example 5.1 For stochastic model (2.5), in order to examine the existence of stationary Markov process and the effect of random fluctuations on viral dynamics numerically, we choose three groups of random noise $(\sigma_1, \sigma_2, \sigma_3, \sigma_4)$ equal to $(0.02, 0.04, 0.4, 0.02)$, $(0.04, 0.08, 0.8, 0.04)$ and $(0.06, 0.12, 1.8, 0.06)$, respectively. The remaining parameter values of system (2.5) are shown in Table 1, then the critical values of $\mathcal{R}_0^s$ corresponding to the three groups of noise are $2.0432$, $1.3381$ and $0.9653$, respectively. Theorem 3.1 is satisfied for the first two groups of random noise.

The numerical simulations show that the stationary Markov process occurs (see Figures 1(a) and 2(a)) and the corresponding histograms of the solution for virus population and CTLs population can be seen in Figures 1(b) and 2(b), respectively. It is observed that, with the increase of noise intensity, the amplitude of virus and CTLs populations becomes large, and small noise intensity may contribute to maintain the existence of stationary Markov process even though the critical condition $\mathcal{R}_0^s$ is less than one.

To further study the effect of random noises on viral dynamics, we assume that there is only one random noise, and observe the effect of this noise on the number of virions and CTLs. From Figures 3 and 4, we find that the the smaller the noise intensity is, the smaller the fluctuation amplitude of virus and CTLs populations number is. With the increase of the noise intensity, the fluctuation amplitude of population increases. This indicates that the noise intensity can affect the fluctuation range of the population.

Example 5.2 Consider model (2.5) with noise intensity $(\sigma_1, \sigma_2, \sigma_3, \sigma_4)$ = $(0.5, 1.0, 6.0, 0.5)$, and all the other parameter values are the same as in Example 5.1. The critical values of $\mathcal{R}_0^s = 0.0599 < 1$. Figure 5 shows that the big noise intensity can make the infected T-cells, virus and CTLs population extinct, while its corresponding deterministic model (2.1) has a attractive positive equilibrium.

Example 5.3 To study the effect of the cell-to-cell infection on model behavior, we compare our stochastic model (2.5) to the stochastic model without cell-to-cell infection. We choose $(\sigma_1, \sigma_2, \sigma_3, \sigma_4)$ = $(0.04, 0.08, 0.8, 0.04)$, and all the other parameter have the same values as in Table 1. Following the definition of critical condition $\mathcal{R}_0^s$ in stochastic model (2.5), we calculate that

Thus, the critical condition of the stochastic model without cell-to-cell infection is $0.6787$, which is less than one. In Figure 6, we can see that the model without cell-to-cell infection could underestimate the number of infected T-cells, virions and CTLs, and overestimate the number of healthy T-cells. Thus, with same noise intensity, the amplitude of each population in the stochastic model without cell-to-cell infection is smaller than that in the stochastic model with cell-to-cell infection.

To study the effect of the delays on model behavior, we compare our stochastic model (2.5) to the stochastic model without delays. We take $r_1 = r_2 = r_3 = 8$, $s_1 = s_2 = s_3 = 2$ for the weak kernels (2.2), the noise intensity $(\sigma_1, \sigma_2, \sigma_3, \sigma_4)$ = $(0.04, 0.08, 0.8, 0.04)$, and all the other parameter values are from Table 1. By computing, for stochastic model (2.5) with distributed delays, we have

for the stochastic model without distributed delays, we have

In Figure 7, we examine that the delays have no significant effect on the number of target T-cells, infected T-cells and virus populations, except for the number of CTLs population. Thus, the stochastic model without delays has no evident impact on the oscillation amplitude of each population.

Examples 5.1 and 5.2 reveal that the small noise intensity can keep the number of virions and CTLs under a certain range, while the big noise intensity can lead to the extinction of the virus even though its corresponding deterministic model has a attractive positive equilibrium. Examples 5.1 and 5.3 indicate that the fluctuation amplitude of population is more sensitive to the noise intensity than the delay, since the fluctuation amplitude of each population changes within a very narrow range with respect to the delay (see Figure 7), and it changes within a wide range with respect to the noise intensity (see Figures 3 and 4). Example 5.3 also demonstrate the stochastic model without cell-to-cell infection could underestimate the number of virions and CTLs.

6.   Conclusions

In this paper, white noises are used to describe the random fluctuations during HIV infection process. We have formulated a stochastic HIV model which includes virus-to-cell infection, cell-to-cell infection, CTL immune response and three distributed delays. For the commonly used gamma distribution delays, we choose the weak kernels form as our study. To my knowledge, few articles have studied the cell-to-cell infection and delays on stochastic HIV model. By transforming the four-dimensional stochastic integro-differential equation into a degenerate seven-dimensional stochastic differential equation, we theoretically obtain three main results: (Ⅰ) The solution of the system is unique and global. (Ⅱ) By constructing suitable Lyapunov functions, we derive the existence of stationary Markov process when the critical condition is greater than one, which implies the persistence of the virus. (Ⅲ) Sufficient conditions are given to ensure the extinction of the virus.

According to the actual parameters obtained in previous references, three main results of system (2.5) are obtained numerically: (Ⅰ) Within the scope of small noise intensity, the smaller the noise is, the smaller the amplitude of the system solution vibration is. Small noise intensity is helpful to keep the number of virions and CTLs fluctuating within some certain range. (Ⅱ) For stochastic model, sufficiently large noise intensity may induce the extinction of virus population even if its corresponding deterministic model has a stable positive equilibrium. (Ⅲ) Cell-to-cell infection can affect the number of each population, while the delay has no significant effect on the number of each population. It indicates that random white noise is more sensitive to the dynamics on the model than the delay.

Compared with HIV stochastic model without distributed delay ^{[33,34,35,36]}, stochastic model with distributed delay can be transformed into a degenerate stochastic differential equation. As far as we know, little work has been done on the theoretical analysis of the degenerate differential equations. Comparing our stochastic model with the model including only one infection mode (virus-to-cell infection or cell-to-cell infection) ^{[33,34,35,36,37,38]}, we find that under the same noise intensity, the model including only one infection mode could underestimate the number of virions and CTLs. Thus, our study can be regard as an extension of the earlier works ^{[33,34,35,36,37]}.

Acknowledgments

The authors would like to thank the referees for their valuable suggestions. This work is supported by National Natural Science Foundation of China (No. 11401589, No. 11501446, No. 11871473, No. 11801566), the Fundamental Research Funds for the Central Universities (No. 17CX02066, No. 18CX02049A), Shandong Provincial Natural Science Foundation (No. ZR2019MA010), the Natural Science Research Fund of Northwest University (14NW17), and the Scientific Research Plan Projects of Education Department of Shaanxi Provincial Government (15JK1765).

Conflict of interest

The authors declare that they have no conflict of interest.