Stationary distribution, extinction and density function for a stochastic HIV model with a Hill-type infection rate and distributed delay

1.   Introduction

Human immunodeficiency virus (HIV) is the one of the most harmful diseases, and thus has always been a hot topic worthy of attention. Considering the crowding effect of the $CD4^+ T$ cells (prey of virus), Bairagi and Adak^[1] proposed the HIV-1 system with a Hill type function based on the following mass action principle:

Here $x(t), y(t)$ and $v(t)$ denote the concentrations of activated $CD4^+ T$ cells, infected $CD4^+ T$ cells and virus particles, respectively. $n\geq 1$ is a Hill constant. Other parameters can be referred to in ^[1]. Besides, from ^[1], the basic reproductive number for System (1.1) is as follows:

And when $R_1 < 1$, the disease-free equilibrium $E_0(\frac{s}{\mu}, 0, 0)$ is globally asymptotically stable, and when $R_1 > 1$, $E_0$ is unstable and there exists a unique endemic equilibrium, which is globally asymptotically stable. Considering that the growth of infected $CD4^+ T$ cells is affected not only by the number of normal cells, but also by the number of previous cells, we investigate the HIV system with distributed delay:

where the kernel $f(t) = \frac{t^k\sigma^{k+1}e^{-\sigma t}}{k!}$ is a gamma distribution initially proposed by Macdonald^[2] and $\sigma > 0$ is a mean delay. For convenience, we choose the weak kernel case ($k = 0$), that is, $f(t) = \sigma e^{-\sigma t}$. Let

Then System (1.3) is transformed into the equivalent form:

Xu ^[3] investigated the global asymptotic stability of an HIV-1 infection model with distributed intracellular delays. Using methods similar to those in Theorems 2.3 and 2.4 of ^[1] and Theorems 3.1 and 3.2 of ^[3], we obtain the following results.

Theorem 1.1. (i) The disease-free equilibrium point $E_0(\frac{s}{\mu}, 0, 0, 0)$ of System (1.4) is globally asymptotically stable when $R_1 < 1$ and unstable when $R_1 > 1$, where $R_1$ is defined in Eq.(1.2);

(ii) The endemic equilibrium point $E^*(x^*, y^*, v^*, w^*)$ of System (1.4) is globally asymptotically stable when $R_1 > 1$, where

with $A = \frac{(\mu+\alpha)\gamma}{c\alpha}$.

In addition, the infectious diseases are inevitably affected by environmental white noises^[4,5,6]. In view of the above, we consider the following stochastic HIV-1 model with distributed delay:

where $B_j(t), \; j = 1, 2, 3$ represent independent Brownion motions whose noise intensities are expressed as $\sigma_j^2, \; j = 1, 2, 3.$ In ^[7,8,9], the existence of a stationary distribution and the extinction of stochastic systems are studied based on the theory of Khasminskii ^[10] by constructing the suitable Lyapunov function, which implies that the diseases will be prevalent or tend to extinction. In ^[11], Guo and Zhang gave the numerical approximation for an HIV infection model incorporating the mean-reverting Ornstein-Uhlenbeck process. In ^[12], the extinction and the existence of a unique invariant probability measure for a stochastic HIV-1 infection model with degenerate diffusion were obtained. In ^[13], a group of stochastic dynamic models of the HIV/AIDS infection in a host population are presented, and global asymptotic and $p-$exponential stability of the disease-free equilibrium in probability was investigated. The most difficulty with our work is determining how to deal with the Hill-type infection rate and distributed delay when constructing the Lyapunov function and proving the positive definiteness.

Similarly, let

Then

Thus, System (1.5) is transformed into the following equivalent form:

In comparison with other existing results, the achieved contributions and innovations can be summarized as follows:

● A stochastic HIV model with a Hill-type infection rate and distributed delay is proposed, which may reflect more reality than the existing ones.

● The existence of a stationary distribution for System (1.6) is obtained by constructing a suitable Lyapunov function, which determines a critical value $R_0^*$ corresponding to the basic reproduction number.

● The exact density function near the endemic quasi-equilibrium is given by solving the Fokker-Planck equation.

● Our main innovation is the development of a technique to deal with the Hill-type infection rate, which is different with the existing ones.

The rest of the article is organized as follows. In Section 2, the sufficient condition of the existence of a stationary distribution for the stochastic system given by Eq.(1.5) is derived; it determines a sharp critical value $R_0^s$. In Section 3, the extinction of the diseases is investigated. In Section 4, the exact probability density function at the quasi-endemic equilibrium is derived. In Section 5, the numerical results are illustrated. Finally, the conclusion is given briefly.

2.   Stationary distribution for the stochastic system given by Eq.(1.6)

The existence and uniqueness of the global positive solution of System (1.6) will be given. Since this is standard, we omit it.

Theorem 2. For any initial value $(x (0), y (0), v(0), w(0))\in R_+^4$, there exists a unique positive solution $(x(t), y(t), v(t), w(t))$ for System (1.6) on $t \ge 0$ and the solution will remain in $R_+^4$ with a probability of one.

Consider the following auxilary Logistic equation:

Similar to Lemma 4.1 in ^[14] and ^[15], we have the following lemma:

Lemma 1. ^[14,15] Eq (2.1) has a unique stationary distribution with the density function $f^*(\cdot)$ defined by

where $a_1 = \frac{2\mu+\sigma_1^2}{\sigma_1^2}, \; b_1 = \frac{2s}{\sigma_1^2}$ and $\varGamma(\cdot)$ is a Gamma function and the following equalities hold:

Proof. Eq.(2.1) has a unique stationary distribution with the density function $f^*(z)$ on $(0, \infty)$, which is defined by Eq.(2.2) and Eq.(2.3) holds. In addition,

Next, we consider the following integral equation:

Lemma 2. ^[10,16] Suppose that System (2.4) has a global positive solution and there exists a non-negative function $V\in C^2(R_+^4, R_+)$ and a bounded closed set $D$ such that $LV\leq-1$ for $R_+^4\setminus D$. Then System (2.4) has a stationary distribution.

Next by applying a method similar to those in ^[7,8,9], which are based on the theory of Khasminskii ^[10] and combining Theorem 2.1 and Lemmas 2.1 and 2.2, we obtain the following main result.

Theorem 3. Assuming that $R_0^s > 1$ and $\mu > \frac{(n-1)\sigma_1^2}{2}$, System (1.6) has a unique stationary distribution $\pi(\cdot)$, where

where $a_1 = \frac{2\mu+\sigma_1^2}{\sigma_1^2}$ and $b_1 = \frac{2s}{\sigma_1^2}$.

Proof. By Theorem 2 and Lemma 2, we only need to construct a Lyapunov function. For convenience, denote $F(x) = \frac{x^n}{a^n+x^n}$. Noting that $\frac{d}{dx}\left(\frac{F(x)}{x}\right) = \frac{x^{n-2}((n-1)a^n-x^n)}{(a^n+x^n)^2},$ we have that $\frac{F(x)}{x}\leq\begin{cases} 1/a, & n = 1\\ (F((n-1)^{\frac{1}{n}}a))/((n-1)^{\frac{1}{n}}a), & n > 1 \end{cases}\triangleq C_1 \; \text{for}\; x\in (0, +\infty).$ By the first equation of Eq.(1.6), we have $dx\leq (s-\mu x)dt+\sigma_1 xdB_1(t)$. By the comparison principle, we have that $x(t)\leq X(t), \; \text{a.s.}$ Define

where

By Ito's formula, Lemma 1 and $\frac{F(x)}{x}\leq C_1$, we have that,

where $R_0^s$ is defined in Eq.(2.5), and

The It$\hat{o}$ formula is applied to

where

Then

and

where $\rho = k-\frac{\theta p}{2} > 0$, and

$A: = \max\limits_{(x, y, v, w)\in R_+^4}\{s\left(x+\frac{y}{2}+\frac{(\alpha+\mu)v}{4c\alpha}+\frac{\beta w}{\sigma}\right)^{\theta}-\frac{\rho}{2}\left(x+\frac{y}{2}+\frac{(\alpha+\mu)v}{4c\alpha}+\frac{\beta w}{\sigma}\right)^{\theta+1}\}$.

Define

where $M$ is a positive constant sufficiently enough satisfying

By Eqs.(2.6)–(2.10), we have that

where

First, we consider the expression $G(x, y, v, w)$ in two cases:

Case 1. If $y\rightarrow 0^+$, then $G(x, y, v, w)\leq -2+(M+1)\frac{\beta C_1 c \alpha y}{\gamma}\rightarrow -2.$

Case 2. If $x\rightarrow \infty$ or $x\rightarrow 0^+$ or $y\rightarrow \infty$ or $v\rightarrow \infty$ or $w\rightarrow \infty$ or $y\nrightarrow 0^+, \; v\rightarrow 0^+$ or $x\nrightarrow 0^+, \; v\nrightarrow 0^+, \; \; w\rightarrow 0^+$, then

where $H = \sup\limits_{y\in (0, +\infty)}\{(M+1)\frac{\beta C_1 c \alpha y}{\gamma}-\frac{\rho}{4}\left(\frac{y}{2}\right)^{\theta+1}\}$.

In light of the above, there exists a sufficiently small constant $\varepsilon > 0$ such that

where $D_{\varepsilon} = \{(x, y, v, w)\in R_+^4|\varepsilon\leq x\leq\frac{1}{\varepsilon}, \; \varepsilon\leq y\leq\frac{1}{\varepsilon}, \; \varepsilon ^2\leq v\leq\frac{1}{\varepsilon ^2}, \; \varepsilon ^3\leq w\leq\frac{1}{\varepsilon ^3} \}.$

By the continuity of $G(x, y, v, w)$, there exists a positive constant $K$ such that

Hence, integrating from $0$ to $t$ and taking the expectation for both sides of Eq.(2.11) give that

Let $t\rightarrow \infty$, then,

which implies that,

where $P(t, (x, y, v, w), \cdot)$ is the transition probability of $(x(t), y(t), v(t), w(t))$. By the invariance of $R_+^4$ and the inequality given by Eq.(2.12), there exists an invariant probability measure $\pi(\cdot)$ on $R_+^4$, (see ^[17,18]). The result is confirmed.

3.   Extinction of the diseases of System (1.6)

Theorem 4. Let $(x(t), y(t), v(t), w(t))$ be a solution of System (1.6) with any initial value $(x(0), y(0), v(0), w(0))\in R_+^4$. Then the following results will hold:

where

where $R_1$ is defined by Eq.(1.2).

Especially, if $\bar{\mu} < 0$, then the infected $CD4^+ T$ cells population $y(t)$ and virus particle $v(t)$ will die out exponentially with a probability of one.

Proof. Define

Applying the $It\hat{o}$ formula,

where

and the monotone increasing property of $F(X)$ with $M_s = \left(\begin{array}{ccc} 0 & 0 & \frac{\beta}{\alpha+\mu} \\ \frac{c\alpha}{\gamma} & 0 & 0 \\ 0 & F(\frac{s}{\mu}) & 0 \\ \end{array} \right),$ and $\lambda = \sqrt[3]{R_1} = \sqrt[3]{\frac{\beta c \alpha F(\frac{s}{\mu})}{(\alpha+\mu)\gamma}}$ satisfies $(\delta_1 (\alpha+\mu), \delta_2\gamma, \delta_3\sigma)M_s = \lambda(\delta_1(\alpha+\mu), \delta_2\gamma, \delta_3\sigma).$ Integrating both sides of Eq.(3.2) yields that,

where $M_1(t) = \int_{0}^{t} \frac{\delta_1\sigma_2y(\tau)}{z(\tau)}\, dB_2(\tau)\; \; \text{and} \; \; M_2(t) = \int_{0}^{t} \frac{\delta_2\sigma_3v(\tau)}{z(\tau)}\, dB_3(\tau).$ Noting that the quadratic variation $\frac{\langle M_1, M_1\rangle_t}{t} = \frac{\sigma_2^2}{t}\int_0^t\frac{\delta_1^2y^2(\tau)}{z^2(\tau)}d\tau\leq \sigma_2^2 < \infty, \; \; \frac{\langle M_2, M_2\rangle_t}{t}\leq\sigma_3^2 < \infty,$

From the ergodic theorem and Eq.(3.3) and Eq.(3.4), it follows that,

If $\bar{\mu} < 0$, it is obvious to get $\lim\limits_{t \to +\infty}z(t) = 0$, a.s., which implies that $\lim\limits_{t \to +\infty}y(t) = \lim\limits_{t \to +\infty}v(t) = \lim\limits_{t \to +\infty}w(t) = 0$, a.s. In the other words, infected cells will exponentially decrease to zero with a probability of one.

4.   Probability density function of System (1.6)

We have obtained the existence of the stationary distribution of System (1.6). Next we give the details about the local probability density function of System (1.6).

First, let $\bar{x} = \ln{x}$, $\bar{y} = \ln{y}$, $\bar{v} = \ln{v}$, $\overline{w} = \ln{w}$, then by It$\hat{o}$'s formula, System (1.6) is transformed into the following form:

where $\mu_i = \mu + \frac{\sigma_i^{2}}{2}, \; i = 1, 2 \; \text{and}\; \mu_3 = \gamma + \frac{\sigma_3^{2}}{2}$.

Define

which is the same as Eq.(1.2) when $\sigma_i = 0(i = 1, 2, 3)$. Define a quasi-endemic equilibrium $\widetilde{E^*} = (\widetilde{x^*}, \widetilde{y^*}, \widetilde{v^*}, \widetilde{w^*}) = (e^{{\overline{x^*}}}, e^{{\overline{y^*}}}, e^{{\overline{v^*}}}, e^{{\overline{w^*}}})$, where

where $\widetilde{A}$ = $\frac{(\mu_2 + \alpha)\mu_3}{c\alpha}$. And $\widetilde{E^*}$ exists if and only if $R_1^s > 1$, and $\widetilde{E^*}$ is the same with the endemic equilibrium $E^*$ of the determined system given by Eq.(1.4) when there is no white noises.

Let $(\hat{x}, \hat{y}, \hat{v}, \hat{w})$ = $(\bar{x}-\overline{x^{*}}, \bar{y}-\overline{y^{*}}, \bar{v}-\overline{v^{*}}, \bar{w}-\overline{w^{*}})$, then the linearized system of Eq.(4.1) at $(\overline{x^*}, \overline{y^*}, \overline{v^*}, \overline{w^*})$ is as follows:

where, combining (4.3),

It is obvious that

Theorem 5. Let $Y = (\hat{x}, \hat{y}, \hat{v}, \hat{w})$ be a solution to Eq.(4.4) with the initial value $(\hat{x}(0), \hat{y}(0), \hat{v}(0), \hat{w}(0)) \in R^4.$ If $R^s_1 > 1,$ then there exists a unique density function $\Phi (Y)$ around the quasi-equilibrium $\widetilde{E^*}$, which can be approximately expressed in the following form

in which $\Sigma = \Sigma_1 + \Sigma_2 + \Sigma_3$ is positive definite, $R_1^s$ is defined by Eq.(4.2), $\Sigma_1 = \rho_1^2 M_1^{-1} \Theta_1 (M_1^{-1})^T$ and $\Sigma_2, \; \Sigma_3$ will be described below for different cases.

Proof. System (4.4) can be rewritten into $dY = AYdt + \Lambda dB(t)$, where $Y = (\hat{x}, \hat{y}, \hat{v}, \hat{w})^T$, $\Lambda = diag(\sigma_1, \sigma_2, \sigma_3, 0)$, $B(t) = (B_1(t), B_2(t), B_3(t), 0)$ and

According to ^[19], the four dimensional Fokker-Planck equation to describe a density function $\Phi (Y) = \Phi(\hat{x}, \hat{y}, \hat{v}, \hat{w})$ of the stationary distribution of Eq.(4.4) around the quasi-equilibrium $\widetilde{E^*}$ is as follows:

which is an approximate representation of the Gaussian distribution $\Phi(Y) = k e^{-\frac{1}{2}(Y-Y^*)Q(Y-Y^*)^T},$ with $Y^* = (0, 0, 0, 0)$; also, Q is a real symmetric matrix, which satisfies $Q\Lambda^2 Q + A^TQ + QA = 0$. If Q is positive-definite, let $Q^{-1} = \Sigma$, then $\Lambda^2 + A\Sigma + \Sigma A^T = 0$. By the finite independent superposition principle,

where $\Sigma = \Sigma_1 + \Sigma_2 + \Sigma_3, \; \Lambda^2 = \Lambda_1^2 + \Lambda_2^2 + \Lambda_3^2, \; \Lambda_1 = diag(\sigma_1, 0, 0, 0), \; \Lambda_2 = diag(0, \sigma_2, 0, 0)$, $\Lambda_3 = diag(0, 0, \sigma_3, 0).$

Step 1. For System (4.4), we consider

Let

then

where

Then, by incorporating Eq.(4.8), Eq.(4.7) is transformed into the following form:

which implies that

where $\Theta_1 = \frac{1}{\rho_1^2} M_1 \Sigma_1 M_1^T$ with $\rho_1 = a_{22}a_{32}a_{41}\sigma_1$.

By tedious and complex computation and the incorporation of (4.5), we get that

By Lemma 3.1 of ^[20], $\Theta_1$ is positive definite. Then $\Sigma_1 = \rho_1^2 M_1^{-1} \Theta_1 (M_1^{-1})^T$ is also a positive definite matrix.

Step 2. For System (4.4), we consider

Let

where $\omega_{1} = a_{43}^2-a_{11}a_{43}+a_{13}a_{41}$.

Case 2.1. If $\omega_{1} \neq 0$, we find the standard matrix $M_{21}$ such that $B_1 = M_{21} T_{1} M_{21}^{-1}$, where

and $B_1$ is defined in Eq.(4.8). Similar to Step 1, Eq.(4.11) is transformed into the following equation:

which is the same with Eq.(4.10), and $\Theta_1$ is positive definite which implies that $\Sigma_2 = \rho_{21}^2(M_{21}M_2)^{-1} \Theta_1 (M_2^{-1}M_{21}^{-1})^T$ with $\rho_{21} = a_{32}|\omega_{1}|\sigma_2$ is also positive definite.

Case 2.2. If $\omega_{1} = 0$, there exists a new standard matrix $M_{22}$ such that $B_{2} = M_{22} T_1 M_{22}^{-1}$, where

Similarly, Eq (4.11) is transformed into the following form

which implies that

where $\Theta_2 = \frac{1}{\rho_{22}^2} (M_{22}M_2)\Sigma_{2}(M_{22}M_2)^T$ with $\rho_{22} = a_{13}a_{32}\sigma_2$. Direct computation induces that

By Lemma 3.2 of ^[20], $\Theta_2$ is semi-positive definite. Then $\Sigma_2 = \rho_{22}^2(M_{22}M_2)^{-1} \Theta_2 ((M_{22}M_2)^{-1})^T$ is also semi-positive definite.

Hence, for System (4.4), $\Sigma_2$ is positive definite for $\omega_{1}\neq 0$ or semi-positive definite for $\omega_{1} = 0$.

Step 3. For System (4.4), we consider the following algebraic equation:

Likewise, let

Next, we discuss them in three cases.

Case 3.1. If $\omega_1 \neq 0$, let $C_3 = P_3 T_2 {P_3}^{-1}$, where

with $\omega_2 = a_{13}a_{41}-a_{11}a_{43}+a_{22}a_{43}.$

Case 3.1.1. If $\omega_2 \neq 0$, there exists a standard matrix $M_{31}$ such that $B_1 = M_{31}C_3(M_{31}^{-1})$, where

with

and $B_1$ is defined in Eq.(4.8). Similar to Step 1, Eq.(4.13) is transformed into the following equation:

which is the same with Eq.(4.10), and $\Theta_1$ is positive definite which implies that $\Sigma_3 = \rho_{31}^2 (M_{31} P_3 M_3)^{-1} \Theta_1 ((M_{31} P_3 M_3)^{-1})^T$ with $\rho_{31} = a_{22}|\omega_2 |\sigma_{3} > 0$ is also positive definite.

Case 3.1.2. If $\omega_{2} = 0$, then

where

Direct computation yields that

Similarly, Eq.(4.13) is transformed into the following equation:

By Lemma 3.2 of ^[20], $\Theta_3$ is semi-positive definite, which implies that $\Sigma_3 = \rho_{32}^2(M_{32}P_{3}M_{3})^{-1} \Theta_2 ((M_{32}P_{3}M_{3})^{-1})^T$ with $\rho_{32} = \frac{a_{13}|\omega_{1}|}{a_{43}}\sigma_3 > 0$ is also semi-positive definite.

Case 3.2. If $\omega_1 = 0$, then

where

and

Similarly, Eq.(4.13) is transformed into the following form:

By Lemma 3.2 of ^[20], $\Theta_4$ is semi-positive definite, which implies that $\Sigma_3 = \rho_{33}^2(M_{33}M_{3})^{-1} \Theta_4 ((M_{33}M_{3})^{-1})^T$ with $\rho_{33} = a_{22}a_{43}\sigma_3 > 0$ is also semi-positive definite.

Hence, for System (4.4), $\Sigma_3$ is positive definite for $\omega_{1}\neq 0, \; \omega_{2}\neq 0$ and semi-positive definite for $\omega_{1}\neq 0$ and $\omega_{2} = 0$ or $\omega_{1} = 0$.

With all the things above, $\Sigma = \Sigma_1 + \Sigma_2 + \Sigma_3$ in Eq.(4.6) is positive definite. Thus, there is a local asymptotic density function $\Phi(\hat{x}, \hat{y}, \hat{v}, \hat{w})$ near the quasi-endemic equilibrium $\widetilde{E}^*$.

5.   Numerical simulations

In this section, we give some numerical simulations to illustrate our theoretical results.

5.1.   Determined system given by Eq.(1.4)

For the ODE system given by Eq.(1.4), the parameters refer to the data in ^[1] as follows

By Eq.(1.2), we obtain $R_1 = \frac{\beta c\alpha s^n}{\gamma(\alpha+\mu)(a^n\mu^n+s^n)}\approx 1.2664 > 1$; then, there exists an endemic equilibrium point $E^*(x^*, y^*, v^*, w^*)$, where $x^* = a\left(\frac{A}{\beta -A}\right)^\frac{1}{n}\approx 148.45, \; \; y^* = \frac{\gamma}{Ac\alpha}\left(s-\mu x^*\right)\approx 20.8763$, $v^* = \frac{c\alpha}{\gamma}y^*\approx 2009.3\; \; \text{and}\; \; w^* = \frac{\alpha+\mu}{\beta}y^*\approx 1503.1,$ with $A = \frac{(\mu+\alpha)\gamma}{c\alpha}\doteq 0.0037$. By Theorem 1.1(ii), the endemic equilibrium point $E^*$ is globally asymptotically stable, as illustrated in Figure 1 (red lines). By decreasing the value $c$ by $c = 400$, we obtain $R_1 = 0.9211 < 1$. By Theorem 1.1(i), the disease-free equilibrium point $E_0 (900, 0, 0, 0)$ is globally asymptotically stable, as illustrated in Figure 1 (blue lines).

5.2.   Stochastic system given by Eq.(1.6)

In the subsection, we consider the effect of white noises on the HIV system and establish the discretized system by Milstein method^[21]:

where $\xi_k, \; \zeta_k, \; \varsigma_k\; (k = 1, 2, \cdots)$ are independent Gaussian random variables, which satisfy the standard normal distribution $N(0, 1)$.

Fixing the same parameters as Eq.(5.1), we choose the noise intensities $\sigma_1 = 0.05, \; \sigma_2 = 0.08, \; \sigma_3 = 0.08$. By Eq.(2.5), we get that

Then by Theorem 3, there exists a stationary distribution for the degenerate system given by Eq.(1.6), which implies the persistence of the disease, as illustrated in the left graph of Figure 2. The right graph of Figure 2 demonstrates the distribution of a density function near the deterministic steady state.

Further increasing the noise intensities $\sigma_1\; = \; \sigma_2\; = \; \sigma_3\; = \; 0.3\;$ and choosing $\mu = 0.07, \; \beta = 0.005, \; c = 50, \; \alpha = 0.5, \; a = 100\; \text{and}\; \sigma = 0.9$, we obtain from Eq.(3.1) that $\bar{\mu} = -0.1715 < 0$. From Theorem 4, we know that the infected $CD4^+ T$ cell population $y(t)$ and virus particle $v(t)$ tend to extinction (see Figure 3).

6.   Discussion and conclusions

Some authors give the numerical simulations and theoretical analysis for an HIV infection model with $CD4^+$ T-cells. For example, Evirgen et al.^[22] gave the existence and uniqueness of the solutions for a fractionalized HIV infection model with the Atangana-Baleanu fractional derivative by using the Arzela-Ascoli theorem. Umar et al. ^[23] provided the numerical outcomes of a nonlinear HIV infection system, which is different from the Runge-Kutta method. Dewasurendra et al.^[24] applied the MDDiM method to an HIV infection model of $CD4^+$ T-cells, which shows the advantages over HAM.

Our main difference is the proposal of an HIV model with a Hill-type infection rate $\frac{x^n}{a^n+x^n}\; (n\geq 1)$ and distributed delay under the disturbance of white noise and proof of the existence of a stationary distribution by constructing a suitable Lyapunov function, which is a vast challenge. More importantly, we have given the exact local probability density function near the quasi-equilibrium by solving the corresponding Fokker-Planck equation. We have given the rigorous mathematical proof by describing the dynamics of the system, not only the numerical simulations.

In this paper, we first demonstrated the global asymptotical stability of the disease-free equilibrium and endemic equilibrium for the deterministic system. Second, the existence of a stationary distribution for an equivalent degenerate stochastic system was derived to obtain the sharp critical value $R_0^s$ by using the theory of Khasminskii. $R_0^s$ is consistent with the basic reproduction number without the white noises. Under a certain condition, the sufficient conditions for the extinction of the diseases have been given. In the part of numerical simulation, by incorporating the experimental data ^[25], we have applied the default parameter values given in Table 1 of ^[1] to verify the effectiveness of a stochastic system with degenerate diffusion.

There are still many interesting and instructive issues worthy of further study. For example, we consider the existence and uniqueness of the positive periodic solutions for the complex periodic system.

Acknowledgments

The work was supported by the Fundamental Research Funds for the Central Universities (No. 22CX03013A) and Shandong Provincial Natural Science Foundation (Nos. ZR2020MA039, ZR202102250288).

Conflict of interest

The authors declare that there is no conflicts of interest.