Threshold Dynamics of a Non-Linear Stochastic Viral Model with Time Delay and CTL Responsiveness

Abstract

This article focuses on a stochastic viral model with distributed delay and CTL responsiveness. It is shown that the viral disease will be extinct if the stochastic reproductive ratio is less than one. However, when the stochastic reproductive ratio is more than one, the viral infection system consists of an ergodic stationary distribution. Furthermore, we obtain the existence and uniqueness of the global positive solution by constructing a suitable Lyapunov function. Finally, we illustrate our results by numerical simulation.

1. Introduction

It is confirmed that approximately 100–250 million people are infected every year by different viruses, especially in regions of Asia and Africa [1]. To control epidemic viral diseases, an epidemic viral model is very useful, which can provide insights into the dynamics of viruses in vivo and offer a better understanding of viral diseases [2,3,4,5,6,7,8,9,10].

CTLs (Cytotoxic T Lymphocytes) play a significant role in antiviral mechanisms. On the one hand, CTLs imply the main immune factor inhibiting cell that limits the development of virus replication in vivo and depends on viral load [11,12,13]; on the other hand, it has recently been demonstrated that infected cells are killed not by the virus but by specific CTLs in some infectious diseases such as hepatitis B [3,14]. Therefore, the dynamics of the epidemic viral model with CTL responsiveness have drawn much attention from researchers in related areas [11,12,13,14,15,16].

The CTL immune response against a single pool of infected viral cells has been considered in [15,16,17], which is described by the following system:

$$\left\{\begin{array}{c}dx\left(t\right)=(\lambda -hx-exy)dt,\hfill \\ dy\left(t\right)=(exy-ay-pyz)dt,\hfill \\ dz\left(t\right)=\left(f\right(y,z)-bz)dt,\hfill \end{array}\right.$$

(1)

where the definitions of the variables are described in the following

Table 1. Variables in the model (adapted from [15,16]).

Variables	Definition
$x\left(t\right)$	the number of susceptible host cells at time t
$y\left(t\right)$	the virus population at time t
$z\left(t\right)$	the number of CTLs at time t
$\lambda$	susceptible host cells generation rate
$hx$	susceptible host cells death rate
$exy$	rate at which susceptible host cells become infected by the virus
$ay$	infected cells death rate
$pyz$	rate at which infected cells are killed by the CTL
$f(y,z)$	rate of immune response due to virus activation

.

Note that time-series data of the immune state of patients look rather irregular. The possibility of dynamics of infinite delay has been introduced into the equations used in mathematical biology models since Volterra [18,19,20,21,22] translated the cumulative effect of the history of a system. K. Wang et al. [13] incorporate a time delay of the immune response in the system (1) to obtain the following system:

$$\left\{\begin{array}{c}dx\left(t\right)=(\lambda -hx(t)-ex(t\left)y\right(t\left)\right)dt,\hfill \\ dy\left(t\right)=\left(ex\right(t\left)y\right(t)-ay(t)-py(t\left)z\right(t\left)\right)dt,\hfill \\ dz\left(t\right)=\left(cy\right(t-\tau )-bz(t\left)\right)dt,\hfill \end{array}\right.$$

(2)

They described the relationship between virus replication and the instantaneous immune response.

In the real world, epidemic models are also affected by escapable environmental white noise [5,6,7,8,9,23,24,25,26,27,28], full of randomness and stochasticity. As we understand it, there is little to no work researching the extinction and the ergodic stationary distribution of a system (2) with stochasticity and distributed delay, which is mainly due to the fact that Equation (2) is of a degenerate type. In this paper, we use an asymptotic approach [6,8,9,25] and give the following system:

$$\left\{\begin{array}{c}dx\left(t\right)=(\lambda -hx\left(t\right)-ex\left(t\right)y\left(t\right))dt-{\sigma }_{1}xd{B}_1\left(t\right),\hfill \\ dy\left(t\right)=(ex\left(t\right)y\left(t\right)-ay\left(t\right)-py\left(t\right)z\left(t\right))dt-{\sigma }_{2}yd{B}_2\left(t\right),\hfill \\ dz\left(t\right)=(cw\left(t\right)-bz\left(t\right))dt-{\sigma }_{3}zd{B}_3\left(t\right),\hfill \\ dw\left(t\right)=\sigma \left(y\right(t)-w(t\left)\right)dt,\hfill \end{array}\right.$$

(3)

The main purpose of this paper is to study the ergodic stationary distribution and extinction of the system (3). The existence and uniqueness of the global positive solution are also introduced.

The remainder of the paper is organized as follows. In Section 2, we introduce some necessary results throughout this paper. In Section 3, we show the uniqueness and positivity of global positive solutions of a stochastic system (3) with any positive initial value. In Section 4, we prove the existence and uniqueness of an ergodic stationary distribution of the solutions to the system (3) by constructing a suitable stochastic Lyapunov function; furthermore, we establish the persistence in the mean of the solutions of the system (3). In Section 5, we establish sufficient conditions for the extinction of the viral model. Some numerical simulations are introduced to demonstrate the theoretical results and reveal the effects of white noise. Finally, some concluding remarks and future directions are presented to close this paper.

2. Preliminaries

Throughout this paper, we first give some basic conceptions, as in [6]. Let $(\Omega ,\mathcal{F},$${\left\{{\mathcal{F}}_t\right\}}_{t\ge 0},\mathbb{P})$ be a complete probability space with a filtration ${\left\{{\mathcal{F}}_t\right\}}_{t\ge 0}$ satisfying the usual conditions (i.e., it is right continuous, and ${\mathcal{F}}_0$ contains all $\mathbb{P}-$null sets). Define ${\mathbb{R}}_{+}^4=\{x\in {\mathbb{R}}^4:x_i>0forall1\le i\le 4\}$, ${\overline{\mathbb{R}}}_{+}^4=\{x\in {\mathbb{R}}^4:x_i\ge 0forall1\le i\le 4\}$. In addition, if $g\left(t\right)$ is an integral function on $t\in [0,\infty )$, define $g^{\mu }=sup\{g\left(t\right)\mid t\ge 0\},$$g^l=inf\{g\left(t\right)\mid t\ge 0\}$.

Firstly, we consider the general 4-dimensional stochastic differential equation

$$\begin{array}{c}dx\left(t\right)=l(x\left(t\right),t)dt+q(x\left(t\right),t)dB\left(t\right),\mathit{for}\mathit{t}\ge {\mathit{t}}_{\mathit{0}}\hfill \end{array}$$

(4)

with initial value $x\left(t_0\right)=x_0\in {\mathbb{R}}^4$, where $B\left(t\right)$ denotes 4-dimensional standard Brownian motion defined on the above probability space. Define the differential operator L associated with Equation (4) by Mao [6] as

$$L=\frac{\partial }{\partial t}+\Sigma l_i(x,t)\frac{\partial }{\partial x_i}+\frac{1}{2}\Sigma {\left[q^T(x,t)q(x,t)\right]}_{ij}\frac{{\partial }^2}{\partial x_i\partial x_j}.$$

If L acts on a function $V\in C^{2,1}({\mathbb{R}}^4\times {\overline{\mathbb{R}}}_{+};{\overline{\mathbb{R}}}_{+})$, then

$$LV(x,t)=V_t(x,t)+V_x(x,t)+\frac{1}{2}trac\left[q^T(x,t)V_{xx}(x,t)q(x,t)\right],$$

where $V_t=\frac{\partial V}{\partial t},V_x=(\frac{\partial V}{\partial x_1},\dots ,\frac{\partial V}{\partial x_4})$ and $V_{xx}={\left(\frac{{\partial }^{2}V}{\partial x_i\partial x_j}\right)}_{4\times 4}$. Let $x\left(t\right)$ be a homogeneous Markov process in ${\mathbb{R}}^4$, which is described as the following stochastic differential equation by It${\widehat{o}}^{\prime }$s formula [6]:

$$dV(x\left(t\right),t)=LV(x\left(t\right),t)dt+V_x(x\left(t\right),t)r(x\left(t\right),t)dB\left(t\right).$$

The diffusion matrix is defined as follows:

$$A\left(x\right)=\left(a_{ij}\left(x\right)\right),a_{ij}=\sum _{r=1}^4{g}_r^i\left(x\right)g_r^j\left(x\right).$$

3. Existence and Uniqueness of the Global Positive Solution

To study the dynamical behaviors of a viral model, where the solution is global and positive, because the coefficients of the system (3) do not satisfy the linear growth, the solutions of the system (3) may explode at a finite time. In this section, we show that there is a unique global positive solution of the system (3) from the idea in [8], the main theorem as follows.

Theorem 1.

For any initial value $(x\left(0\right),y\left(0\right),z\left(0\right),w\left(0\right))\in {\mathbb{R}}_{+}^4$, there is a unique positive solution $\left(x\right(t),y(t),z(t),w(t\left)\right)$ of system (3) on $t\ge 0$, and the solution will remain in ${\mathbb{R}}_{+}^4$ with a probability of one; that is to say, $(x\left(t\right),y\left(t\right),z\left(t\right),w\left(t\right))\in {\mathbb{R}}_{+}^4$ for all $t\ge 0$, almost surely.

Proof. 

Our proof is based on the works of Mao et al. [9]. We know the coefficients of system (3) are locally Lipschitz continuous; thus, there is a unique solution $\left(x\right(t),y(t),z(t),w(t\left)\right)$ on $[0,{\tau }_0)$ for any initial value $(x\left(0\right),y\left(0\right),z\left(0\right),w\left(0\right))\in {\mathbb{R}}_{+}^4$, where ${\tau }_0$ is an explosion time. If ${\tau }_0=\infty$ a.s, we can determine that the local solution is global. Let $n_0$ be a sufficiently large positive number for every component of $\left(x\right(0),y(0),z(0),w(0\left)\right)$ lying in $[\frac{1}{n_0},n_0]$; for any $n\ge n_0$, the stopping time

$$\begin{array}{c}{\tau }_n=inf\{t\in [0,{\tau }_0):min\{x\left(t\right),y\left(t\right),z\left(t\right),w\left(t\right)\}\le \frac{1}{n}ormax\{x\left(t\right),y\left(t\right),z\left(t\right),w\left(t\right)\}\ge n\}.\hfill \end{array}$$

(5)

We set $inf\varphi =\infty$. Obviously, ${\tau }_n$ is increasing as $n\to \infty$. Set ${\tau }_{\infty }=\underset{n\to \infty }{lim}{\tau }_n$; thus, ${\tau }_{\infty }\le {\tau }_0$ a.s. If we can prove that ${\tau }_{\infty }=\infty$ a. s., then ${\tau }_0=\infty$, which implies $(x\left(t\right),y\left(t\right),z\left(t\right),w\left(t\right))\in {\mathbb{R}}_{+}^4$ a. s. for all $t\ge 0$. If ${\tau }_{\infty }<\infty$ a.s., it is the same. If we want to complete the proof, we should verify that ${\tau }_{\infty }=\infty$ a. s. If this assertion is false, there are two constants $T\ge 0$ and $\epsilon \in (0,1)$, such that

$$P\{{\tau }_{\infty }\le T\}\ge \epsilon .$$

In addition, there is an integer $n_1\ge n_0$, such that

$$P\{{\tau }_n\le T\}\ge \epsilon foralln\ge n_1.$$

We define a fundamental $C^2$-function $U:{\mathbb{R}}_{+}^4\to {\overline{\mathbb{R}}}_{+}$, which is

$$\begin{array}{c}U=(x\left(t\right)-\alpha -\alpha ln\frac{x\left(t\right)}{\alpha })+(y\left(t\right)-\beta -\beta ln\frac{y\left(t\right)}{\beta })+(z\left(t\right)-1-lnz\left(t\right))+\gamma (w\left(t\right)-1-lnw\left(t\right)),\hfill \end{array}$$

(6)

where $\alpha ,\beta ,\gamma$ are positive constants, which will be determined in the following text. The non-negativity of the function U can be seen from $x-1-lnx\ge 0$ for any $x>0$.

Applying It${\widehat{o}}^{\prime }$s formula [6], we obtain

$$\begin{array}{c}dU(x,y,z,w)=LUdt-{\sigma }_1(x\left(t\right)-\alpha )d{B}_1\left(t\right)-{\sigma }_2(y\left(t\right)-\beta )d{B}_2\left(t\right)-{\sigma }_3(z\left(t\right)-1)d{B}_3\left(t\right),\hfill \end{array}$$

(7)

where

$$\begin{array}{cc}LU(x,y,z,w)\hfill & =(1-\frac{\alpha }{x})dx+(1-\frac{\beta }{y})dy+(1-\frac{1}{z})dz+\gamma (1-\frac{1}{w})dw+\frac{1}{2}(\alpha {\sigma }_1^2+\beta {\sigma }_2^2+{\sigma }_3^2\hfill \\ & =\lambda -dx-exy-\frac{\alpha \lambda }{x}+d\alpha +\alpha ey+exy-ay-pyz-\beta ex\hfill \\ & +\beta a+\beta pz+cw-bz-\frac{cw}{z}+b+\gamma \sigma y-\gamma \sigma w-\frac{\gamma \sigma y}{w}+\gamma \sigma +\frac{1}{2}(\alpha {\sigma }_1^2+\beta {\sigma }_2^2+{\sigma }_3^2.\hfill \end{array}$$

(8)

Choosing

$$\alpha =\frac{a-c}{e},\beta =\frac{b}{p},\gamma =\frac{c}{\sigma },$$

such that $\alpha e+\gamma \sigma =a,\beta p=b,\gamma \sigma =c$, we can then obtain

$$\begin{array}{ccc}LU(x,y,z,w)\le \lambda +ad+\beta a+b+\gamma \sigma \le \lambda +ad+\frac{ab}{p}+b+c+\frac{1}{2}(\frac{a-c}{e}{\sigma }_1^2+\frac{b}{p}{\sigma }_2^2+{\sigma }_3^2)\hfill & :=\hfill & K,\hfill \end{array}$$

(9)

where K is a positive constant. The remainder of the proof is similar to Theorem 3.1 in Mao [9]. Hence, we omit it here. □

4. Ergodic Stationary Distribution

Here, we present some theories about the stationary distribution in this section. Although there is no endemic equilibrium point of the stochastic system (3), we want to obtain the existence of an ergodic stationary distribution, which indicates the persistence of the disease. Firstly, we define $R_0^{*}$ as a stochastic reproductive ratio of the system (3), such as

$$R_0^{*}=\frac{\lambda ea}{{(a+\frac{1}{2}{\sigma }_2^2)}^2(d+\frac{1}{2}{\sigma }_1^2)},$$

which is equal to $R_0=\frac{\lambda e}{ad}$ when ${\sigma }_1={\sigma }_2=0$ [13]. Some known results about the theory of Has’Minskii are found in [28].

Lemma 1

([28]). The Markov process $X\left(t\right)$ has a unique ergodic stationary distribution $\mu (·)$ if there exists a bounded domain $U\subset E_l$ with regular boundary Γ, and

(A.1)

there is a positive number M such that ${\sum }_{i,j=1}^l{a}_{ij}\left(x\right){\xi }_i{\xi }_j\ge M{\left|\xi \right|}^2,x\in U,\xi \in R^l.$

(A.2)

there exists a non-negative $C^2$ function V, such that $LV$ is negative for any $E_l\backslash U$. Then,

$$P_x\left\{\underset{T\to \infty }{lim}\frac{1}{T}{\int }_0^{T}f\left(X\left(t\right)\right)dt={\int }_{E_l}f\left(x\right)\mu \left(dx\right)\right\}=1,$$

for all $x\in E_l$, where $f(·)$ is a function integrable with respect to the measure μ.

Based on the theory of Has’minskii [28], we will give conditions which guarantee the existence of an ergodic stationary distribution.

Theorem 2.

Assume that $R_0^{*}>1$, then the solution $\left(x\right(t),y(t),z(t),w(t\left)\right)$ of system (3) has an ergodic unique stationary distribution for any initial value $(x\left(0\right),y\left(0\right),z\left(0\right),w\left(0\right))\in {\mathbb{R}}_{+}^4.$

Proof. 

The proof of Theorem 2 should satisfy the conditions of Lemma 1. Verify that $(A.1)$ holds. Apparently, the corresponding diffusion matrix of system (3) is given by

$$A=\left(\begin{array}{cccc}{\sigma }_1^2{x}^2& 0& 0& 0\\ 0& {\sigma }_2^2{y}^2& 0& 0\\ 0& 0& {\sigma }_3^2{z}^2& 0\\ 0& 0& 0& w^2\end{array}\right).$$

Choosing $\tilde{M}=\underset{(x,y,z,w)\in D_{\epsilon }}{min}\{{\sigma }_1^2{x}^2,{\sigma }_2^2{y}^2,{\sigma }_3^2{z}^2,w^2\}>0$, we obtain

$$\sum _{i,j=1}^4{a}_{ij}(x,y,z,w){\xi }_i{\xi }_j={\sigma }_1^2{x}^2{\xi }_1^2+{\sigma }_2^2{y}^2{\xi }_2^2+{\sigma }_3^2{z}^2{\xi }_3^2+w^2{\xi }_4^2\ge \tilde{M}\mid \xi {\mid }^2,$$

for all $(x,y,z,w)\in D_{\epsilon },\xi =({\xi }_1,{\xi }_2,{\xi }_3,{\xi }_4)\in {\mathbb{R}}_{+}^4$, where

$$D_{\epsilon }=\{(x,y,z,w)\in {\mathbb{R}}_{+}^4:\epsilon <x<\frac{1}{\epsilon },\epsilon <y<\frac{1}{\epsilon },\epsilon <z<\frac{1}{\epsilon },\epsilon <w<\frac{1}{\epsilon }\}$$

is a bounded closed set, and $\epsilon >0$ is a sufficiently small number. Thus, condition $(A.1)$ is completed.

Now, we construct a $C^2$-function $V:{\mathbb{R}}_{+}^4\to \mathbb{R}$ as follows:

$$V(x,y,z,w)=M(-lnx-c_{0}lny)-c_{1}lny-lnz-lnw+\frac{1}{\theta +1}{(x+y+\frac{a}{4c}z+\frac{a}{2\sigma }w)}^{\theta +1},$$

where $\theta \in (0,min\frac{a}{a+{\sigma }_2^2},\frac{b}{b+{\sigma }_3^2},\frac{d}{d+{\sigma }_1^2})$ is a positive constant, and $c_1\le -\frac{2+b+\sigma +\frac{1}{2}({\sigma }_2^2+{\sigma }_3^2)}{a}.$

When $R_0^{*}>1$, the constant M satisfies the following condition:

$$0<M\le \frac{-2}{(d+\frac{1}{2}{\sigma }_1^2)(1-R_0^{*})}.$$

Applying It${\widehat{o}}^{\prime }$s formula to the function $V(x,y,z,w)$, denote

$$V_1=-lnx-c_{0}lny,V_2=-c_{1}lny,$$

$$V_3=-lnz,V_4=-lnw,$$

$$V_5=\frac{1}{\theta +1}{(x+y+\frac{a}{4c}z+\frac{a}{2\sigma }w)}^{\theta +1}.$$

We can apply the differential operator L to the above functions, respectively,

$$\begin{array}{cc}L{V}_1\hfill & =-\frac{1}{x}{x}^{\prime }-c_0\frac{1}{y}{y}^{\prime }+\frac{1}{2}({\sigma }_1^2+c_0{\sigma }_2^2)\hfill \\ & =-\frac{\lambda }{x}+d+ey-c_{0}ex+c_{0}a+c_{0}pz+\frac{1}{2}({\sigma }_1^2+c_0{\sigma }_2^2)\hfill \\ & \le -2\sqrt{\lambda c_{0}e}+c_{0}a+\frac{1}{2}{c}_0{\sigma }_2^2+d+\frac{1}{2}{\sigma }_1^2+ey+c_{0}pz.\hfill \end{array}$$

(10)

Supposing

$$g\left(c_0\right)=-2\sqrt{\lambda c_{0}e}+c_{0}a+\frac{1}{2}{c}_0{\sigma }_2^2,$$

$$g^{\prime }\left(c_0\right)=-\frac{\lambda e}{\sqrt{\lambda c_{0}e}}+a+\frac{1}{2}{\sigma }_2^2=0,$$

we can obtain

$$c_0=\frac{\lambda e}{{(a+\frac{1}{2}{\sigma }_2^2)}^2}.$$

Hence,

$$\begin{array}{cc}L{V}_1\hfill & =g\left(c_0\right)+d+\frac{1}{2}{\sigma }_1^2+ey+c_{0}pz\hfill \\ & \le -\frac{\lambda ea}{{(a+\frac{1}{2}{\sigma }_2^2)}^2}+d+\frac{1}{2}{\sigma }_1^2+ey+c_{0}pz\hfill \\ & =(d+\frac{1}{2}{\sigma }_1^2)(1-\frac{\lambda ea}{{(a+\frac{1}{2}{\sigma }_2^2)}^2(d+\frac{1}{2}{\sigma }_1^2)})+ey+\frac{\lambda ep}{{(a+\frac{1}{2}{\sigma }_2^2)}^2}z\hfill \\ & =(d+\frac{1}{2}{\sigma }_1^2)(1-R_0^{*})+ey+\frac{\lambda ep}{{(a+\frac{1}{2}{\sigma }_2^2)}^2}z;\hfill \end{array}$$

(11)

$$L{V}_2=-c_1\frac{1}{y}{y}^{\prime }+\frac{1}{2}{\sigma }_1^2=-c_{1}ex+c_{1}a+c_{1}pz+\frac{1}{2}{\sigma }_2^2;$$

(12)

$$L{V}_3=-\frac{1}{z}{z}^{\prime }+\frac{1}{2}{\sigma }_3^2=-\frac{cw}{z}+b+\frac{1}{2}{\sigma }_3^2;$$

(13)

$$L{V}_4=-\frac{1}{w}{w}^{\prime }=-\frac{\sigma y}{w}+b+\sigma ;$$

(14)

$$\begin{array}{cc}L{V}_5\hfill & ={(x+y+\frac{a}{4c}z+\frac{a}{2\sigma }w)}^{\theta }(x^{\prime }+y^{\prime }+\frac{a}{4c}{z}^{\prime }+\frac{a}{2\sigma }{w}^{\prime })\hfill \\ & +\frac{\theta }{2}{(x+y+\frac{a}{4c}z+\frac{a}{2\sigma }w)}^{\theta -1}({\sigma }_1^2{x}^2+{\sigma }_2^2{y}^2+{\left(\frac{a}{4c}\right)}^2{\sigma }_3^2{z}^2)\hfill \\ & \le {(x+y+\frac{a}{4c}z+\frac{a}{2\sigma }w)}^{\theta }(\lambda -dx-\frac{a}{2}y-\frac{a}{4}w-\frac{ab}{4c}z)\hfill \\ & +\frac{\theta }{2}{(x+y+\frac{a}{4c}z+\frac{a}{2\sigma }w)}^{\theta -1}({\sigma }_1^2{x}^2+{\sigma }_2^2{y}^2+{\left(\frac{a}{4c}\right)}^2{\sigma }_3^2{z}^2)\hfill \\ & \le \lambda {(x+y+\frac{a}{4c}z+\frac{a}{2\sigma }w)}^{\theta }-(dx+\frac{a}{2}y+\frac{ab}{4c}z+\frac{a}{4}w)(x^{\theta }+y^{\theta }+{\left(\frac{a}{4c}\right)}^{\theta }{z}^{\theta }+{\left(\frac{a}{2\sigma }\right)}^{\theta }{w}^{\theta })\hfill \\ & +\frac{\theta }{2}{(x+y+\frac{a}{4c}z+\frac{a}{2\sigma }w)}^{\theta -1}({\sigma }_1^2{x}^2+{\sigma }_2^2{y}^2+{\left(\frac{a}{4c}\right)}^2{\sigma }_3^2{z}^2)\hfill \\ & \le \lambda {(x+y+\frac{a}{4c}z+\frac{a}{2\sigma }w)}^{\theta }-d{x}^{\theta +1}-\frac{a}{2}{y}^{\theta +1}-{\left(\frac{a}{4c}\right)}^{\theta +1}b{z}^{\theta +1}-\frac{a}{4}{\left(\frac{a}{2\sigma }\right)}^{\theta }{w}^{\theta +1}\hfill \\ & +\frac{\theta }{2}{(x+y+\frac{a}{4c}z+\frac{a}{2\sigma }w)}^{\theta -1}({\sigma }_1^2{x}^2+{\sigma }_2^2{y}^2+{\left(\frac{a}{4c}\right)}^2{\sigma }_3^2{z}^2)\hfill \\ & \le B-d\theta x^{\theta +1}-\frac{a}{2}\theta y^{\theta +1}-{\left(\frac{a}{4c}\right)}^{\theta +1}b\theta z^{\theta +1}-\frac{a^{\theta +1}}{2^{\theta +3}{\sigma }^{\theta }}{w}^{\theta +1},\hfill \end{array}$$

(15)

where

$$\begin{gathered} B=sup\{\lambda {(x+y+\frac{a}{4c}z+\frac{a}{2\sigma }w)}^{\theta }+\frac{\theta }{2}{(x+y+\frac{a}{4c}z+\frac{a}{2\sigma }w)}^{\theta -1}({\sigma }_1^2{x}^2+{\sigma }_2^2{y}^2+{(\frac{a}{4c})}^2{\sigma }_3^2{z}^2), \\ -d(1-\theta )x^{\theta +1}-\frac{a}{2}(1-\theta )y^{\theta +1}-{(\frac{a}{4c})}^{\theta +1}b(1-\theta )z^{\theta +1}-\frac{a^{\theta +1}}{2^{\theta +4}{\sigma }^{\theta }}{w}^{\theta +1}\}, \end{gathered}$$

where $\theta \in (0,min\frac{a}{a+{\sigma }_2^2},\frac{b}{b+{\sigma }_3^2},\frac{d}{d+{\sigma }_1^2})$ is a positive constant, and $B<0$. From the above analysis, we have

$$\begin{array}{cc}LV(x,y,z,w)\hfill & =ML{V}_1+L{V}_2+L{V}_3+L{V}_4+L{V}_5\hfill \\ & \le M((d+\frac{1}{2}{\sigma }_1^2)(1-\frac{\lambda ea}{{(a+\frac{1}{2}{\sigma }_2^2)}^2(d+\frac{1}{2}{\sigma }_1^2)})+ey+\frac{\lambda ep}{{(a+\frac{1}{2}{\sigma }_2^2)}^2}z)\hfill \\ & -c_{1}ex+c_{1}a+c_{1}pz+\frac{1}{2}{\sigma }_2^2\hfill \\ & -\frac{1}{z}{z}^{\prime }+\frac{1}{2}{\sigma }_3^2=-\frac{cw}{z}+b+\frac{1}{2}{\sigma }_3^2-\frac{1}{w}{w}^{\prime }=-\frac{\sigma y}{w}+b+\sigma \hfill \\ & +B-d\theta x^{\theta +1}-\frac{a}{2}\theta y^{\theta +1}-{\left(\frac{a}{4c}\right)}^{\theta +1}b\theta z^{\theta +1}-\frac{a^{\theta +1}}{2^{\theta +3}{\sigma }^{\theta }}{w}^{\theta +1},\hfill \end{array}$$

(16)

and we define

$$f_1\left(x\right)=c_{1}a+b+\sigma +\frac{1}{2}({\sigma }_2^2+{\sigma }_3^2)-c_{1}ex-d\theta x^{\theta +1},$$

$$f_2\left(y\right)=M((d+\frac{1}{2}{\sigma }_1^2)(1-\frac{\lambda ea}{{(a+\frac{1}{2}{\sigma }_2^2)}^2(d+\frac{1}{2}{\sigma }_1^2)})+ey)-\frac{a}{2}\theta y^{\theta +1},$$

$$f_3\left(z\right)=M\frac{\lambda ep}{{(a+\frac{1}{2}{\sigma }_2^2)}^2}z+c_{1}pz-\frac{cw}{z}-{\left(\frac{a}{4c}\right)}^{\theta +1}b\theta z^{\theta +1},$$

$$f_3\left(w\right)=B-\frac{\sigma y}{w}-\frac{a^{\theta +1}}{2^{\theta +3}{\sigma }^{\theta }}{w}^{\theta +1}.$$

We can divide ${\mathbb{R}}_{+}^4\backslash D_{\epsilon }$ into the following eight domains:

$$D_1=\{(x,y,z,w)\in {\mathbb{R}}_{+}^4:0<x<\epsilon \};D_2=\{(x,y,z,w)\in {\mathbb{R}}_{+}^4:0<y<\epsilon \};$$

$$D_3=\{(x,y,z,w)\in {\mathbb{R}}_{+}^4:0<z<\epsilon ,\epsilon <w<\frac{1}{\epsilon }\};D_4=\{(x,y,z,w)\in {\mathbb{R}}_{+}^4:\epsilon <z<\frac{1}{\epsilon },0<w<\epsilon \};$$

$$D_5=\{(x,y,z,w)\in {\mathbb{R}}_{+}^4:x>\frac{1}{\epsilon }\};D_6=\{(x,y,z,w)\in {\mathbb{R}}_{+}^4:y>\frac{1}{\epsilon }\};$$

$$D_7=\{(x,y,z,w)\in {\mathbb{R}}_{+}^4:z>\frac{1}{\epsilon }\};D_8=\{(x,y,z,w)\in {\mathbb{R}}_{+}^4:w>\frac{1}{\epsilon }\}.$$

Clearly, $D_{\epsilon }={\bigcup }_{j=1}^8{D}_j$. In the following text, we will show that $LV(x,y,z,w)\le -1$ on ${\mathbb{R}}_{+}^4\backslash D_{\epsilon }$, which is equivalent to prove it on the above eight domains.

Case 1. If $(x,y,z,w)\in D_1$, one can choose

$$c_1\le -\frac{2+b+\sigma +\frac{1}{2}({\sigma }_2^2+{\sigma }_3^2)}{a},$$

and

$$LV(x,y,z,w)\le c_{1}a+b+\sigma +\frac{1}{2}({\sigma }_2^2+{\sigma }_3^2)-c_{1}ex-d\theta x^{\theta +1}\le -2;$$

Case 2. If $(x,y,z,w)\in D_2$, one can choose $R_0^{*}>1$ and

$$0<M\le \frac{-2}{(d+\frac{1}{2}{\sigma }_1^2)(1-R_0^{*})},$$

$$LV(x,y,z,w)\le M((d+\frac{1}{2}{\sigma }_1^2)(1-R_0^{*})+ey)-\frac{a}{2}\theta y^{\theta +1}\le -2;$$

Case 3. If $(x,y,z,w)\in D_3$,

$$LV(x,y,z,w)\le M\frac{\lambda ep}{{(a+\frac{1}{2}{\sigma }_2^2)}^2}z+c_{1}pz-\frac{cw}{z}-{(\frac{a}{4c})}^{\theta +1}b\theta z^{\theta +1}\le -\frac{cw}{z}\le -2;$$

Case 4. If $(x,y,z,w)\in D_4$,

$$LV(x,y,z,w)\le B-\frac{\sigma y}{w}-\frac{a^{\theta +1}}{2^{\theta +3}{\sigma }^{\theta }}{w}^{\theta +1}\le -\frac{\sigma y}{w}\le -2;$$

Case 5. If $(x,y,z,w)\in D_5$,

$$LV(x,y,z,w)\le c_{1}a+b+\sigma +\frac{1}{2}({\sigma }_2^2+{\sigma }_3^2)-c_{1}ex-d\theta x^{\theta +1}\le -c_{1}ex-d\theta x^{\theta +1}\le -2;$$

Case 6. If $(x,y,z,w)\in D_6$,

$$\begin{array}{cc}LV(x,y,z,w)\hfill & \le M((d+\frac{1}{2}{\sigma }_1^2)(1-R_0^{*})+ey)-\frac{a}{2}\theta y^{\theta +1}\hfill \\ & \le M((d+\frac{1}{2}{\sigma }_1^2)(1-R_0^{*})+ey)-\frac{a}{2}\theta {\left(\frac{1}{\epsilon }\right)}^{\theta +1}\hfill \\ & \le -2;\hfill \end{array}$$

(17)

Case 7. If $(x,y,z,w)\in D_7$, because $M\frac{\lambda ep}{{(a+\frac{1}{2}{\sigma }_2^2)}^2}+c_{1}p-{\left(\frac{a}{4c}\right)}^{\theta +1}b\theta {\left(\frac{1}{\epsilon }\right)}^{\theta }<0$, we determine that

$$\begin{array}{cc}LV(x,y,z,w)\hfill & \le M\frac{\lambda ep}{{(a+\frac{1}{2}{\sigma }_2^2)}^2}z+c_{1}pz-\frac{cw}{z}-{\left(\frac{a}{4c}\right)}^{\theta +1}b\theta z^{\theta +1}\hfill \\ & \le z(M\frac{\lambda ep}{{(a+\frac{1}{2}{\sigma }_2^2)}^2}+c_{1}p)-cw\epsilon -{\left(\frac{a}{4c}\right)}^{\theta +1}b\theta z^{\theta +1}\hfill \\ & \le z(M\frac{\lambda ep}{{(a+\frac{1}{2}{\sigma }_2^2)}^2}+c_{1}p-{\left(\frac{a}{4c}\right)}^{\theta +1}b\theta {\left(\frac{1}{\epsilon }\right)}^{\theta })\hfill \\ & \le -\frac{1}{\epsilon }\hfill \\ & \le -2;\hfill \end{array}$$

(18)

Case 8. If $(x,y,z,w)\in D_8$, because $B<0$, we determine that

$$LV(x,y,z,w)\le B-\frac{\sigma y}{w}-\frac{a^{\theta +1}}{2^{\theta +3}{\sigma }^{\theta }}{w}^{\theta +1}\le B-\frac{a^{\theta +1}}{2^{\theta +3}{\sigma }^{\theta }}{(\frac{1}{\epsilon })}^{\theta +1}\le -2.$$

Therefore, for all $(x,y,z,w)\in {\mathbb{R}}_{+}^4\backslash D_{\epsilon }$, $V(x,y,z,w)\le -1,$ which indicates that assumption $(A.2)$ holds.

We can know that the system (3) is ergodic and has a unique stationary distribution. This completes the proof. □

5. Extinction of the System 3

For the dynamical behavior of epidemic viral models, the main concern is finding the condition in which the virus will be eradicated in a long time when $R_0^{*}>1$. In this section, we shall consider the extinction of the system (3).

According to the results in [9], we can obtain the following lemma.

Lemma 2.

For any initial value, the solution of the stochastic model satisfies

$$\underset{t\to \infty }{lim}\frac{lnx\left(t\right)}{t}\le 0,\underset{t\to \infty }{lim}\frac{lny\left(t\right)}{t}\le 0,\underset{t\to \infty }{lim}\frac{lnz\left(t\right)}{t}\le 0,\underset{t\to \infty }{lim}\frac{lnw\left(t\right)}{t}\le 0a.s.$$

(19)

$$\underset{t\to \infty }{lim}\frac{x\left(t\right)+y\left(t\right)+z\left(t\right)+w\left(t\right)}{t}=0,a.s..$$

(20)

Moreover,

$$\underset{t\to 0}{lim}\frac{1}{t}{\int }_0^{t}x\left(m\right)d{B}_1\left(m\right)=0,\underset{t\to 0}{lim}\frac{1}{t}{\int }_0^{t}y\left(m\right)d{B}_2\left(m\right)=0,\underset{t\to 0}{lim}\frac{1}{t}{\int }_0^{t}z\left(m\right)d{B}_3\left(m\right)=0a.s..$$

(21)

Theorem 3.

Let $(x,y,z,w)$ be the solution of system (3) with any initial value $(x\left(0\right),y\left(0\right),z\left(0\right),w\left(0\right))\in {\mathbb{R}}_{+}^4$. If $R_0^{*}<1$, then the solution $(x,y,z,w)$ of system (3) satisfies

$$\underset{t\to \infty }{lim\; sup}\frac{lny\left(t\right)}{t}\le (a+\frac{1}{2}{\sigma }_2^2)(R_0^{*}-1)<0a.s..$$

Namely, the disease will be eradicated in the long term.

Proof. 

Applying Itô’s formula to $lny\left(t\right)$, we obtain

$$\begin{array}{cc}dlny\left(t\right)\hfill & =(ex-a-pz-\frac{1}{2}{\sigma }_2^2)dt+{\sigma }_{2}d{B}_2\left(t\right)\hfill \\ & \le [ex-(a+\frac{1}{2}{\sigma }_2^2)]dt+{\sigma }_{2}d{B}_2\left(t\right).\hfill \end{array}$$

(22)

Integrating the above formula from 0 to t on both sides, then

$$lny\left(t\right)-lny\left(0\right)\le {\int }_0^t[ex-(a+\frac{1}{2}{\sigma }_2^2)]ds+{\int }_0^t{\sigma }_{2}d{B}_2\left(s\right).$$

According to the strong law of large numbers [29], we have

$$\underset{t\to 0}{lim}\frac{1}{t}{\int }_0^{t}d{B}_2\left(s\right)=0a.s..$$

For,

$$\begin{array}{cc}dlnx\left(t\right)\hfill & =\frac{1}{x}(\lambda -dx\left(t\right)-ex\left(t\right)y\left(t\right)-x{\sigma }_1^2)dt+{\sigma }_{1}d{B}_1\left(t\right)\hfill \\ & \le [\frac{\lambda }{x}-(d+\frac{1}{2}{\sigma }_1^2)]dt+{\sigma }_{1}d{B}_1\left(t\right).\hfill \end{array}$$

(23)

From (23), we can obtain

$$dlnx\left(t\right)\ge -(d+\frac{1}{2}{\sigma }_1^2)dt,$$

and

$$\begin{array}{c}\frac{1}{dlnx\left(t\right)}\le -\frac{1}{(d+\frac{1}{2}{\sigma }_1^2)dt}.\hfill \end{array}$$

(24)

Meanwhile,

$$\frac{\lambda }{x}\ge dlnx\left(t\right)+(d+\frac{1}{2}{\sigma }_1^2)+{\sigma }_{1}d{B}_1\left(t\right),$$

which indicates that

$$\begin{array}{c}x\le \frac{\lambda }{dlnx\left(t\right)+(d+\frac{1}{2}{\sigma }_1^2)+{\sigma }_{1}d{B}_1\left(t\right)}\le \frac{\lambda }{d+\frac{1}{2}{\sigma }_1^2}+\frac{\lambda }{dlnx\left(t\right)+{\sigma }_{1}d{B}_1\left(t\right)}.\hfill \end{array}$$

(25)

There exists a small number ${\sigma }_2$, such that

$$\begin{array}{cc}{\int }_0^{t}x\left(s\right)ds\hfill & \le \frac{\lambda }{d+\frac{1}{2}{\sigma }_1^2}t+{\int }_0^t\frac{\lambda }{dlnx\left(s\right)+{\sigma }_{1}d{B}_1\left(s\right)}ds\hfill \\ & \le \frac{\lambda }{d+\frac{1}{2}{\sigma }_1^2}t+\lambda {\int }_0^t\frac{1}{dlnx\left(s\right)}ds\hfill \\ & \le \frac{\lambda }{d+\frac{1}{2}{\sigma }_1^2}t-\lambda {\int }_0^t\frac{1}{(d+\frac{1}{2}{\sigma }_1^2)}ds\hfill \\ & \le \frac{\lambda }{d+\frac{1}{2}{\sigma }_1^2}t\hfill \\ & \le \frac{\lambda a}{(d+\frac{1}{2}{\sigma }_1^2)(a+\frac{1}{2}{\sigma }_2^2)}t.\hfill \end{array}$$

(26)

Taking the superior limit and using the stochastic comparison theorem, combining (26), we obtain

$$\begin{array}{cc}\underset{t\to \infty }{lim\; sup}\frac{lny\left(t\right)}{t}\hfill & =\underset{t\to \infty }{lim\; sup}e\frac{1}{t}{\int }_0^{t}x\left(s\right)ds-(a+\frac{1}{2}{\sigma }_2^2)\hfill \\ & \le e\frac{\lambda a}{(d+\frac{1}{2}{\sigma }_1^2)(a+\frac{1}{2}{\sigma }_2^2)}-(a+\frac{1}{2}{\sigma }_2^2)\hfill \\ & =(a+\frac{1}{2}{\sigma }_2^2)(\frac{\lambda ea}{{(a+\frac{1}{2}{\sigma }_2^2)}^2(d+\frac{1}{2}{\sigma }_1^2)}-1)\hfill \\ & =(a+\frac{1}{2}{\sigma }_2^2)(R_0^{*}-1)\hfill \\ & <0a.s..\hfill \end{array}$$

(27)

Therefore, this indicates that

$$\underset{t\to \infty }{lim}y\left(t\right)=0a.s..$$

Consequently, this means that the virus will be eradicated in a long time. This completes the proof. □

6. Examples and Numerical Simulations

In this section, we will introduce some examples and numerical simulations to demonstrate the above theoretical results. Using the Milstein higher-order method developed in [23], we obtain the discretization equation of the system (3).

$$\left\{\begin{array}{c}x(k+1)=x\left(k\right)+[\lambda -dx\left(k\right)-ex\left(k\right)y\left(k\right)]▵t+{\sigma }_{1}x\left(k\right)\sqrt{▵t}{\xi }_k+\frac{{\sigma }_1^{2}x\left(k\right)}{2}▵t({\xi }_k^2-1),\hfill \\ y(k+1)=y\left(k\right)+[ex\left(k\right)y\left(k\right)-ay\left(k\right)-py\left(k\right)z\left(k\right)]▵t+{\sigma }_{2}y\left(k\right)\sqrt{▵t}{\eta }_k+\frac{{\sigma }_2^{2}y\left(k\right)}{2}▵t({\eta }_k^2-1),\hfill \\ z(k+1)=z\left(k\right)+[cw\left(k\right)-bz\left(k\right)]▵t+{\sigma }_{3}z\left(k\right)\sqrt{▵t}{\zeta }_k+\frac{{\sigma }_3^{2}z\left(k\right)}{2}▵t({\zeta }_k^2-1),\hfill \\ w(k+1)=w\left(k\right)+\sigma \left[y\right(k)-w(k\left)\right]▵t.\hfill \end{array}\right.$$

(28)

where the time increment $▵t$ is positive, and ${\xi }_k,{\eta }_k,{\zeta }_k$ are the Gaussian random variables which follow the distribution $N(0,1),k=1,2,3.$

In system (3), according to [27]:

Example 1.

In order to check the existence of an ergodic stationary distribution, we choose the values of the system parameters as follows: $({\sigma }_1,{\sigma }_2,{\sigma }_3)=(0.1,0.1,0.1)$, $\lambda =1000,d=0.1,e=0.002,a=5,p=0.2,c=0.2,b=0.3,\sigma =0.2$, then $R_0^{*}=3.208>1,$ where $R_0^{*}$ is defined before Theorem 2 and choosing $c_1=-15,M=7$. In other words, the conditions of Theorem 2 hold, and there is an ergodic stationary distribution $\mu (·)$ of system (3), which will persist for a long time. Figure 1 confirms this.

Example 2.

In order to check the extinction of the system (3), we choose the values of the system parameters as follows: $({\sigma }_1,{\sigma }_2,{\sigma }_3)=(0.1,0.1,0.1)$, $\lambda =250,d=0.1,e=0.002,a=5,p=0.2,c=0.2,b=0.3,\sigma =0.2$; then, $R_0^{*}=0.849<1,$ where $R_0^{*}$ is defined before Theorem 3 and choosing $c_1=-20,M=5$. In other words, the conditions of Theorem 3 hold, which will be extinct in a long time. Figure 2 confirms this.

7. Discussion

This paper considers the parameters d, a, and b disturbed by the white noise and assumes the dynamics of the corresponding stochastic system (3) with the time delay and CTL responsiveness. The reason for choosing the three disturbed parameters is that these three parameters are important in controlling the viral disease. Of course, in the following research, we will focus on the general situation to investigate the influence of white noise.