Dynamics Analysis of a Nonlinear Stochastic SEIR Epidemic System with Varying Population Size

Abstract

This paper considers a stochastic susceptible exposed infectious recovered (SEIR) epidemic model with varying population size and vaccination. We aim to study the global dynamics of the reduced nonlinear stochastic proportional differential system. We first investigate the existence and uniqueness of global positive solution of the stochastic system. Then the sufficient conditions for the extinction and permanence in mean of the infectious disease are obtained. Furthermore, we prove that the solution of the stochastic system has a unique ergodic stationary distribution under appropriate conditions. Finally, the discussion and numerical simulation are given to demonstrate the obtained results.

1. Introduction

Since the pioneering work of Kermack and Mckendrick [1], mathematical modeling for the dynamics of epidemic transmission has a realistic significance in predicting and controlling the spread of infectious diseases in the field of epidemiological research [2,3,4,5,6,7,8]. Recently, stochastic differential equations have been widely applied to physics, engineering, chemistry, and biology [9,10,11,12,13,14,15,16,17,18,19,20,21], which have obtained some novel results.

In fact, with the development of modern medicine, vaccination has become an important strategy for disease control. Then numerous scholars have investigated the effect of vaccination on disease [22,23,24,25,26,27]. The epidemic model with a constant population size is relatively effective for diseases with a low mortality and short duration. However, it is clearly untenable for diseases with a high mortality and varying populations. Thus epidemic models with varying population size seem to be more reasonable, which have attracted much interest from the research scientists [28,29,30]. Moreover, many infectious diseases incubate inside the hosts for a period of time before becoming infectious, so it is very meaningful to consider the effect of the incubation period. Based on the above considerations, Sun et al. [28] studied an SEIR model with varying population size and vaccination. The system can be described by

$$\left\{\begin{array}{cc}\hfill \dot{S}=& bN-\frac{(1-\delta )\beta SI}{N}-\frac{\delta (1-p)\beta SI}{N}-\delta pS-\mu S,\hfill \\ \hfill \dot{E}=& \frac{(1-\delta )\beta SI}{N}+\frac{\delta (1-p)\beta SI}{N}-\alpha E-\mu E,\hfill \\ \hfill \dot{I}=& \alpha E-(\epsilon +\gamma +\mu )I,\hfill \\ \hfill \dot{R}=& \delta pS+\gamma I-\mu R,\hfill \end{array}\right.$$

(1)

where $S\left(t\right)$, $E\left(t\right)$, $I\left(t\right)$ and $R\left(t\right)$, respectively, stand for the densities of the susceptible, the exposed, the infective and recovered individuals at time t, the total population size is denoted by $N\left(t\right)=S\left(t\right)+E\left(t\right)+I\left(t\right)+R\left(t\right)$. b represents the inflow rate (including birth and immigration), $\mu$ denotes the outflow rate (including natural death and emigration). The function $\frac{\beta SI}{N}$ stands for the standard incidence rate, here $\beta$ represents the transmission rate of disease. $\delta$$(0\le \delta <1)$ is the vaccine coverage rate of susceptible individuals, p$(0\le p\le 1)$ is the vaccine efficacy, $\alpha$ represents the rate at which the exposed individuals become infectious, $\epsilon$ is the rate of disease-related death and $\gamma$ stands for the recovery rate of infective individuals. The parameters $\delta$ and p are all non-negative constants and b, $\mu$, $\beta$, $\alpha$, $\epsilon$ and $\gamma$ are positive constants. Moreover, the differential equation of total population size $N\left(t\right)$ is given by $\dot{N}=(b-\mu )N-\epsilon I$. The authors [28] explored the proportions of individuals in the four epidemiological classes, namely

$$\tilde{s}=\frac{S}{N},\tilde{e}=\frac{E}{N},\tilde{i}=\frac{I}{N},\tilde{r}=\frac{R}{N}.$$

(2)

It is easy to get that the variables $\tilde{s}$, $\tilde{e}$, $\tilde{i}$ and $\tilde{r}$ satisfy the following system of differential equations

$$\left\{\begin{array}{cc}\hfill \dot{\tilde{s}}=& b-(1-\delta p)\beta \tilde{s}\tilde{i}-(\delta p+b)\tilde{s}+\epsilon \tilde{s}\tilde{i},\hfill \\ \hfill \dot{\tilde{e}}=& (1-\delta p)\beta \tilde{s}\tilde{i}-(\alpha +b)\tilde{e}+\epsilon \tilde{e}\tilde{i},\hfill \\ \hfill \dot{\tilde{i}}=& \alpha \tilde{e}-(\epsilon +\gamma +b)\tilde{i}+\epsilon {\tilde{i}}^2,\hfill \\ \hfill \dot{\tilde{r}}=& \delta p\tilde{s}+\gamma \tilde{i}-b\tilde{r}+\epsilon \tilde{i}\tilde{r}.\hfill \end{array}\right.$$

Since variable $\tilde{r}$ does not appear in the first, second, third equations of the above system. Then the above system becomes the following reduced system

$$\left\{\begin{array}{cc}\hfill \dot{\tilde{s}}=& b-(1-\delta p)\beta \tilde{s}\tilde{i}-(\delta p+b)\tilde{s}+\epsilon \tilde{s}\tilde{i},\hfill \\ \hfill \dot{\tilde{e}}=& (1-\delta p)\beta \tilde{s}\tilde{i}-(\alpha +b)\tilde{e}+\epsilon \tilde{e}\tilde{i},\hfill \\ \hfill \dot{\tilde{i}}=& \alpha \tilde{e}-(\epsilon +\gamma +b)\tilde{i}+\epsilon {\tilde{i}}^2\hfill \end{array}\right.$$

(3)

which is subject to the constraint $\tilde{r}=1-\tilde{s}-\tilde{e}-\tilde{i}$. In the region $△=\left\{\left(\tilde{s},\tilde{e},\tilde{i}\right)\in {\mathbb{R}}_{+}^3|0\le \tilde{s}+\tilde{e}+\tilde{i}\le 1\right\}$, they established the epidemiological threshold condition $R_0$, which determines disease extinction or permanence, where

$$R_0=\frac{b\alpha \beta (1-\delta p)}{(\alpha +b)(\delta p+b)(\epsilon +\gamma +b)}.$$

Meanwhile, they analyzed the global dynamics of system (3) and derived the equilibria (including the disease-free equilibrium and the endemic equilibrium) and their global stability. In addition, the parameter restrictions for uniform permanence were obtained.

Nevertheless, the biological populations in the ecosystem are inevitably subjected to uncertain environmental perturbations. It is worth noting that this phenomenon is ubiquitous in the natural environment. So various stochastic epidemic models have been proposed and studied [31,32,33,34,35,36]. To the best of our knowledge, there are not too many researches on global dynamics of the stochastic SEIR epidemic model with varying population size and vaccination yet. In this paper, to make this epidemic model (1) more reasonable and realistic, we suppose the stochastic perturbations are directly proportional to $\tilde{s}$, $\tilde{e}$, $\tilde{i}$ and $\tilde{r}$ under the influence of white noise type, influenced on the $\dot{\tilde{s}}\left(t\right)$, $\dot{\tilde{e}}\left(t\right)$, $\dot{\tilde{i}}\left(t\right)$ and $\dot{\tilde{r}}\left(t\right)$ in system (1), respectively. This implies the stochastic effects of white noise on the birth and death rates of $S,E,I,R$. Then corresponding to system (1), a stochastic version can be reached by

$$\left\{\begin{array}{cc}\mathrm{d}S=\hfill & \left[bN-\frac{(1-\delta )\beta SI}{N}-\frac{\delta (1-p)\beta SI}{N}-\delta pS-\mu S\right]\mathrm{d}t+{\sigma }_{1}S\mathrm{d}{B}_1\left(t\right),\hfill \\ \mathrm{d}E=\hfill & \left[\frac{(1-\delta )\beta SI}{N}+\frac{\delta (1-p)\beta SI}{N}-\alpha E-\mu E\right]\mathrm{d}t+{\sigma }_{2}E\mathrm{d}{B}_2\left(t\right),\hfill \\ \mathrm{d}I=\hfill & \left[\alpha E-(\epsilon +\gamma +\mu )I\right]\mathrm{d}t+{\sigma }_{3}I\mathrm{d}{B}_3\left(t\right),\hfill \\ \mathrm{d}R=\hfill & \left(\delta pS+\gamma I-\mu R\right)\mathrm{d}t+{\sigma }_{4}R\mathrm{d}{B}_4\left(t\right),\hfill \end{array}\right.$$

(4)

where $B_i\left(t\right)(i=1,2,3,4)$ is the standard Wiener processes with $B_i\left(0\right)=0$ a.s. ${\sigma }_i\left(t\right)(i=1,2,3,4)$ stands for a continuous and bounded function for any $t\ge 0$ and ${\sigma }_i^2\left(t\right)(i=1,2,3,4)$ represents the intensities of Wiener processes. Furthermore, the differential equation of total population size $N\left(t\right)$ is given by the following form

$$\mathrm{d}N=[(b-\mu )N-\epsilon I]\mathrm{d}t+{\sigma }_{1}S\mathrm{d}{B}_1\left(t\right)+{\sigma }_{2}E\mathrm{d}{B}_2\left(t\right)+{\sigma }_{3}I\mathrm{d}{B}_3\left(t\right)+{\sigma }_{4}R\mathrm{d}{B}_4\left(t\right).$$

From (2), the system (4) becomes the following proportional system

$$\left\{\begin{array}{cc}\hfill \mathrm{d}\tilde{s}=& \left[b-(1-\delta p)\beta \tilde{s}\tilde{i}-(\delta p+b)\tilde{s}+\epsilon \tilde{s}\tilde{i}-{\sigma }_1^2{\tilde{s}}^2+\tilde{s}\left({\sigma }_1^2{\tilde{s}}^2+{\sigma }_2^2{\tilde{e}}^2+{\sigma }_3^2{\tilde{i}}^2+{\sigma }_4^2{\tilde{r}}^2\right)\right]\mathrm{d}t\hfill \\ & +{\sigma }_1\tilde{s}(1-\tilde{s})\mathrm{d}{B}_1\left(t\right)-{\sigma }_2\tilde{s}\tilde{e}\mathrm{d}{B}_2\left(t\right)-{\sigma }_3\tilde{s}\tilde{i}\mathrm{d}{B}_3\left(t\right)-{\sigma }_4\tilde{s}\tilde{r}\mathrm{d}{B}_4\left(t\right),\hfill \\ \hfill \mathrm{d}\tilde{e}=& \left[(1-\delta p)\beta \tilde{s}\tilde{i}-(\alpha +b)\tilde{e}+\epsilon \tilde{e}\tilde{i}-{\sigma }_2^2{\tilde{e}}^2+\tilde{e}\left({\sigma }_1^2{\tilde{s}}^2+{\sigma }_2^2{\tilde{e}}^2+{\sigma }_3^2{\tilde{i}}^2+{\sigma }_4^2{\tilde{r}}^2\right)\right]\mathrm{d}t\hfill \\ & -{\sigma }_1\tilde{s}\tilde{e}\mathrm{d}{B}_1\left(t\right)+{\sigma }_2\tilde{e}(1-\tilde{e})\mathrm{d}{B}_2\left(t\right)-{\sigma }_3\tilde{e}\tilde{i}\mathrm{d}{B}_3\left(t\right)-{\sigma }_4\tilde{e}\tilde{r}\mathrm{d}{B}_4\left(t\right),\hfill \\ \hfill \mathrm{d}\tilde{i}=& \left[\alpha \tilde{e}-(\epsilon +\gamma +b)\tilde{i}+\epsilon {\tilde{i}}^2-{\sigma }_3^2{\tilde{i}}^2+\tilde{i}\left({\sigma }_1^2{\tilde{s}}^2+{\sigma }_2^2{\tilde{e}}^2+{\sigma }_3^2{\tilde{i}}^2+{\sigma }_4^2{\tilde{r}}^2\right)\right]\mathrm{d}t\hfill \\ & -{\sigma }_1\tilde{s}\tilde{i}\mathrm{d}{B}_1\left(t\right)-{\sigma }_2\tilde{e}\tilde{i}\mathrm{d}{B}_2\left(t\right)+{\sigma }_3\tilde{i}\left(1-\tilde{i}\right)\mathrm{d}{B}_3\left(t\right)-{\sigma }_4\tilde{i}\tilde{r}\mathrm{d}{B}_4\left(t\right),\hfill \\ \hfill \mathrm{d}\tilde{r}=& \left[\delta p\tilde{s}+\gamma \tilde{i}-b\tilde{r}+\epsilon \tilde{i}\tilde{r}-{\sigma }_4^2{\tilde{r}}^2+\tilde{r}\left({\sigma }_1^2{\tilde{s}}^2+{\sigma }_2^2{\tilde{e}}^2+{\sigma }_3^2{\tilde{i}}^2+{\sigma }_4^2{\tilde{r}}^2\right)\right]\mathrm{d}t\hfill \\ & -{\sigma }_1\tilde{s}\tilde{r}\mathrm{d}{B}_1\left(t\right)-{\sigma }_2\tilde{e}\tilde{r}\mathrm{d}{B}_2\left(t\right)-{\sigma }_3\tilde{i}\tilde{r}\mathrm{d}{B}_3\left(t\right)+{\sigma }_4\tilde{r}(1-\tilde{r})\mathrm{d}{B}_4\left(t\right).\hfill \end{array}\right.$$

(5)

It is worthy to note that, the variables $\tilde{s}$, $\tilde{e}$, $\tilde{i}$ and $\tilde{r}$ satisfy the relation $\tilde{r}=1-\tilde{s}-\tilde{e}-\tilde{i}$, we can omit analysis of the fourth equation of system (5) and explore the following reduced system

$$\left\{\begin{array}{cc}\mathrm{d}\tilde{s}=\hfill & [b-(1-\delta p)\beta \tilde{s}\tilde{i}-(\delta p+b)\tilde{s}+\epsilon \tilde{s}\tilde{i}-{\sigma }_1^2{\tilde{s}}^2+\tilde{s}({\sigma }_1^2{\tilde{s}}^2+{\sigma }_2^2{\tilde{e}}^2+{\sigma }_3^2{\tilde{i}}^2\hfill \\ & +{\sigma }_4^2{\left(1-\tilde{s}-\tilde{e}-\tilde{i}\right)}^2\left)\right]\mathrm{d}t+{\sigma }_1\tilde{s}(1-\tilde{s})\mathrm{d}{B}_1\left(t\right)-{\sigma }_2\tilde{s}\tilde{e}\mathrm{d}{B}_2\left(t\right)-{\sigma }_3\tilde{s}\tilde{i}\mathrm{d}{B}_3\left(t\right)\hfill \\ & -{\sigma }_4\tilde{s}\left(1-\tilde{s}-\tilde{e}-\tilde{i}\right)\mathrm{d}{B}_4\left(t\right),\hfill \\ \mathrm{d}\tilde{e}=\hfill & [(1-\delta p)\beta \tilde{s}\tilde{i}-(\alpha +b)\tilde{e}+\epsilon \tilde{e}\tilde{i}-{\sigma }_2^2{\tilde{e}}^2+\tilde{e}({\sigma }_1^2{\tilde{s}}^2+{\sigma }_2^2{\tilde{e}}^2+{\sigma }_3^2{\tilde{i}}^2\hfill \\ & +{\sigma }_4^2{\left(1-\tilde{s}-\tilde{e}-\tilde{i}\right)}^2\left)\right]\mathrm{d}t-{\sigma }_1\tilde{s}\tilde{e}\mathrm{d}{B}_1\left(t\right)+{\sigma }_2\tilde{e}(1-\tilde{e})\mathrm{d}{B}_2\left(t\right)-{\sigma }_3\tilde{e}\tilde{i}\mathrm{d}{B}_3\left(t\right)\hfill \\ & -{\sigma }_4\tilde{e}\left(1-\tilde{s}-\tilde{e}-\tilde{i}\right)\mathrm{d}{B}_4\left(t\right),\hfill \\ \mathrm{d}\tilde{i}=\hfill & \left[\alpha \tilde{e}-(\epsilon +\gamma +b)\tilde{i}+\epsilon {\tilde{i}}^2-{\sigma }_3^2{\tilde{i}}^2+\tilde{i}\left({\sigma }_1^2{\tilde{s}}^2+{\sigma }_2^2{\tilde{e}}^2+{\sigma }_3^2{\tilde{i}}^2+{\sigma }_4^2{\left(1-\tilde{s}-\tilde{e}-\tilde{i}\right)}^2\right)\right]\mathrm{d}t\hfill \\ & -{\sigma }_1\tilde{s}\tilde{i}\mathrm{d}{B}_1\left(t\right)-{\sigma }_2\tilde{e}\tilde{i}\mathrm{d}{B}_2\left(t\right)+{\sigma }_3\tilde{i}\left(1-\tilde{i}\right)\mathrm{d}{B}_3\left(t\right)-{\sigma }_4\tilde{i}\left(1-\tilde{s}-\tilde{e}-\tilde{i}\right)\mathrm{d}{B}_4\left(t\right)\hfill \end{array}\right.$$

(6)

with the initial value $\left(\tilde{s}\left(0\right),\tilde{e}\left(0\right),\tilde{i}\left(0\right)\right)\in {\mathbb{R}}_{+}^3$ and $\tilde{s}\left(0\right)+\tilde{e}\left(0\right)+\tilde{i}\left(0\right)<1$.

Since system (6) is a three-dimensional stochastic system with many high-order nonlinear terms, this makes the stochastic analysis novel and more complex than [34,36].

Throughout this article, unless otherwise specified, let $(\mathsf{\Omega },\mathcal{F},{\left\{\mathcal{F}\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 increasing and right continuous while ${\mathcal{F}}_0$ contains all $\mathbb{P}$-null sets). Further suppose $B_i\left(t\right)(i=1,2,3,4)$ stands for the mutually independent standard Wiener processes defined on the complete probability space $\mathsf{\Omega }$. For an integrable function $x\left(t\right)$ on $[0,+\infty )$, let us define $\langle x\left(t\right)\rangle =\frac{1}{t}{\int }_0^{t}x\left(r\right)\mathrm{d}r$.

2. Global Positive Solution

The following Itô’s formula will be used frequently in the sequel.

Lemma 1.

[37] Assume that $X\left(t\right)\in R^{+}$ is an Itô’s process of the form

$$dX\left(t\right)=F(X\left(t^{-}\right),t^{-})dt+G(X\left(t^{-}\right),t^{-})dB\left(t\right),$$

where $F:R^n\times R_{+}\times S\to R^n$ and $G:R^n\times R_{+}\times S\to R^n$ are measurable functions.

Given $V\in C^{2,1}(R^n\times R_{+}\times S;R_{+})$, we define the operator $LV$ by

$$\begin{array}{cc}\hfill LV(X,t)=& V_t(X,t)+V_X(X,t)F(X,t)+\frac{1}{2}\mathit{trace}\left[G^T(X,t)V_{XX}(X,t)G(X,t)\right],\hfill \end{array}$$

where

$$V_t(X,t)=\frac{\partial V_X(X,t)}{\partial t},V_X(X,t)=\left(\frac{\partial V_X(X,t)}{\partial X_1},\cdots ,\frac{\partial V_X(X,t)}{\partial X_n}\right),V_{XX}(X,t)={\left(\frac{{\partial }^2{V}_X(X,t)}{\partial X_i\partial X_j}\right)}_{n\times n}.$$

Then the generalized Itô’s formula is given by

$$\begin{array}{cc}\hfill dV(X,t)=& LV(X,t)dt+V_X(X,t)G(X,t)dB\left(t\right).\hfill \end{array}$$

To explore the dynamical behaviors of a population system, we first concern the global existence and positivity of the solutions of system (6).

Lemma 2.

For any given initial value $\left(\tilde{s}\left(0\right),\tilde{e}\left(0\right),\tilde{i}\left(0\right)\right)\in {\mathbb{R}}_{+}^3$ and $\tilde{s}\left(0\right)+\tilde{e}\left(0\right)+\tilde{i}\left(0\right)<1$, the system (6) has a unique positive local solution $(\tilde{s}\left(t\right),\tilde{e}\left(t\right),\tilde{i}\left(t\right))$ for $t\in [-\omega ,{\tau }_e)$, where ${\tau }_e$ is the explosion time [37].

Theorem 1.

For any given initial value $\left(\tilde{s}\left(0\right),\tilde{e}\left(0\right),\tilde{i}\left(0\right)\right)\in {\mathbb{R}}_{+}^3$ and $\tilde{s}\left(0\right)+\tilde{e}\left(0\right)+\tilde{i}\left(0\right)<1$, the system (6) has a unique positive solution $\left(\tilde{s}\left(t\right),\tilde{e}\left(t\right),\tilde{i}\left(t\right)\right)\in {\mathbb{R}}_{+}^3$ on $t>0$ a.s.

Proof. 

The following proof is divided into two parts.

Part I. Since the coefficients of the system (6) satisfy local Lipschitz condition, from Lemma 2, it is easy to see that the system (6) has a unique positive local solution $\left(\tilde{s}\left(t\right),\tilde{e}\left(t\right),\tilde{i}\left(t\right)\right)$ for any given initial value $\left(\tilde{s}\left(0\right),\tilde{e}\left(0\right),\tilde{i}\left(0\right)\right)\in {\mathbb{R}}_{+}^3$ and $\tilde{s}\left(0\right)+\tilde{e}\left(0\right)+\tilde{i}\left(0\right)<1$.

Part II. Now we prove that the positive solution is global, that is ${\tau }_e=\infty$ a.s. Let $k_0\ge 0$ be sufficiently large such that $\tilde{s}\left(0\right)$, $\tilde{e}\left(0\right)$ and $\tilde{i}\left(0\right)$ all lie in $\left[\frac{1}{k_0},k_0\right]$. For each integer $k\ge k_0$, let us define the stopping time

$${\tau }_k=\inf \left\{t\in [-\omega ,{\tau }_e):\tilde{s}\left(t\right)\notin \left(\frac{1}{k},k\right),\tilde{e}\left(t\right)\notin \left(\frac{1}{k},k\right)\mathrm{or}\tilde{i}\left(t\right)\notin \left(\frac{1}{k},k\right)\right\},$$

where we define $\inf \varnothing =\infty$ (∅ stands for the empty set). Evidently, ${\tau }_k$ is strictly increasing when $k\to \infty$. Let ${\displaystyle {\tau }_{\infty }=\underset{k\to \infty }{\lim }{\tau }_k}$, thus ${\tau }_{\infty }\le {\tau }_e$ a.s. So we just need to show that ${\tau }_{\infty }=\infty$ a.s. If ${\tau }_{\infty }=\infty$ is untrue, then there exist two constants $T>0$ and $ϵ\in (0,1)$ such that $\mathbb{P}\{{\tau }_{\infty }\le T\}>ϵ$. Hence, there exists $k_1\ge k_0(k_1\in N_{+})$ such that

$$\begin{array}{c}\hfill \mathbb{P}\{{\tau }_k\le T\}\ge ϵ,k\ge k_1.\end{array}$$

(7)

Define a $C^2$-function V: ${\mathbb{R}}_{+}^3\to {\mathbb{R}}_{+}$ by

$$\begin{array}{cc}\hfill V\left(\tilde{s},\tilde{e},\tilde{i}\right)=& -\ln \left(1-\tilde{s}-\tilde{e}-\tilde{i}\right)-\ln \tilde{s}-\ln \tilde{e}-\ln \tilde{i}-3.\hfill \end{array}$$

The non-negativity of $V(\tilde{s},\tilde{e},\tilde{i})$ can be obtained from $m-1-\ln m\ge 0,m>0$.

In terms of the multi-dimensional Itô’s formula and system (6), we have

$$\begin{array}{cc}\mathrm{d}V=\hfill & LV\mathrm{d}t+{\sigma }_1\left(4\tilde{s}-1\right)\mathrm{d}{B}_1\left(t\right)+{\sigma }_2\left(4\tilde{e}-1\right)\mathrm{d}{B}_2\left(t\right)+{\sigma }_3\left(4\tilde{i}-1\right)\mathrm{d}{B}_3\left(t\right)\hfill \\ & +{\sigma }_4\left(3-4\tilde{s}-4\tilde{e}-4\tilde{i}\right)\mathrm{d}{B}_4\left(t\right),\hfill \end{array}$$

where $LV$ is given in Appendix A in detail. Then we have

$$\begin{array}{cc}\hfill LV\le & \beta +\delta p+4b+\epsilon +\alpha +\gamma +\frac{1}{2}{\sigma }_1^2+\frac{1}{2}{\sigma }_2^2+\frac{1}{2}{\sigma }_3^2+\frac{1}{2}{\sigma }_4^2\hfill \\ \hfill :=& M_0,\hfill \end{array}$$

where $M_0$ is a positive constant.

So we get

$$\begin{array}{cc}\mathrm{d}V\le \hfill & M_0\mathrm{d}t+{\sigma }_1\left(4\tilde{s}-1\right)\mathrm{d}{B}_1\left(t\right)+{\sigma }_2\left(4\tilde{e}-1\right)\mathrm{d}{B}_2\left(t\right)+{\sigma }_3\left(4\tilde{i}-1\right)\mathrm{d}{B}_3\left(t\right)\hfill \\ & +{\sigma }_4\left(3-4\tilde{s}-4\tilde{e}-4\tilde{i}\right)\mathrm{d}{B}_4\left(t\right).\hfill \end{array}$$

(8)

Integrating both sides of (8) from 0 to ${\tau }_k\wedge T$ and then taking the expectation yield

$$\begin{array}{cc}\hfill \mathbb{E}V\left(\tilde{s}({\tau }_k\wedge T),\tilde{e}({\tau }_k\wedge T),\tilde{i}({\tau }_k\wedge T)\right)\le & V\left(\tilde{s}\left(0\right),\tilde{e}\left(0\right),\tilde{i}\left(0\right)\right)+\mathbb{E}{\int }_0^{{\tau }_k\wedge T}{M}_0\mathrm{d}t\hfill \\ \hfill \le & V\left(\tilde{s}\left(0\right),\tilde{e}\left(0\right),\tilde{i}\left(0\right)\right)+M_{0}T.\hfill \end{array}$$

(9)

Let ${\mathsf{\Omega }}_k=\{{\tau }_k\le T\},k\ge k_1$ and from (7), we have $P\left({\mathsf{\Omega }}_k\right)\ge ϵ$. Notice that for every $\omega \in {\mathsf{\Omega }}_k$, there exists $\tilde{s}({\tau }_k,\omega ),\tilde{e}({\tau }_k,\omega )$ or $\tilde{i}({\tau }_k,\omega )$ equals either $\frac{1}{k}$ or k. Thus

$$\begin{array}{c}\hfill V\left(\tilde{s}({\tau }_k,\omega ),\tilde{e}({\tau }_k,\omega ),\tilde{i}({\tau }_k,\omega )\right)\ge \left(\frac{1}{k}-1-\ln \frac{1}{k}\right)\wedge (k-1-\ln k).\end{array}$$

(10)

By virtue of (9) and (10), we have

$$\begin{array}{cc}\hfill V\left(\tilde{s}\left(0\right),\tilde{e}\left(0\right),\tilde{i}\left(0\right)\right)+M_{0}T\ge & \mathbb{E}\left[1_{{\mathsf{\Omega }}_k\left(\omega \right)}V\left(\tilde{s}({\tau }_k,\omega ),\tilde{e}({\tau }_k,\omega ),\tilde{i}({\tau }_k,\omega )\right)\right]\hfill \\ \hfill \ge & ϵ\left[\left(\frac{1}{k}-1-\ln \frac{1}{k}\right)\wedge (k-1-\ln k)\right],\hfill \end{array}$$

here $1_{{\mathsf{\Omega }}_k\left(\omega \right)}$ represents the indicator function of ${\mathsf{\Omega }}_k\left(\omega \right)$.

Let $k\to \infty$, which implies

$$\infty >V\left(\tilde{s}\left(0\right),\tilde{e}\left(0\right),\tilde{i}\left(0\right)\right)+M_{0}T=\infty$$

is a contradiction. Obviously, we get that ${\tau }_{\infty }=\infty$. The proof of Theorem 1 is complete.  ☐

3. Extinction

For a population system, the parameter conditions of disease extinction and permanence have become an important issue that attracts more and more attention in real life. In this section, we mainly investigate the extinction of disease and leave the argument of permanence to the next section.

Theorem 2.

Let $\left(\tilde{s}\left(t\right),\tilde{e}\left(t\right),\tilde{i}\left(t\right)\right)$ be the solution of system (6) with the initial value $\left(\tilde{s}\left(0\right),\tilde{e}\left(0\right),\tilde{i}\left(0\right)\right)$$\in {\mathbb{R}}_{+}^3$ and $\tilde{s}\left(0\right)+\tilde{e}\left(0\right)+\tilde{i}\left(0\right)<1$. If the parameter conditions

$$M_1<2(1-\varrho )\alpha ,\epsilon \le (1-\delta p)\beta <\epsilon +\gamma +b$$

hold, then

$$\underset{t\to \infty }{\lim \sup }\frac{\ln \left(\tilde{e}\left(t\right)+\varrho \tilde{i}\left(t\right)\right)}{t}\le (\varrho -1)\alpha +\frac{M_1}{2}<0a.s.,$$

where

$$\varrho =\frac{-(\epsilon +\gamma +b-\alpha )+\sqrt{{(\epsilon +\gamma +b-\alpha )}^2+4(1-\delta p)\alpha \beta }}{2\alpha }$$

and

$$M_1=\max \left\{{\sigma }_1^2,\left(1-{\varrho }^2\right){\sigma }_3^2,{\sigma }_4^2\right\},$$

namely, $\tilde{e}\left(t\right)$ and $\tilde{i}\left(t\right)$ tend to zero exponentially a.s. That is to say, the exposed and infective individuals go to extinction almost surely.

Proof. 

Let us define a differentiable function V by

$$V=\ln \left(\tilde{e}\left(t\right)+\varrho \tilde{i}\left(t\right)\right),$$

here $\varrho$ is a positive constant to be determined later. According to the Itô’s formula and system (6), we have

$$\begin{array}{cc}\mathrm{d}V=\hfill & LV\mathrm{d}t-{\sigma }_1\tilde{s}\mathrm{d}{B}_1\left(t\right)+\frac{{\sigma }_2\tilde{e}\left(1-\tilde{e}-\varrho \tilde{i}\right)}{\tilde{e}+\varrho \tilde{i}}\mathrm{d}{B}_2\left(t\right)+\frac{{\sigma }_3\tilde{i}\left(\varrho -\tilde{e}-\varrho \tilde{i}\right)}{\tilde{e}+\varrho \tilde{i}}\mathrm{d}{B}_3\left(t\right)\hfill \\ & -{\sigma }_4\left(1-\tilde{s}-\tilde{e}-\tilde{i}\right)\mathrm{d}{B}_4\left(t\right),\hfill \end{array}$$

(11)

where $LV$ is given in Appendix B in detail. One can derive that

$$\begin{array}{cc}\mathrm{d}V\le \hfill & \left[(\varrho -1)\alpha +\frac{M_1}{2}\right]\mathrm{d}t-{\sigma }_1\tilde{s}\mathrm{d}{B}_1\left(t\right)+\frac{{\sigma }_2\tilde{e}\left(1-\tilde{e}-\varrho \tilde{i}\right)}{\tilde{e}+\varrho \tilde{i}}\mathrm{d}{B}_2\left(t\right)+\frac{{\sigma }_3\tilde{i}\left(\varrho -\tilde{e}-\varrho \tilde{i}\right)}{\tilde{e}+\varrho \tilde{i}}\mathrm{d}{B}_3\left(t\right)\hfill \\ & -{\sigma }_4\left(1-\tilde{s}-\tilde{e}-\tilde{i}\right)\mathrm{d}{B}_4\left(t\right),\hfill \end{array}$$

(12)

here $M_1=\max \left\{{\sigma }_1^2,\left(1-{\varrho }^2\right){\sigma }_3^2,{\sigma }_4^2\right\}$. Then, integrating from 0 to t and dividing by t on both sides of (12) yield

$$\begin{array}{c}\hfill \frac{\ln \left(\tilde{e}\left(t\right)+\varrho \tilde{i}\left(t\right)\right)}{t}\le (\varrho -1)\alpha +\frac{M_1}{2}+\frac{\ln \left(\tilde{e}\left(0\right)+\varrho \tilde{i}\left(0\right)\right)}{t}+\frac{\tilde{M}\left(t\right)}{t},\end{array}$$

here

$$\begin{array}{cc}\tilde{M}\left(t\right)=\hfill & -{\sigma }_1{\int }_0^t\tilde{s}\left(r\right)\mathrm{d}{B}_1\left(r\right)+{\sigma }_2{\int }_0^t\frac{\tilde{e}\left(r\right)\left(1-\tilde{e}\left(r\right)-\varrho \tilde{i}\left(r\right)\right)}{\tilde{e}\left(r\right)+\varrho \tilde{i}\left(r\right)}\mathrm{d}{B}_2\left(r\right)\hfill \\ & +{\sigma }_3{\int }_0^t\frac{\tilde{i}\left(r\right)\left(\varrho -\tilde{e}\left(r\right)-\varrho \tilde{i}\left(r\right)\right)}{\tilde{e}\left(r\right)+\varrho \tilde{i}\left(r\right)}\mathrm{d}{B}_3\left(r\right)-{\sigma }_4{\int }_0^t\left(1-\tilde{s}\left(r\right)-\tilde{e}\left(r\right)-\tilde{i}\left(r\right)\right)\mathrm{d}{B}_4\left(r\right).\hfill \end{array}$$

In a similar way as [38], making use of the strong law of large numbers [37] yields

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

Therefore,

$$\underset{t\to \infty }{\lim \sup }\frac{\ln \left(\tilde{e}\left(t\right)+\varrho \tilde{i}\left(t\right)\right)}{t}\le (\varrho -1)\alpha +\frac{M_1}{2}<0a.s.,$$

which shows that

$$\underset{t\to \infty }{\lim }\tilde{e}\left(t\right)=0,\underset{t\to \infty }{\lim }\tilde{i}\left(t\right)=0a.s.$$

The proof of Theorem 2 is complete. ☐

4. Permanence in Mean

Theorem 3.

Let $\left(\tilde{s}\left(t\right),\tilde{e}\left(t\right),\tilde{i}\left(t\right)\right)$ be the solution of system (6) with the initial value $\left(\tilde{s}\left(0\right),\tilde{e}\left(0\right),\tilde{i}\left(0\right)\right)$$\in {\mathbb{R}}_{+}^3$ and $\tilde{s}\left(0\right)+\tilde{e}\left(0\right)+\tilde{i}\left(0\right)<1$. If the parameter condition

$$\sqrt[3]{b(1-\delta p)\beta \alpha }>\frac{\delta p+3b+\alpha +\epsilon +\gamma +\frac{1}{2}{\sigma }_1^2+\frac{1}{2}{\sigma }_2^2+\frac{1}{2}{\sigma }_3^2}{3}$$

holds, then

$$\underline{\tilde{i}}\le \underset{t\to \infty }{\lim \inf }〈\tilde{i}\left(t\right)〉\le \underset{t\to \infty }{\lim \sup }〈\tilde{i}\left(t\right)〉\le \overline{\tilde{i}}a.s.,$$

where

$$\underline{\tilde{i}}=\frac{3\sqrt[3]{b(1-\delta p)\beta \alpha }-\left(\delta p+3b+\alpha +\epsilon +\gamma +\frac{1}{2}{\sigma }_1^2+\frac{1}{2}{\sigma }_2^2+\frac{1}{2}{\sigma }_3^2\right)}{(1-\delta p)\beta }$$

and

$$\overline{\tilde{i}}=\frac{2b+\alpha +\epsilon +\gamma +\frac{1}{2}{\sigma }_2^2+\frac{1}{2}{\sigma }_3^2}{2\epsilon },$$

that is to say, the infective individuals $\tilde{i}\left(t\right)$ are permanent in mean almost surely.

Proof. 

The following proof is divided into two steps.

Step I. According to the Itô’s formula and system (6), we have

$$\begin{array}{cc}\mathrm{d}\left(\ln \tilde{s}+\ln \tilde{e}+\ln \tilde{i}\right)=\hfill & [\frac{b}{\tilde{s}}-(1-\delta p)\beta \tilde{i}+3\epsilon \tilde{i}+\frac{(1-\delta p)\beta \tilde{s}\tilde{i}}{\tilde{e}}+\frac{\alpha \tilde{e}}{\tilde{i}}+\frac{3}{2}{\sigma }_1^2{\tilde{s}}^2+\frac{3}{2}{\sigma }_2^2{\tilde{e}}^2+\frac{3}{2}{\sigma }_3^2{\tilde{i}}^2\hfill \\ & +\frac{3}{2}{\sigma }_4^2{\left(1-\tilde{s}-\tilde{e}-\tilde{i}\right)}^2-(\delta p+3b+\alpha +\epsilon +\gamma +\frac{1}{2}{\sigma }_1^2+\frac{1}{2}{\sigma }_2^2\hfill \\ & +\frac{1}{2}{\sigma }_3^2\left)\right]\mathrm{d}t+{\sigma }_1\left(1-3\tilde{s}\right)\mathrm{d}{B}_1\left(t\right)+{\sigma }_2\left(1-3\tilde{e}\right)\mathrm{d}{B}_2\left(t\right)\hfill \\ & +{\sigma }_3\left(1-3\tilde{i}\right)\mathrm{d}{B}_3\left(t\right)-3{\sigma }_4\left(1-\tilde{s}-\tilde{e}-\tilde{i}\right)\mathrm{d}{B}_4\left(t\right).\hfill \end{array}$$

(13)

Integrating from 0 to t and dividing by t on both sides of (13) lead to

$$\begin{array}{cc}\frac{\ln \tilde{s}\left(t\right)}{t}+\frac{\ln \tilde{e}\left(t\right)}{t}+\frac{\ln \tilde{i}\left(t\right)}{t}\ge \hfill & 3\sqrt[3]{b(1-\delta p)\beta \alpha }-(1-\delta p)\beta 〈\tilde{i}〉-(\delta p+3b+\alpha +\epsilon +\gamma \hfill \\ & +\frac{1}{2}{\sigma }_1^2+\frac{1}{2}{\sigma }_2^2+\frac{1}{2}{\sigma }_3^2)+\frac{\ln \tilde{s}\left(0\right)}{t}+\frac{\ln \tilde{e}\left(0\right)}{t}+\frac{\ln \tilde{i}\left(0\right)}{t}+\frac{\underline{M}\left(t\right)}{t},\hfill \end{array}$$

here

$$\begin{array}{cc}\underline{M}\left(t\right)=\hfill & {\sigma }_1{\int }_0^t\left(1-3\tilde{s}\left(r\right)\right)\mathrm{d}{B}_1\left(r\right)+{\sigma }_2{\int }_0^t\left(1-3\tilde{e}\left(r\right)\right)\mathrm{d}{B}_2\left(r\right)+{\sigma }_3{\int }_0^t\left(1-3\tilde{i}\left(r\right)\right)\mathrm{d}{B}_3\left(r\right)\hfill \\ & -3{\sigma }_4{\int }_0^t\left(1-\tilde{s}\left(r\right)-\tilde{e}\left(r\right)-\tilde{i}\left(r\right)\right)\mathrm{d}{B}_4\left(r\right).\hfill \end{array}$$

The detail derivation process for the above inequality is given in Appendix C.

In a similar way as [38], making use of the strong law of large numbers [37] leads to

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

Then, by virtue of $-\infty <\ln \tilde{s}\left(t\right)<0$, $-\infty <\ln \tilde{e}\left(t\right)<0$, $-\infty <\ln \tilde{i}\left(t\right)<0$$(\tilde{s}+\tilde{e}+\tilde{i}+\tilde{r}=1)$ and $\delta p<1$, it is easy to get that

$$\begin{array}{c}\hfill \underset{t\to \infty }{\lim \inf }〈\tilde{i}\left(t\right)〉\ge \underline{\tilde{i}}>0a.s.\end{array}$$

Step II. Similarly, using the Itô’s formula and system (6), we have

$$\begin{array}{cc}\mathrm{d}\left(\ln \tilde{e}+\ln \tilde{i}\right)=\hfill & [\frac{(1-\delta p)\beta \tilde{s}\tilde{i}}{\tilde{e}}+\frac{\alpha \tilde{e}}{\tilde{i}}+2\epsilon \tilde{i}+{\sigma }_1^2{\tilde{s}}^2+{\sigma }_2^2{\tilde{e}}^2+{\sigma }_3^2{\tilde{i}}^2+{\sigma }_4^2{\left(1-\tilde{s}-\tilde{e}-\tilde{i}\right)}^2\hfill \\ & -\left(2b+\alpha +\epsilon +\gamma +\frac{1}{2}{\sigma }_2^2+\frac{1}{2}{\sigma }_3^2\right)]\mathrm{d}t-2{\sigma }_1\tilde{s}\mathrm{d}{B}_1\left(t\right)+{\sigma }_2\left(1-2\tilde{e}\right)\mathrm{d}{B}_2\left(t\right)\hfill \\ & +{\sigma }_3\left(1-2\tilde{i}\right)\mathrm{d}{B}_3\left(t\right)-2{\sigma }_4\left(1-\tilde{s}-\tilde{e}-\tilde{i}\right)\mathrm{d}{B}_4\left(t\right).\hfill \end{array}$$

(14)

Integrating from 0 to t and dividing by t on both sides of (14) result in

$$\begin{array}{cc}\frac{\ln \tilde{e}\left(t\right)}{t}+\frac{\ln \tilde{i}\left(t\right)}{t}=\hfill & (1-\delta p)\beta 〈\frac{\tilde{s}\tilde{i}}{\tilde{e}}〉+\alpha 〈\frac{\tilde{e}}{\tilde{i}}〉+2\epsilon 〈\tilde{i}〉+〈{\sigma }_1^2{\tilde{s}}^2+{\sigma }_2^2{\tilde{e}}^2+{\sigma }_3^2{\tilde{i}}^2+{\sigma }_4^2{\left(1-\tilde{s}-\tilde{e}-\tilde{i}\right)}^2〉\hfill \\ & -\left(2b+\alpha +\epsilon +\gamma +\frac{1}{2}{\sigma }_2^2+\frac{1}{2}{\sigma }_3^2\right)+\frac{\ln \tilde{e}\left(0\right)}{t}+\frac{\ln \tilde{i}\left(0\right)}{t}+\frac{\widehat{M}\left(t\right)}{t}\hfill \\ \hfill \ge & 2\epsilon 〈\tilde{i}〉-\left(2b+\alpha +\epsilon +\gamma +\frac{1}{2}{\sigma }_2^2+\frac{1}{2}{\sigma }_3^2\right)+\frac{\ln \tilde{e}\left(0\right)}{t}+\frac{\ln \tilde{i}\left(0\right)}{t}+\frac{\widehat{M}\left(t\right)}{t},\hfill \end{array}$$

here

$$\begin{array}{cc}\widehat{M}\left(t\right)=\hfill & -2{\sigma }_1{\int }_0^t\tilde{s}\left(r\right)\mathrm{d}{B}_1\left(r\right)+{\sigma }_2{\int }_0^t\left(1-2\tilde{e}\left(r\right)\right)\mathrm{d}{B}_2\left(r\right)+{\sigma }_3{\int }_0^t\left(1-2\tilde{i}\left(r\right)\right)\mathrm{d}{B}_3\left(r\right)\hfill \\ & -2{\sigma }_4{\int }_0^t\left(1-\tilde{s}\left(r\right)-\tilde{e}\left(r\right)-\tilde{i}\left(r\right)\right)\mathrm{d}{B}_4\left(r\right).\hfill \end{array}$$

In a similar way as [38], using the strong law of large numbers [37], we have

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

Therefore,

$$\begin{array}{c}\hfill \underset{t\to \infty }{\lim \sup }〈\tilde{i}\left(t\right)〉\le \overline{\tilde{i}}>0a.s.\end{array}$$

The proof of Theorem 3 is complete. ☐

5. Stationary Distribution and Ergodicity

Recently, the stationary distribution attract deep research interests of many authors [32,33,34,35]. The ergodicity is one of the most important properties for the stochastic system, and geometric ergodicity for finite-dimensional systems has been shown in detail and well-developed in many earlier works [39,40]. In this section, based on the theory of Khasminskii [41] and the Lyapunov function method, we explore the conditions of the existence of an ergodic stationary distribution, which shows that the epidemic disease will prevail.

Assume $X\left(t\right)$ be a time-homogeneous Markov process in ${\mathbb{D}}_n\subset {\mathbb{R}}^n$, which is described by the stochastic differential equation

$$\mathrm{d}X\left(t\right)=b\left(X\right)\mathrm{d}t+\sum _{\eta =1}^n{\sigma }_{\eta }\left(X\right)\mathrm{d}{B}_{\eta }\left(t\right),$$

here ${\mathbb{D}}_n$ stands for a n-dimensional Euclidean space.

The diffusion matrix is as follows:

$$A\left(x\right)=\left(a_{ij}\left(x\right)\right),a_{ij}\left(x\right)=\sum _{\eta =1}^n{\sigma }_{\eta }^i\left(x\right){\sigma }_{\eta }^j\left(x\right).$$

Assumption 1.

Assume that there exists a bounded domain $U\subset {\mathbb{D}}_n$ with regular boundary Γ such that $\overline{U}\subset {\mathbb{D}}_n($$\overline{U}$ is the closure of $U)$, satisfying the following properties:

(i) In the domain U and some neighborhood thereof, the smallest eigenvalue of the diffusion matrix $A\left(x\right)$ is bounded away from zero.

(ii) If $x\in {\mathbb{D}}_n\backslash U$, the mean time τ at which a path issuing from x reaches the set U is finite, and ${\displaystyle \underset{x\in \mathsf{\Theta }}{\sup }{\mathbb{E}}_x\tau <\infty }$ for every compact subset $\mathsf{\Theta }\subset {\mathbb{D}}_n$.

Lemma 3.

[41] When Assumption 1 holds, then the Markov process $X\left(t\right)$ has a stationary distribution $\pi (·)$. In addition, when $f(·)$ is a function integrable with respect to the measure π, then

$$\begin{array}{c}\hfill {\mathbb{P}}_x\left\{\underset{T\to \infty }{\lim }\frac{1}{T}{\int }_0^{T}f\left(X\left(t\right)\right)\mathrm{d}t={\int }_{{\mathbb{D}}_n}f\left(x\right)\pi \left(\mathrm{d}x\right)\right\}=1\end{array}$$

for all $x\in {\mathbb{D}}_n$.

Remark 1.

To prove Assumption 1$\left(i\right)$ [42], it suffices to demonstrate that $\mathrm{F}$ is uniformly elliptical in any bounded domain H, here

$$\mathrm{F}u=b\left(x\right)u_x+\frac{1}{2}\mathbf{trace}\left(A\left(x\right)u_{xx}\right),$$

namely, there exists a positive number Z such that

$$\sum _{i,j=1}^n{a}_{ij}\left(x\right){\xi }_i{\xi }_j\ge Z{\left|\xi \right|}^2,x\in \overline{H},\xi \in {\mathbb{R}}^n.$$

To prove Assumption 1$\left(ii\right)$ [43], it suffices to demonstrate that there exist some neighborhood U and a nonnegative $C^2$-function V such that $\forall x\in {\mathbb{D}}_n\backslash U$, $LV\left(x\right)<0$.

Making use of the Lemma 3, we can obtain the main results as follows.

Theorem 4.

Let $\left(\tilde{s}\left(t\right),\tilde{e}\left(t\right),\tilde{i}\left(t\right)\right)$ be the solution of system (6) with the initial value $\left(\tilde{s}\left(0\right),\tilde{e}\left(0\right),\tilde{i}\left(0\right)\right)$$\in {\mathbb{R}}_{+}^3$ and $\tilde{s}\left(0\right)+\tilde{e}\left(0\right)+\tilde{i}\left(0\right)<1$. If the parameter condition

$$\sqrt[3]{b(1-\delta p)\beta \alpha }>\frac{2\delta p+5b+\alpha +2\epsilon +2\gamma +2{\sigma }_1^2+2{\sigma }_2^2+2{\sigma }_3^2+\frac{1}{2}{\sigma }_4^2}{3}$$

holds, then the system (6) has a unique stationary distribution $\pi (·)$ and it has ergodic property.

Proof. 

Now let us define a positive-definite function V by

$$V=-\ln \left(\tilde{s}+\tilde{e}+\tilde{i}\right)-\ln \tilde{s}-\ln \tilde{e}-\ln \tilde{i}-\ln \tilde{r}.$$

Using the Itô’s formula yields

$$LV\le -\frac{b}{\tilde{s}}-\frac{(1-\delta p)\beta \tilde{s}\tilde{i}}{\tilde{e}}-\frac{\alpha \tilde{e}}{\tilde{i}}-\frac{b}{\tilde{s}+\tilde{e}+\tilde{i}}-\frac{\delta p\tilde{s}}{\tilde{r}}-\frac{\gamma \tilde{i}}{\tilde{r}}+M_1,$$

(15)

here

$$M_1=(1-\delta p)\beta +2\delta p+5b+\alpha +2\epsilon +2\gamma +2{\sigma }_1^2+2{\sigma }_2^2+2{\sigma }_3^2+\frac{1}{2}{\sigma }_4^2.$$

The detail derivation process for the above inequality of $LV$ is given in Appendix D.

Next let us construct the following compact subset U:

$$U=\left\{\left(\tilde{s},\tilde{e},\tilde{i}\right)\in \tilde{U}:{\psi }_1\le \tilde{s}<1,{\psi }_2\le \tilde{e}<1,{\psi }_3\le \tilde{i}<1,{\psi }_4\le \tilde{s}+\tilde{e}+\tilde{i}\le 1-{\psi }_4\right\},$$

where

$$\tilde{U}=\left\{0<\tilde{s}<1,0<\tilde{e}<1,0<\tilde{i}<1,0<\tilde{s}+\tilde{e}+\tilde{i}<1\right\}$$

and ${\psi }_i\in (0,1)(i=1,2,3,4)$ is a sufficiently small constant satisfying the following conditions:

$${\psi }_2={{\psi }_1}^2{\psi }_3,{\psi }_4={{\psi }_1}^2={{\psi }_3}^2,$$

(16)

$$-\frac{b}{{\psi }_1}+M_1\le -1,$$

(17)

$$-\frac{(1-\delta p)\beta }{{\psi }_1}+M_1\le -1,$$

(18)

$$(1-\delta p)\beta {\psi }_3-b<0,$$

(19)

$$-\frac{b}{{\psi }_4}+M_1\le -1,$$

(20)

$$-\frac{\delta p}{{\psi }_1}-\frac{\gamma }{{\psi }_3}+M_1\le -1.$$

(21)

Then

$$\tilde{U}\backslash U=U_1\cup U_2\cup U_3\cup U_4\cup U_5,$$

with

$$U_1=\left\{\left(\tilde{s},\tilde{e},\tilde{i}\right)\in \tilde{U}:0<\tilde{s}<{\psi }_1\right\},U_2=\left\{\left(\tilde{s},\tilde{e},\tilde{i}\right)\in \tilde{U}:{\psi }_1\le \tilde{s}<1,0<\tilde{e}<{\psi }_2,{\psi }_3\le \tilde{i}<1\right\},$$

$$U_3=\left\{\left(\tilde{s},\tilde{e},\tilde{i}\right)\in \tilde{U}:0<\tilde{i}<{\psi }_3\right\},U_4=\left\{\left(\tilde{s},\tilde{e},\tilde{i}\right)\in \tilde{U}:0<\tilde{s}+\tilde{e}+\tilde{i}<{\psi }_4\right\},$$

$$U_5=\left\{\left(\tilde{s},\tilde{e},\tilde{i}\right)\in \tilde{U}:{\psi }_1\le \tilde{s}<1,{\psi }_3\le \tilde{i}<1,1-{\psi }_4<\tilde{s}+\tilde{e}+\tilde{i}<1\right\}.$$

Now we prove the negativity of $LV$ for any $\tilde{U}\backslash U$.

Case I. If $\left(\tilde{s},\tilde{e},\tilde{i}\right)\in U_1$, it follows from (A1) and (17) that

$$\begin{array}{c}\hfill LV\le -\frac{b}{\tilde{s}}+M_1\le -\frac{b}{{\psi }_1}+M_1\le -1.\end{array}$$

Case II. If $\left(\tilde{s},\tilde{e},\tilde{i}\right)\in U_2$, (16) and (18) derive that

$$\begin{array}{c}\hfill LV\le -\frac{(1-\delta p)\beta \tilde{s}\tilde{i}}{\tilde{e}}+M_1\le -\frac{(1-\delta p)\beta {\psi }_1{\psi }_3}{{\psi }_2}+M_1=-\frac{(1-\delta p)\beta }{{\psi }_1}+M_1\le -1.\end{array}$$

Case III. If $\left(\tilde{s},\tilde{e},\tilde{i}\right)\in U_3$, (A1) and (19) yield that

$$\begin{array}{cc}\hfill LV\le & -b-\frac{b}{\tilde{s}}-\frac{(1-\delta p)\beta \tilde{s}\tilde{i}}{\tilde{e}}-\frac{\alpha \tilde{e}}{\tilde{i}}+(1-\delta p)\beta {\psi }_3+2\delta p+5b+\alpha +2\epsilon +2\gamma \hfill \\ & +2{\sigma }_1^2+2{\sigma }_2^2+2{\sigma }_3^2+\frac{1}{2}{\sigma }_4^2\hfill \\ \hfill \le & -3\sqrt[3]{b(1-\delta p)\beta \alpha }+2\delta p+5b+\alpha +2\epsilon +2\gamma +2{\sigma }_1^2+2{\sigma }_2^2+2{\sigma }_3^2+\frac{1}{2}{\sigma }_4^2<0.\hfill \end{array}$$

Case IV. If $\left(\tilde{s},\tilde{e},\tilde{i}\right)\in U_4$, (A1) and (20) imply that

$$\begin{array}{c}\hfill LV\le -\frac{b}{\tilde{s}+\tilde{e}+\tilde{i}}+M_1\le -\frac{b}{{\psi }_4}+M_1\le -1.\end{array}$$

Case V. If $\left(\tilde{s},\tilde{e},\tilde{i}\right)\in U_5$, it follows from (A1), (16) and (21) that

$$\begin{array}{c}\hfill LV\le -\frac{\delta p\tilde{s}}{\tilde{r}}-\frac{\gamma \tilde{i}}{\tilde{r}}+M_1\le -\frac{\delta p{\psi }_1}{{\psi }_4}-\frac{\gamma {\psi }_3}{{\psi }_4}+M_1=-\frac{\delta p}{{\psi }_1}-\frac{\gamma }{{\psi }_3}+M_1\le -1.\end{array}$$

Define

$$\varphi =\max \left\{-1,-3\sqrt[3]{b(1-\delta p)\beta \alpha }+2\delta p+5b+\alpha +2\epsilon +2\gamma +2{\sigma }_1^2+2{\sigma }_2^2+2{\sigma }_3^2+\frac{1}{2}{\sigma }_4^2\right\}<0.$$

Obviously, one can see that $LV\le \varphi <0$ for all $\left(\tilde{s},\tilde{e},\tilde{i}\right)\in \tilde{U}\backslash U$, which shows that Assumption 1$\left(ii\right)$ is satisfied. On the other hand, there exists a positive number

$$\begin{array}{cc}Z=\hfill & \min \{\left({\sigma }_1^2{\left(1-\tilde{s}\right)}^2+{\sigma }_2^2{\tilde{e}}^2+{\sigma }_3^2{\tilde{i}}^2+{\sigma }_4^2{\left(1-\tilde{s}-\tilde{e}-\tilde{i}\right)}^2\right){\tilde{s}}^2,({\sigma }_1^2{\tilde{s}}^2+{\sigma }_2^2{\left(1-\tilde{e}\right)}^2+{\sigma }_3^2{\tilde{i}}^2+{\sigma }_4^2\hfill \\ & \times {\left(1-\tilde{s}-\tilde{e}-\tilde{i}\right)}^2){\tilde{e}}^2,\left({\sigma }_1^2{\tilde{s}}^2+{\sigma }_2^2{\tilde{e}}^2+{\sigma }_3^2{\left(1-\tilde{i}\right)}^2+{\sigma }_4^2{\left(1-\tilde{s}-\tilde{e}-\tilde{i}\right)}^2\right){\tilde{i}}^2,\left(\tilde{s},\tilde{e},\tilde{i}\right)\in \tilde{U}\}\hfill \end{array}$$

such that

$$\begin{array}{cc}\hfill \sum _{i,j=1}^3{a}_{ij}{\xi }_i{\xi }_j=& \left({\sigma }_1^2{\left(1-\tilde{s}\right)}^2+{\sigma }_2^2{\tilde{e}}^2+{\sigma }_3^2{\tilde{i}}^2+{\sigma }_4^2{\left(1-\tilde{s}-\tilde{e}-\tilde{i}\right)}^2\right){\tilde{s}}^2{\xi }_1^2+({\sigma }_1^2{\tilde{s}}^2+{\sigma }_2^2{\left(1-\tilde{e}\right)}^2+{\sigma }_3^2{\tilde{i}}^2\hfill \\ & +{\sigma }_4^2{\left(1-\tilde{s}-\tilde{e}-\tilde{i}\right)}^2){\tilde{e}}^2{\xi }_2^2+\left({\sigma }_1^2{\tilde{s}}^2+{\sigma }_2^2{\tilde{e}}^2+{\sigma }_3^2{\left(1-\tilde{i}\right)}^2+{\sigma }_4^2{\left(1-\tilde{s}-\tilde{e}-\tilde{i}\right)}^2\right){\tilde{i}}^2{\xi }_3^2\hfill \\ \hfill \ge & Z{\left|\xi \right|}^2,\left(\tilde{s},\tilde{e},\tilde{i}\right)\in \tilde{U},\xi \in {\mathbb{R}}^3,\hfill \end{array}$$

which shows that Assumption 1$\left(i\right)$ is satisfied. Consequently, the system (6) has a unique stationary distribution $\pi (·)$ and it has ergodic property. The proof of Theorem 4 is complete. ☐

6. Simulations and Conclusions

6.1. Simulations

Next, in order to support the results of the above theorems, we carry out some computer simulations.

In Figure 1, take $\tilde{s}\left(0\right)=0.3$, $\tilde{e}\left(0\right)=0.25$, $\tilde{i}\left(0\right)=0.15$, $b=0.15$, $\beta =0.5$, $\gamma =0.3$, $\alpha =0.2$, $\delta =0.25$, $p=0.2$, $\epsilon =0.15$ and ${\sigma }_1={\sigma }_2={\sigma }_3={\sigma }_4=0.25$. Then

$$M_1=\max \left\{{\sigma }_1^2,\left(1-{\varrho }^2\right){\sigma }_3^2,{\sigma }_4^2\right\}=0.0625<2(1-\varrho )\alpha =0.0652$$

and

$$\epsilon =0.15<(1-\delta p)\beta =0.475<\epsilon +\gamma +b=0.6$$

satisfy the parameter conditions in Theorem 2, we can get that the exposed and infective individuals go to extinction almost surely. Obviously, Figure 1 validates our results of the Theorem 2.

In Figure 2, take $\tilde{s}\left(0\right)=0.15$, $\tilde{e}\left(0\right)=0.2$, $\tilde{i}\left(0\right)=0.15$, $b=0.02$, $\beta =0.9$, $\gamma =0.01$, $\alpha =0.1$, $\delta =0.02$, $p=0.02$, $\epsilon =0.16$, ${\sigma }_1=0.05$, ${\sigma }_2=0.05$, ${\sigma }_3=0.05$ and ${\sigma }_4=0.1$. Obviously,

$$0.1216=\sqrt[3]{b(1-\delta p)\beta \alpha }>\frac{\delta p+3b+\alpha +\epsilon +\gamma +\frac{1}{2}{\sigma }_1^2+\frac{1}{2}{\sigma }_2^2+\frac{1}{2}{\sigma }_3^2}{3}=0.1114$$

satisfies the parameter condition in Theorem 3, then

$$0.0342=\underline{\tilde{i}}\le \underset{t\to \infty }{\lim \inf }〈\tilde{i}\left(t\right)〉\le \underset{t\to \infty }{\lim \sup }〈\tilde{i}\left(t\right)〉\le \overline{\tilde{i}}=0.9766,$$

we can get that the infective individuals $\tilde{i}\left(t\right)$ are permanent in mean almost surely. As expected, Figure 2 confirms our results of the Theorem 3.

From Figure 2 and Figure 3, a set of large stochastic parameter values ${\sigma }_1={\sigma }_2={\sigma }_3={\sigma }_4=0.25$ can lead to infective individuals go to extinction (see Figure 2), while infective individuals can be permanent in mean under the condition of a set of small stochastic parameter values ${\sigma }_1={\sigma }_2={\sigma }_3=0.05$ and ${\sigma }_4=0.1$ (see Figure 3).

In Figure 3, take $\tilde{s}\left(0\right)=0.15$, $\tilde{e}\left(0\right)=0.2$, $\tilde{i}\left(0\right)=0.15$, $b=0.02$, $\beta =2.1$, $\gamma =0.01$, $\alpha =0.2$, $\delta =0.02$, $p=0.02$, $\epsilon =0.1$ and ${\sigma }_1={\sigma }_2={\sigma }_3={\sigma }_4=0.01$. Then

$$0.2033=\sqrt[3]{b(1-\delta p)\beta \alpha }>\frac{2\delta p+5b+\alpha +2\epsilon +2\gamma +2{\sigma }_1^2+2{\sigma }_2^2+2{\sigma }_3^2+\frac{1}{2}{\sigma }_4^2}{3}=0.1738$$

satisfies the parameter condition in Theorem 4, we can get that the stochastic system (6) has a unique stationary distribution $\pi (·)$ and it has ergodic property. Figure 3 indicates that the solution of system (6) swings up and down in a small neighborhood. According to the density functions in Figure 3d–f, we can see that there exists a stationary distribution. As expected, Figure 3 supports our results of the Theorem 4.

The Figure 1, Figure 2 and Figure 3 above show that the large white noise value can lead to infectious diseases to go to extinction, which implies that stochastic fluctuations can suppress the disease outbreak, while the small white noise value can cause infectious diseases to be persistent. In addition, The Figure 3 also shows the stochastic system (6) has a unique ergodic stationary distribution under appropriate conditions. Therefore, the numerical simulation examples are completely consistent with the theoretical results of the Theorems 2–4.

6.2. Conclusions

In this paper, we apply stochastic analysis methods to study the global dynamics of a high-dimensional stochastic reduced proportional SEIR epidemic system which makes the analysis novel and complex. We obtain the existence of a unique global positive solution and parameter conditions of extinction or permanence in mean. Furthermore, the solution of the stochastic system has a unique ergodic stationary distribution under certain sufficient parameter conditions. Cubic terms of $\tilde{s},\tilde{e},\tilde{i}$ and multiple stochastic terms for $d{B}_i\left(t\right)(i=1,2,3,4)$ in system (6) make the analysis more difficult and complex than the models in [34,36]. Some ingenious inequality techniques are used to deal with cubic terms of $\tilde{s},\tilde{e},\tilde{i}$ of system (6). Therefore, compare with previous methods and research results, we develop previous methods and improve the main results of previous studies.

We summarize the main conclusions as follows:

(I) When

$$M_1<2(1-\varrho )\alpha ,\epsilon \le (1-\delta p)\beta <\epsilon +\gamma +b$$

hold, then

$$\underset{t\to \infty }{\lim \sup }\frac{\ln \left(\tilde{e}\left(t\right)+\varrho \tilde{i}\left(t\right)\right)}{t}<0a.s.$$

That is to say, the exposed and infective individuals go to extinction almost surely.

(II) When

$$\sqrt[3]{b(1-\delta p)\beta \alpha }>\frac{\delta p+3b+\alpha +\epsilon +\gamma +\frac{1}{2}{\sigma }_1^2+\frac{1}{2}{\sigma }_2^2+\frac{1}{2}{\sigma }_3^2}{3}$$

holds, then

$$\underline{\tilde{i}}\le \underset{t\to \infty }{\lim \inf }〈\tilde{i}\left(t\right)〉\le \underset{t\to \infty }{\lim \sup }〈\tilde{i}\left(t\right)〉\le \overline{\tilde{i}}a.s.$$

That is to say, the infective individuals $\tilde{i}\left(t\right)$ are permanent in mean almost surely.

(III) When

$$\sqrt[3]{b(1-\delta p)\beta \alpha }>\frac{2\delta p+5b+\alpha +2\epsilon +2\gamma +2{\sigma }_1^2+2{\sigma }_2^2+2{\sigma }_3^2+\frac{1}{2}{\sigma }_4^2}{3}$$

holds, then the system (6) has a unique stationary distribution $\pi (·)$ and it has ergodic property.

By comparing the above conclusions (II) and (III), we can see that when system (6) has a ergodic stationary distribution, then the infective individuals $\tilde{i}\left(t\right)$ are permanent in mean almost surely. However, it is not applicable in reverse. The above results of Theorems 2–4 show a large stochastic disturbance can cause infectious diseases to go to extinction, in other words, the persistent infectious disease of a deterministic system can become extinct due to the white noise stochastic disturbance. This implies that stochastic fluctuations can suppress the disease outbreak.