Stationary distribution for a three-dimensional stochastic viral infection model with general distributed delay

1.   Introduction

In the past few decades, there has been a lot of interest in mathematical models of viral dynamics and epidemic dynamics. Since viruses can directly reproduce inside of their hosts, a suitable model can shed light on the dynamics of the viral load population in vivo. In fact, by attacking infected cells, cytotoxic T lymphocytes (CTLs) play a crucial part in antiviral defense in the majority of virus infections. As a result, recent years have seen an enormous quantity of research into the population dynamics of viral infection with CTL response (see ^[1,2,3,4]). On the other hand, Bartholdy et al. ^[3] and Wodarz et al. ^[4] found that the turnover of free virus is much faster than that of infected cells, which allowed them to make a quasi-steady-state assumption, that is, the amount of free virus is simply proportional to the number of infected cells. In addition, the most basic models only consider the source of uninfected cells but ignore proliferation of the target cells. Therefore, a reasonable model for the population dynamics of target cells should take logistic proliferation term into consideration. Furthermore, in many biological models, time delay cannot be disregarded. A length of time $\tau$ may be required for antigenic stimulation to produce CTLs, and the CTL response at time $t$ may rely on the antigen population at time $t-\tau$. Xie et al. ^[4] present a model of delayed viral infection with immune response

where $x(t), y(t)$ and $z(t)$ represent the number of susceptible host cells, viral population and CTLs, respectively. At a rate of $\lambda$, susceptible host cells are generated, die at a rate of $dx$ and become infected by the virus at a rate of $\beta xy$. According to the lytic effector mechanisms of the CTL response, infected cells die at a rate of $ay$ and are killed by the CTL response at a rate of $pyz$. The CTL response occurs proportionally to the number of infected cells at a given time $cy(t-z)(t-z)$ and exponentially decays according to its level of activity $bz$. Additionally, the CTL response time delay is $\tau$.

The dynamical behavior of infectious diseases model with distributed delay has been studied by many researchers (see ^[5,6,7,8]) . Similar to ^[5] , in this paper, we will mainly consider the following viral infection model with general distribution delay

The delay kernel $F:[0, \infty)\rightarrow [0, \infty)$ takes the form $F(s) = \frac{s^n \alpha^{n+1} e^{-\alpha s}}{n!}$ for constant $\alpha > 0$ and integer $n\geq0$. The kernel with $n = 0$, i.e., $F(s) = \alpha \, e^{-\alpha s}$ is called the weak kernel which is the case to be considered in this paper.

However, in the real world, many unavoidable factors will affect the viral infection model. As a result, some authors added white noise to deterministic systems to demonstrate how environmental noise affects infectious disease population dynamics (see ^[9,10,11,12]). Linear perturbation, which is the simplest and most common assumption to introduce stochastic noise into deterministic models, is extensively used for species interactions and disease transmission. Here, we establish the stochastic infection model with distributed delay by taking into consideration the two factors mentioned above.

In our literature, we will consider weight function is weak kernel form. Let

Based on the linear chain technique, the equations for system (1.2) are transformed as follows

For the purpose of later analysis and comparison, we need to introduce the corresponding deterministic system of model (1.3), namely,

Using the similar method of Ma ^[13], the basic reproduction of system (1.4) can be expressed as $R_0 = \lambda\beta/ ad.$ If $R_0\leq 1,$ system (1.4) has an infection-free equilibrium $E_0 = (\frac{\lambda}{d}, 0, 0, 0)$ and is globally asymptotically stable. If $1 < R_0\leq 1+b\beta/cd,$ in addition to the infection-free equilibrium $E_0$, then system (1.4) has another unique equilibrium $E_1 = (\bar{x}, \bar{y}, \; \bar{z}, \; \bar{w}) = (\frac{a}{\beta}, \frac{\beta \lambda-ad}{a\beta}, 0, 0)$ and is globally asymptotically stable. If $R_0 > 1+\frac{b\beta}{cd}$, in addition $E_0$ and $E_1$, then system $(1.4)$ still has another unique infected equilibrium $E_2 = (x^+, \; y^+, \; z^+, \; w^+) = \bigg(\frac{c\lambda}{cd+b\beta}, \; \frac{b}{c}, \; \frac{c\beta\lambda-acd-ab\beta}{cdp+bp\beta}, \; \frac{b(c\beta\lambda-acd-ab\beta)}{c(cdp+bp\beta)}\bigg).$

We shall focus on the existence and uniqueness of a stable distribution of the positive solutions to model (1.3) in this paper. The stability of positive equilibrium state plays a key role in the study of the dynamical behavior of infectious disease systems. Compared with model (1.4), stochastic one (1.3) has no positive equilibrium to investigate its stability. Since stationary distribution means weak stability in stochastic sense, we focus on the existence of stationary distribution for model (1.3). The main effort is to construct the suitable Lyapunov function. As far as we comprehend, it is very challenging to create the proper Lyapunov function for system (1.3). This encourages us to work in this area. The remainder of this essay is structured as follows. The existence and uniqueness of a global beneficial solution to the system (1.3) are demonstrated in Section 2. In Section 3, several suitable Lyapunov functions are constructed to illustrate that the global solution of system (1.3) is stationary.

2.   Existence and uniqueness of a global positive solution

Theorem 2.1: For any initial value $(x(0), y(0), z(0), w(0))\in \mathbb{R }_+^4$, there is a unique solution $(x(t), y(t), z(t), w(t))$ of system (1.3) on $t\geq 0$ and the solution will remain in $\mathbb{R }_+^4$ with probability 1, i.e., $(x(t), y(t), z(t), w(t))\in \mathbb{R }_+^4$ for $t\geq 0$ almost surely $(a.s.).$

Proof. In light of the similarity with ^[14], the beginning of the proof is omitted. We only present the key stochastic Lyapunov function.

Define a $C^2$-function $Q(x, y, z, w)$ by

where $c_1, c_2, c_3$ are positive constant to be determined later. The nonnegativity of this function can be seen from

Using Itô's formula, we get

where

Let $c_1 = \frac{a}{\beta}, \, c_2 = \frac{b}{p},$ $0 < c_3\leq \min\{\frac{p}{\alpha}, \frac{c}{\alpha}\}$ such that $c_1\beta-a = 0$, $c_2p-b = 0$, $c_3\alpha-p\leq0$, $c-c_3\alpha\leq0.$ Then,

Obviously, $k_0$ is a positive constant which is independent of $x, y, z\, \, {\text{and }}\, w$. Hence, we omit the rest of the proof of Theorem 2.1 since it is mostly similar to Wang ^[14]. This completes the proof.

3.   Existence of stationary distribution

We need the following lemma to prove our main result. Consider the integral equation:

Lemma 3.1(^[15]). Suppose that the coefficients of (3.1) are independent of $t$ and satisfy the following conditions for some constant $B$:

in $D_\rho \in \mathbb{R }_+^d$ for every $\rho > 0,$ and that there exists a nonnegative $C^2$-function $V(x)$ in $\mathbb{R }_+^d$ such that

Then, system (3.1) has a solution, which is a stationary Markov process.

Here, we present a stationary distribution theorem. Define

where $\Lambda = c\bar{y}-(b+\frac{\sigma_3^2}{2}) > 0,$ $r = d\wedge a,$ and we denote $a\wedge b = \min\{a, \, b\}, \, a\vee b = \max\{a, \, b\}.$

Theorem 3.1. Assume $R^s > 0, \, d-\sigma_1^2 > 0$ and $r-\sigma_2^2 > 0.$ Then there exists a positive solution $(x(t), y(t), z(t), w(t))$ of system (1.3) which is a stationary Markov process.

Proof. We can substitute the global existence of the solutions of model (1.3) for condition (3.2) in Lemma 3.1, based on Remark 5 of Xu et al. ^[16]. We have established that system (1.3) has a global solution by Theorem 2.1. Thus condition (3.2) is satisfied. We simply need to confirm that condition (3.3) holds. This means that for any $(x, y, z, w)\in \mathbb{R }_+^4 \backslash D_\epsilon$, $LV(x, y, z, w)\leq-1$, we only need to construct a nonnegative $C^2$-function $V$ and a closed set $D_\epsilon$. As a convenience, we define

where $e_1$ is a positive constant to be determined later, $l = \frac{8r^2c^2}{(d-\sigma_1^2)(r-\sigma_2^2)\Lambda}$, $U_1 = \frac{(x-\bar{x})^2}{2}$, $U_2 = y-\bar{y}-\bar{y}\ln\frac{y}{\bar{y}}$, $U_3 = \frac{(x-\bar{x}+y-\bar{y})^2}{2}$, $Q_1 = -\ln z-\frac{e_1}{\alpha}\ln w$, $Q_2 = U_1+\frac{a}{\beta}U_2+\frac{(d-\sigma_1^2)(r-\sigma_2^2)}{2r^2}U_3$, $Q_3 = \frac{ap\bar{y}}{b\beta}z+\frac{ap\bar{y}\, c}{\alpha\beta b}w$.

Since $\lambda-d\bar{x} = \beta\bar{x}\, \overline{y} = a\bar{y}$, we apply Itô's formula to obtain

and

where $r = d\land a$, we also use the basic inequality $(a+b)^2\leq2(a^2+b^2)$ and Young inequality. It follows from (3.4)–(3.6) that

Making use of Itô's formula to $Q_3$ yields

Therefore,

In addition,

Letting $e_1 = c\cdot\bar{y}$, by virtue of Young inequality, one gets

Together with (3.7), this results in

where

Define

Then, we obtain

and

Define

where $\theta$ is a constant satisfying $0 < \theta < \min\{ \frac{d-\frac{\sigma_1^2}{2}}{d+\frac{\sigma_1^2}{2}}, \frac{a-\frac{\sigma_2^2}{2}}{a+\frac{\sigma_2^2}{2}}, \frac{b-\frac{\sigma_3^2}{2}}{b+\frac{\sigma_3^2}{2}}\}$. Then,

in which

Construct

where $M > 0$, satisfies

and

Note that $G$ is a continuous function and $\liminf\limits_{n\to\infty, (x, y, z, w)\in \mathbb{R } _+^4\backslash Q_n} G(x, y, z, w) = +\infty$, where $Q_n = (\frac{1}{n}, n)\times(\frac{1}{n}, n)\times(\frac{1}{n}, n)\times(\frac{1}{n}, n)$. Thus, $G(x, y, z, w)$ has a minimum point $(x_0, y_0, z_0, w_0)$ in the interior of $\mathbb{R }^4_+$. Define a nonnegative $C^2$ -function by

In view of (3.8)–(3.10) and (3.11), we get

One can easily see from (3.12) that, if $y\to0^+\, {\text{or}}\, \, z\to 0^+,$ then

if $x\to 0^+$ or $w\to 0^+$ or $x\to+\infty$ or $y\to+\infty$ or\, $z\to+\infty$ or $w\to+\infty$, then

In other words,

where $D_\epsilon = \{(x, y, z, w)\in \mathbb{R }_+^4:\epsilon\leq x \leq\frac{1}{\epsilon}, \epsilon\leq y < \frac{1}{\epsilon}, \epsilon\leq z\leq\frac{1}{\epsilon}, \epsilon^3\leq w\leq\frac{1}{\epsilon^3}\}$ and $\epsilon$ is a sufficiently small constant. The proof is completed.

Remark 3.1. In the proof of above theorem, the construction of $V_1$ is one of the difficulties. The term $Q_3$ in $V_1$ is is constructed to eliminate $\frac{ap\, \bar{y}}{\beta}\, z$ in $LQ_2$. The item $l(Q_2+Q_3)$ is used to eliminate $\frac{2c^2}{\Lambda}(y-\bar{y})^2$ in $LQ_1$.

Remark 3.2. From the expression of $R^s$, we can see that if there is no white noise, $R^s > 0$ is equivalent to $R_1 > 1+\frac{b\beta}{cd}$.

4.   Numerical simulations

Using the well-known numerical method of Milstein ^[17], we get the discretization equation for system (1.3)

where the time increment $\triangle t > 0$ and $\eta_{i, k} \, \, (i = 1, 2, 3)$ are three independent Gaussian random variables which follow the distribution $N(0, 1),$ equivalently, they come from the three independent from each other components of a three dimensional Wiener process with zero mean and variance $\triangle t,$ for $k = 1, 2, \cdots.$ According to Xie et al. ^[4], the corresponding biological parameters of system (1.3) are assumed: $\lambda = 255, \alpha = 1, d = 0.1, \beta = 0.002, a = 5, p = 0.1, c = 0.2, b = 0.1, r = d\land a = 0.1, \,$ $\sigma_1 = \sigma_2 = \sigma_3 = 0.0001.$ The initial condition is $(x_0, y_0, z_0, w_0) = (2600, 0.5, 0.5, 0.25).$ Then, we calculate that $R^s = 0.05 > 0.$ Based on Theorems 2.1 and 3.1, we can conclude that system (1.3) admits a global positive stationary solution on $\mathbb{R } _+^4$, see the left-hand figures of Figure 1 and the corresponding histograms of each population can be seen in right-hand column.

5.   Conclusions

In this paper, we consider a special kernel function $F(t) = \alpha e^{-\alpha t}$ to investigate the continuous delay effect on the population of stochastic viral infection systems. We derived the sufficient conditions for the existence of stationary distribution by constructing a suitable stochastic Lyapunov function. In addition, we only consider the effect of white noise on the dynamics of the viral infection system with distributed delays. It is interesting to consider the effect of Lévy jumps. Some researchers ^[18,19] studied the persistence and extinction of the stochastic systems with Lévy jumps. Furthermore, it should be noted that the system may be analytically solved by using the Lie algebra method ^[20,21]. In our further research, we will study the existence of a unique stationary distribution of the stochastic systems with distributed delay and Lévy jumps. Also, it may be possible to solve the stochastic model via the Lie algebra method.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgement

This work is supported by Hainan Provincial Natural Science Foundation of China (No.121RC554,122RC679) and Talent Program of Hainan Medical University (No. XRC2020030).

Conflicts of interest

No potential conflict of interest.