Dynamic Properties for a Second-Order Stochastic SEIR Model with Infectivity in Incubation Period and Homestead-Isolation of the Susceptible Population

Abstract

In this article, we analyze a second-order stochastic SEIR epidemic model with latent infectious and susceptible populations isolated at home. Firstly, by putting forward a novel inequality, we provide a criterion for the presence of an ergodic stationary distribution of the model. Secondly, we establish sufficient conditions for extinction. Thirdly, by solving the corresponding Fokker–Plank equation, we derive the probability density function around the quasi-endemic equilibrium of the stochastic model. Finally, by using the epidemic data of the corresponding deterministic model, two numerical tests are presented to illustrate the validity of the theoretical results. Our conclusions demonstrate that nations should persevere in their quarantine policies to curb viral transmission when the COVID-19 pandemic proceeds to spread internationally.

1. Introduction

As everyone knows, Corona-virus 2019 (COVID-19) has become a serious threat to human health and lives in recent years. On 25 January 2022, the World Health Organization (WHO) declared that nearly 50,000 new deaths were reported. As of 23 January 2022, 634 million confirmed cases and 6.6 million deaths have been reported worldwide [1]. It is worth emphasizing that social isolation has had a direct effect on the physical transmission of infectious diseases [2]. The COVID-19 outbreak clearly showed that physical protection and social isolation play a key role in controlling epidemics when vaccines or antiviral drugs for the virus are lacking. Mathematical modeling is an important promoting factor in discussing the dynamic behaviors of a disease’s spread [3]. Based on the importance of social isolation, Jiao et al. [4] developed an SEIR model with infectivity in the incubation period and homestead-isolation of the susceptible population. The model is given as follows:

$$\left\{\begin{array}{c}\frac{dS}{dt}=\Lambda -\beta (1-{\theta }_1)S[I+{\theta }_{2}E]-\mu S,\hfill \\ \frac{dE}{dt}=\beta (1-{\theta }_1)S[I+{\theta }_{2}E]-(\delta +\mu )E,\hfill \\ \frac{dI}{dt}=\delta E-(\gamma +\sigma +\mu )I,\hfill \\ \frac{dR}{dt}=(\gamma +{\theta }_3\sigma )I-\mu R,\hfill \end{array}\right.$$

(1)

where S is the numbers of susceptible people, E denotes the numbers of the exposed population, I stands for the numbers of the infected population, $R\left(t\right)$ denotes the numbers of the recovered population at time t. $\Lambda >0$ signifies the enrolling rate, $\beta >0$ refers to the infection rate from S to E, $0<{\theta }_1<1$ indicates the homestead-isolation rate of the susceptible population, $0<{\theta }_2<1$ represents the infectious effect of the exposed population in incubation period, $\mu >0$ is the natural death rate, $\delta >0$ is the transition rate from E to I, $\gamma >0$ is the transition rate from I to R, $\sigma >0$ represents the hospitalized rate of I for the disease, ${\theta }_3>0$ depicts the recurring rate of I, and $\delta >{\theta }_2(\gamma +\sigma +\mu )$.

It is obvious that the individual $R\left(t\right)$ does not adversely affect the dynamics of SEIR. Therefore, we will neglect the last equation of model (1), and it can be reducible to:

$$\left\{\begin{array}{c}\frac{dS}{dt}=\Lambda -\beta (1-{\theta }_1)S[I+{\theta }_{2}E]-\mu S,\hfill \\ \frac{dE}{dt}=\beta (1-{\theta }_1)S[I+{\theta }_{2}E]-(\delta +\mu )E,\hfill \\ \frac{dI}{dt}=\delta E-(\gamma +\sigma +\mu )I.\hfill \end{array}\right.$$

(2)

Then, Jiao et al. [4] derived that model (1) demonstrates a disease-free equilibrium $P^0\left(S^0,0,0\right)$ with ${\displaystyle S^0=\frac{\Lambda }{\mu }}$ and an interior endemic equilibrium $P^{+}\left(S^{+},E^{+},I^{+}\right)$, where

$${\displaystyle S^{+}=\frac{(\gamma +\sigma +\mu )(\delta +\mu )}{\beta (1-\theta )\left[\delta +{\theta }_2(\gamma +\sigma +\mu )\right]},E^{+}=\frac{\Lambda \beta \left(1-{\theta }_1\right)\left[\delta +{\theta }_2(\gamma +\sigma +\mu )\right]-\mu (\gamma +\sigma +\mu )(\delta +\mu )}{\beta \left(1-{\theta }_1\right)(\delta +{\theta }_2(\gamma +\sigma +\mu ))(\delta +\mu )},}$$

${\displaystyle I^{+}=\frac{\delta E^{+}}{\gamma +\sigma +\mu }}$ with $\Lambda \beta (1-{\theta }_1)[\delta +{\theta }_2(\gamma +\sigma +\mu )]>\mu (\gamma +\sigma +\mu )(\delta +\mu )$, and the basic reproduction number is defined as ${\displaystyle {\mathfrak{R}}_0=\frac{\Lambda \beta \left(1-{\theta }_1\right)\left[\delta +{\theta }_2(\gamma +\sigma +\mu )\right]}{\mu (\gamma +\sigma +\mu )(\delta +\mu )}}$.

In practical terms, the process of disease transmission is frequently disturbed by external factors. Mao et al. have claimed that the incorporation of a small amount of random fluctuation in the deterministic model of infectious disease can actively prevent outbreaks in the latent population [5]. Recently, several authors have considered the random fluctuation factor into biological mathematical models of infectious diseases and environmental pollution [6,7,8,9,10]. In [11], in view of nonlinear stochastic disturbance, Mu et al. discussed the extinction and persistence for a stochastic micro-organism flocculation model. In [12], for a stochastic staged model of the progression of HIV/AIDS with nonlinear disturbances, Jiang et al. first presented a new approach called the “stochastic p-threshold method” to remove nonlinear terms caused by second-order disturbances and obtained the unique ergodic stationary distribution of the model. In [13], Zhou et al. made use of the “p-stochastic criterion technique” together with some inequalities they presented to eliminate the effect of third-order perturbation of a stochastic model. Motivated by these works, we introduce nonlinear stochastic disturbances into model (2), which is then expressed by:

$$\left\{\begin{array}{c}dS=[\Lambda -\beta (1-{\theta }_1)S(I+{\theta }_{2}E)-\mu S]dt+({\sigma }_{11}+{\sigma }_{12}S)Sd{B}_1\left(t\right),\hfill \\ dE=[\beta (1-{\theta }_1)S(I+{\theta }_{2}E)-(\delta +\mu )E]dt+({\sigma }_{21}+{\sigma }_{22}E)Ed{B}_2\left(t\right),\hfill \\ dI=[\delta E-(\gamma +\sigma +\mu )I]dt+({\sigma }_{31}+{\sigma }_{32}I)Id{B}_3\left(t\right),\hfill \end{array}\right.$$

(3)

where ${\left\{B_i\left(t\right)\right\}}_{t\ge 0}(i=1,2)$ denote mutually independent standard Brownian motions and are defined on the complete probability space $\{\Omega ,\Gamma ,{\left\{{\Gamma }_t\right\}}_{t\ge 0},\mathbb{P}\}$ with an increasing and right continuous filtration ${\left\{{\Gamma }_t\right\}}_{t\ge 0}$; the positive constants ${\sigma }_{ij}(i=1,2,3,j=1,2)$ are the white noise intensities.

Notably, relatively few discussions exist deriving the explicit expression of the probability density function as a result of the difficulty of solving the high-dimensional Fokker–Plank equation. In [14], Zhou et al. investigated the density function analysis of a stochastic SVI epidemic model with a half-saturated incidence rate for the first time. In light of this work, some researchers are looking into the properties of stochastic epidemic models [2,6,15,16,17,18,19,20]. However, no authors have studied the exact expression of the density function of model (3) yet, which is the key work of the present paper. In addition, the aforementioned “stochastic p-threshold method” is not completely suitable for the corresponding deterministic case with a complicated basic reproduction number when researching the unique ergodic stationary distribution of a stochastic epidemic disease model. Therefore, another aim of this paper is to improve the “stochastic p-threshold method” and explore the unique ergodic stationary distribution and extinction for model (3).

The remainder of this paper is arranged in the following manner. In Section 2, we investigate the stationary distribution for model (3). In Section 3, we give adequate conditions of extinction for model (3). Section 4 derives a precise formula for the probability density function of model (3). In Section 5, two examples are presented to demonstrate the practicality of our conclusions.

2. Existence of Ergodic Stationary Distribution

In the whole paper, let ${\mathbb{R}}^n$ be an n-dimensional Euclidean space: ${\mathbb{R}}_{+}^n=\{(y_1,\dots ,y_n)\in {\mathbb{R}}^n|y_k>0,0\le k\le n\},l_1\wedge l_2\wedge \dots \wedge l_n=min\{l_1,l_2,\dots ,l_n\},h_1\vee h_2\vee \dots \vee h_n=max\{h_1,h_2,\dots ,h_n\}$. If $B_{s\times t}$ is a real matrix, let $B^T$ be the transpose matrix of B. When $s=t$, $B^{-1}$ represents the inverse matrix of B.

Lemma 1.

For any initial data $(S_0,E_0,I_0)\in {\mathbb{R}}_{+}^3$, there is a unique positive solution $\left(S\right(t),E(t),I(t\left)\right)$ of model (3) on $t\ge 0$, and the solution will remain in ${\mathbb{R}}_{+}^3$ with a probability of one (a.s.).

The proof of Lemma 1 is given in the following Appendix A.

Consider $Y\left(t\right)$ to be a regular time-homogeneous Markov process in ${\mathbb{R}}^d$ satisfying the following SDE:

$$\begin{array}{c}\hfill dY\left(t\right)=f\left(Y\left(t\right)\right)dt+\sum _{l=1}^k{g}_l\left(Y\left(t\right)\right)d{B}_l\left(t\right).\end{array}$$

(4)

Its diffusion matrix is:

$$\begin{array}{c}\hfill A\left(y\right)=\left({\widehat{a}}_{ij}\left(y\right)\right),{\widehat{a}}_{ij}\left(y\right)=\sum _{l=1}^k{g}_l^i\left(y\right)g_l^j\left(y\right).\end{array}$$

Lemma 2

([17,21,22]). Assume that there exists a bounded domain $\mathbb{D}\subset {\mathbb{R}}^d$ with the regular boundary Γ and

$\left(A_1\right)$∃$L>0$, such that ${\Sigma }_{i,j}^d={\widehat{a}}_{ij}\left(y\right){\widehat{\xi }}_i{\widehat{\xi }}_j,y\in \mathbb{D},\widehat{\xi }\in {\mathbb{R}}^d$;

$\left(A_2\right)$∃$C^2$-function $V\left(y\right)\ge 0$ and $\exists \varpi >0$ such that $\mathcal{L}V\le -\varpi$, where $\mathcal{L}$ denotes the Itô’s differential operator defined in (45) of Ref. [13].

Then, system (4) exists a unique ergodic stationary distribution $\pi (·)$ for any ${\mathbb{R}}_{+}^d$, such that

$${\mathbb{P}}_y\left\{\underset{T\to \infty }{lim}\frac{1}{T}{\int }_0^{T}f\left(Y\left(t\right)\right)dt={\int }_{{\mathbb{R}}^d}f\left(y\right)\pi \left(dy\right)\right\}=1,$$

and let $f(·)$ be an integrable function with respect to the measure $\pi (·)$.

Lemma 3.

If $x\ge 0$, then

$$x^n\ge \left(\frac{n-1}{2(n-2)}{x}^{n-2}-\frac{1}{2(n-2)}\right)(x^2+1),$$

where $n\ge 3$ and the sign of inequality holds if and only if $x=1$.

The proof of Lemma 3 is given in Appendix A.

Theorem 1.

Assume that ${\mathfrak{R}}_0^S>1$, where

$$\begin{array}{cc}\hfill {\mathfrak{R}}_0^S=& \frac{\beta (1-{\theta }_1)\Lambda \delta }{\left(\mu +\frac{{\sigma }_{11}^2}{2}+2\sqrt{\Lambda {\sigma }_{11}{\sigma }_{12}}+2\sqrt[3]{{\Lambda }^2{\sigma }_{12}^2}\right)\left(\gamma +\sigma +\mu +\frac{{\sigma }_{31}^2}{2}\right)\left(\delta +\mu +\frac{{\sigma }_{21}^2}{2}+2\sqrt[3]{{\Lambda }^2{\sigma }_{22}^2}\right)}\hfill \\ \hfill & +\frac{\beta (1-{\theta }_1){\theta }_2\Lambda }{\left(\mu +\frac{{\sigma }_{11}^2}{2}+2\sqrt{\Lambda {\sigma }_{11}{\sigma }_{12}}+2\sqrt[3]{{\Lambda }^2{\sigma }_{12}^2}\right)\left(\delta +\mu +\frac{{\sigma }_{21}^2}{2}+2\sqrt[3]{{\Lambda }^2{\sigma }_{22}^2}\right)};\hfill \end{array}$$

then, for any initial data $(S\left(0\right),E\left(0\right),I\left(0\right))\in {\mathbb{R}}_{+}^3$, model (3) possesses a unique ergodic stationary distribution.

The proof of Theorem is given in Appendix A.

Remark 1.

By combining setting up constant $C_1+C_2$ in $V_1$ with “stochastic p-threshold method”, we successfully proved Theorem 1. The new method can make up for the deficiency of the common “stochastic p-threshold method” to a certain extent.

Remark 2.

We can see that ${\mathfrak{R}}_0^S={\mathfrak{R}}_0$ when ${\sigma }_{ij}^2\equiv 0(i=1,2,3,j=1,2)$ in model (3).

Remark 3.

From ${\mathfrak{R}}_0^S$ in Theorem 1, we can come to the conclusion that ${\theta }_1$ plays an active role in the prevention and control of an infectious disease.

3. Extinction

Note that

$$\begin{array}{c}\hfill d\tilde{A}=\left(\Lambda -\mu \tilde{A}\right)dt+\left({\sigma }_{11}\tilde{A}+{\sigma }_{12}{\tilde{A}}^2\right)d{B}_1\left(t\right),\end{array}$$

(5)

with the same initial data $\tilde{A}\left(0\right)=S\left(0\right)$. According to Ref. [23], model (5) exists as a stationary solution with the density

$$\begin{array}{cc}\hfill \pi \left(x\right)=& Q{x}^{-2-\frac{2\left(2\Lambda {\sigma }_{12}+\mu {\sigma }_{11}\right)}{{\sigma }_{11}^3}}{\left({\sigma }_{11}+{\sigma }_{12}x\right)}^{-2+\frac{2\left(2\Lambda {\sigma }_{12}+\mu {\sigma }_{11}\right)}{{\sigma }_{11}^3}}exp\left(-\frac{2\Lambda }{{\sigma }_{11}({\sigma }_{11}+{\sigma }_{12}x)}\right)\left(\frac{\Lambda }{x}\right.\hfill \\ \hfill & \left.+\frac{2\Lambda {\sigma }_{12}+\mu {\sigma }_{11}}{{\sigma }_{11}}\right),x\in (0,+\infty ),\hfill \end{array}$$

where

$${\displaystyle Q=\frac{1}{{\int }_0^{\infty }{x}^{-2-\frac{2\left(2\Lambda {\sigma }_{12}+\mu {\sigma }_{11}\right)}{{\sigma }_{11}^3}}{\left({\sigma }_{11}+{\sigma }_{12}x\right)}^{-2+\frac{2\left(2\Lambda {\sigma }_{12}+\mu {\sigma }_{11}\right)}{{\sigma }_{11}^3}}exp\left(-\frac{2\Lambda }{{\sigma }_{11}({\sigma }_{11}+{\sigma }_{12}x)}\right)\left(\frac{\Lambda }{x}\right.\left.+\frac{2\Lambda {\sigma }_{12}+\mu {\sigma }_{11}}{{\sigma }_{11}}\right)dx}}$$

is a constant such that ${\displaystyle {\int }_0^{\infty }\pi \left(x\right)dx=1.}$ Define

$${\mathfrak{R}}_0^E={\mathfrak{R}}_0{\rho }_0{\int }_0^{\infty }\left|x-\frac{\Lambda }{\mu }\right|\pi \left(x\right)dx-\frac{{\sigma }_{22}^2\wedge {\sigma }_{32}^2}{4},$$

where ${\displaystyle {\rho }_0=\frac{\beta (1-{\theta }_1)}{{\rho }_1\wedge {\rho }_2}}$.

Theorem 2.

Let $\left(S\right(t),E(t),I(t\left)\right)$ be the solution to model (3) with any initial data $(S_0,E_0,I_0)\in {\mathbb{R}}_{+}^3$. If ${\mathfrak{R}}_0<1,{\mathfrak{R}}_0^E<0$, then

$$\underset{t\to \infty }{lim}sup\frac{1}{t}ln({\rho }_{1}E\left(t\right)+{\rho }_{2}I\left(t\right))\le {\mathfrak{R}}_0^E<0,a.s.,$$

and

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

where the positive constants ${\rho }_1,{\rho }_2$ are defined as follows:

$$\begin{array}{c}\hfill {\rho }_1={\mathfrak{R}}_0,{\rho }_2=\frac{\Lambda \beta (1-{\theta }_1)}{\mu (\gamma +\sigma +\mu )}.\end{array}$$

The proof of Theorem 2 is given in Appendix A.

Remark 4.

In view of ${\mathfrak{R}}_0^E$ in Theorem 2, a larger ${\theta }_1$ is favorable to restraining the number of infected individuals and eradicating the disease at the beginning of its outbreak.

4. Density Function Analysis of Model (3)

Here, we discuss the representation of the local probability density function of model (3) under only linear stochastic perturbations.

Let $x_1=lnS,x_2=lnE,x_3=lnI$. For Equation (3), we have:

$$\left\{\begin{array}{cc}\hfill d{x}_1=& \left(\Lambda e^{-x_1}-\beta (1-{\theta }_1)(e^{x_3}+{\theta }_2{e}^{x_2})-{\mu }_1\right)dt+{\sigma }_{11}d{B}_1\left(t\right),\hfill \\ \hfill d{x}_2=& \left(\beta (1-{\theta }_1)e^{x_1}(e^{x_3}+{\theta }_2{e}^{x_2})e^{-x_2}-(\delta +{\mu }_2)\right)dt+{\sigma }_{21}d{B}_2\left(t\right),\hfill \\ \hfill d{x}_3=& (\delta e^{x_2-x_3}-(\gamma +\sigma +{\mu }_3))dt+{\sigma }_{31}d{B}_3\left(t\right),\hfill \end{array}\right.$$

(6)

where ${\displaystyle {\mu }_i=\mu +\frac{{\sigma }_{i1}^2}{2},i=1,2,3}$. Assume that $\widehat{{\mathfrak{R}}_0^S}>1$; then, the unique quasi-endemic equilibrium $E^{*}=(S^{*},E^{*},I^{*})=(e^{x_1^{*}},e^{x_2^{*}},e^{x_3^{*}})$ satisfies:

$$\begin{array}{cc}\hfill & \left\{\begin{array}{c}\Lambda e^{-x_1^{*}}-\beta (1-{\theta }_1)(e^{x_3^{*}}+{\theta }_2{e}^{x_2^{*}})-{\mu }_1=0,\hfill \\ \beta (1-{\theta }_1)e^{x_1}(e^{x_3^{*}}+{\theta }_2{e}^{x_2^{*}})e^{-x_2^{*}}-(\delta +{\mu }_2)=0,\hfill \\ \delta e^{x_2^{*}}{e}^{-x_3^{*}}-(\gamma +\sigma +{\mu }_3)=0,\hfill \end{array}\right.\hfill \end{array}$$

where

$$\begin{array}{cc}\hfill & S^{*}=\frac{\left(\gamma +\sigma +{\mu }_3\right)\left(\delta +{\mu }_2\right)}{\beta \left(1-{\theta }_1\right)\left[\delta +{\theta }_2\left(\gamma +\sigma +{\mu }_3\right)\right]}>0,E^{*}=\frac{{\mu }_1\left(\delta +{\mu }_2\right)\left(\gamma +\sigma +{\mu }_3\right)(\widehat{{\mathfrak{R}}_0^S}-1)}{\beta \left(1-{\theta }_1\right)\left[\delta +{\theta }_2\left(\gamma +\sigma +{\mu }_3\right)\right]\left(\delta +{\mu }_2\right)}>0,\hfill \\ \hfill & I^{*}=\frac{\delta E^{*}}{\gamma +\sigma +{\mu }_3}>0,\widehat{{\mathfrak{R}}_0^S}=\frac{\Lambda \beta (1-{\theta }_1)[\delta +{\theta }_2(\gamma +\sigma +{\mu }_3)]}{{\mu }_1(\delta +{\mu }_2)(\gamma +\sigma +{\mu }_3)}.\hfill \end{array}$$

Evidently, $(S^{*},E^{*},I^{*})$ coincides with the stable endemic equilibrium point $(S^{+},E^{+},I^{+})$ of the deterministic model (1) when ${\sigma }_{i1}=0,i=1,2,3$, $\widehat{{\mathfrak{R}}_0^S}={\mathfrak{R}}_0^S$ if ${\sigma }_{11}={\sigma }_{22}=0$.

Let $(y_1,y_2,y_3)=(x_1-x_1^{*},x_2-x_2^{*},x_3-x_3^{*})$, where $x_1^{*}=ln{S}^{*},x_2^{*}=ln{E}^{*},x_3^{*}=ln{I}^{*}$. The linearized model of Equation (6) is as follows:

$$\left\{\begin{array}{cc}\hfill d{y}_1=& (-a_{11}{y}_1-a_{12}{y}_2-a_{13}{y}_3)dt+{\sigma }_{1}d{B}_1\left(t\right),\hfill \\ \hfill d{y}_2=& (a_{21}{y}_1-a_{22}{y}_2+a_{23}{y}_3)dt+{\sigma }_{2}d{B}_2\left(t\right),\hfill \\ \hfill d{y}_3=& (a_{32}{y}_2-a_{33}{y}_3)dt+{\sigma }_{3}d{B}_3\left(t\right),\hfill \end{array}\right.$$

(7)

where $a_{11}=\Lambda e^{-x_1^{*}}>0,a_{12}=\beta (1-{\theta }_1){\theta }_2{e}^{x_2^{*}}>0,a_{13}=\beta (1-{\theta }_1)e^{x_3^{*}}>0,a_{21}=\beta (1-{\theta }_1)e^{x_1^{*}}{e}^{x_3^{*}-x_2^{*}}+\beta (1-{\theta }_1){\theta }_2{e}^{x_1^{*}}>0,a_{22}=\beta (1-{\theta }_1)e^{-x_2^{*}}{e}^{x_1^{*}+x_3^{*}}>0,a_{23}=\beta (1-{\theta }_1)e^{x_1^{*}-x_2^{*}}{e}^{x_3^{*}}>0,a_{32}=\delta e^{x_2^{*}-x_3^{*}}>0,a_{33}=\delta e^{x_2^{*}-x_3^{*}}>0$. Furthermore, letting $Y={(y_1,y_2,y_3)}^T,\Lambda =diag({\sigma }_1,{\sigma }_2,{\sigma }_3)$, and $B\left(t\right)={(B_1\left(t\right),B_2\left(t\right),B_3\left(t\right))}^T$, then model (7) can be rewritten into $dY=AYdt+\Lambda dB\left(t\right)$, where

$$A=\left(\begin{array}{ccc}-a_{11}& -a_{12}& -a_{13}\\ a_{21}& -a_{22}& a_{23}\\ 0& a_{32}& -a_{33}\end{array}\right).$$

Lemma 4

([15]). For ${\mathrm{G}}_0^2+A_0{\Sigma }_0+{\Sigma }_0{A}_0^T=0$, where $G_0=diag(1,0,0)$,

$$A_0=\left(\begin{array}{ccc}-a_1& -a_2& -a_3\\ 1& 0& 0\\ 0& 1& 0\end{array}\right),$$

(8)

$${\Sigma }_0=\left(\begin{array}{ccc}\frac{a_2}{2(a_1{a}_2-a_3)}& 0& -\frac{1}{2(a_1{a}_2-a_3)}\\ 0& \frac{1}{2(a_1{a}_2-a_3)}& 0\\ -\frac{1}{2(a_1{a}_2-a_3)}& 0& \frac{a_1}{2(a_1{a}_2-a_3)}\end{array}\right).$$

If $a_1>0,a_3>0$ and $a_1{a}_2-a_3>0$, then matrix ${\Sigma }_0$ is positive definite.

Lemma 5

([15]). For $G_0^2+B_0{\theta }_0+{\theta }_0{B}_0^T=0$, where $G_0=diag(1,0,0)$,

$$B_0=\left(\begin{array}{ccc}-b_1& -b_2& -b_3\\ 1& 0& 0\\ 0& 0& b_{33}\end{array}\right).$$

(9)

If $b_1>0$ and $b_2>0$, then ${\theta }_0$ is semi-positive definite, which takes the form

$${\theta }_0=\left(\begin{array}{ccc}\frac{1}{2b_1}& 0& 0\\ 0& \frac{1}{2b_1{b}_2}& 0\\ 0& 0& 0\end{array}\right).$$

(10)

Theorem 3.

Let $Y=(y_1,y_2,y_3)$ be a solution to model (3) with any initial data $(y_1\left(0\right),y_2\left(0\right),y_3\left(0\right))\in {\mathbb{R}}_3$. If $\widehat{{\mathfrak{R}}_0^S}>1,r_1{r}_2-r_3>0$, where $r_1=a_{11}+a_{22}+a_{32},r_2=a_{11}{a}_{32}+a_{11}{a}_{22}+a_{12}{a}_{21},r_3=a_{13}{a}_{21}{a}_{32}+a_{32}{a}_{12}{a}_{21}$, then there exists a unique density function $\phi (S,E,I)$ around the quasi-stable equilibrium point $E^{*}$, where

$$\phi (S,E,I)={\left(2\pi \right)}^{-\frac{3}{2}}{|\Sigma |}^{-\frac{1}{2}}{e}^{-\frac{1}{2}(ln\frac{S}{S^{*}},ln\frac{E}{E^{*}},ln\frac{I}{I^{*}}){\Sigma }^{-1}{(ln\frac{S}{S^{*}},ln\frac{E}{E^{*}},ln\frac{I}{I^{*}})}^T},$$

in which $\Sigma ={\Sigma }_1+{\Sigma }_2+{\Sigma }_3$ is positive definite. ${\Sigma }_1$ is described in (A28), ${\Sigma }_2$ is depicted in (A30) and, for different cases, ${\Sigma }_3$ is described in (A32) and (A33).

The proof of Theorem 3 is given in Appendix A.

Remark 5.

If $\widehat{{\mathfrak{R}}_0^S}>1$, Theorem 3 shows that all of the distributions of the subpopulations $S\left(t\right)$, $E\left(t\right)$ and $I\left(t\right)$ will separately converge to the corresponding stationary marginal distributions ${\mu }_1\left(S\right),{\mu }_2\left(E\right),{\mu }_3\left(I\right)$ of $\varpi (·)$ as $t\to \infty$. By letting $\Sigma ={\left({\rho }_{ij}\right)}_{3\times 3}$, we then combine Theorem 3 to obtain that the distribution ${\mu }_1\left(S\right)$ around $S^{*}$ has an approximately log-normal density function $I_1\left(S\right)$. Moreover, the distribution ${\mu }_2\left(E\right)$ around $E^{*}$ has an approximately log-normal density function $I_2\left(E\right)$, and the distribution ${\mu }_3\left(I\right)$ around $I^{*}$ has an approximately log-normal density function $I_3\left(I\right)$, where

$$I_1\left(S\right)=\frac{1}{S\sqrt{2\pi {\rho }_{11}}}{e}^{-\frac{{(lnS-ln{S}^{*})}^2}{2{\rho }_{11}}},I_2\left(E\right)=\frac{1}{E\sqrt{2\pi {\rho }_{22}}}{e}^{-\frac{{(lnE-ln{E}^{*})}^2}{2{\rho }_{22}}},$$

$$I_3\left(I\right)=\frac{1}{I\sqrt{2\pi {\rho }_{33}}}{e}^{-\frac{{(lnI-ln{I}^{*})}^2}{2{\rho }_{33}}}.$$

5. Numerical Tests

Now we present two numerical tests and several figures to support our findings.

Numerical Test 1. To make use of model (3) to investigate the spread of COVID-19, in system (1), we use the same data as in Ref. [4]: $\Lambda =10,\mu =0.3,\sigma =0.2,\beta =0.2,\gamma =0.2,\delta =0.3,{\theta }_1=0.7,{\theta }_2=0.1$. This implies that the $P^{+}$ of model (1) is globally asymptotically stable. Now, let ${\sigma }_{11}=0.01$, ${\sigma }_{12}=0.0008$, ${\sigma }_{21}=0.01$, ${\sigma }_{22}=0.002$, ${\sigma }_{31}=0.019$ and ${\sigma }_{32}=0.02$. According to our calculation, ${\mathfrak{R}}_0^S=1.0661>1$. From Theorem 1, system (3) has a unique ergodic stationary distribution (see Figure 1).

Additionally, let ${\sigma }_{11}=0.01,{\sigma }_{21}=0.01,{\sigma }_{31}=0.019,{\sigma }_{12}={\sigma }_{22}={\sigma }_{32}=0$; we derive the equilibrium point $P^{*}=(18.925,7.201,3.086),\widehat{{\mathfrak{R}}_0^S}=1.761,r_1{r}_2-r_3=1.024>0$. Then, the criteria of Theorem 3 hold. According to Theorem 3, one can see that

$$\Sigma =\left(\begin{array}{ccc}0.0001567& -0.0000231& 0.0000905\\ -0.0000231& 0.0005228& -0.0005173\\ 0.0000905& -0.0005173& 0.0007306\end{array}\right).$$

$\phi (S,E,I)$ of model (3) is deduced as

$$\phi (S,E,I)=16174e^{-\frac{1}{2}(ln\frac{S}{18.925},ln\frac{E}{7.201},ln\frac{I}{3.086}){\Sigma }^{-1}{(ln\frac{S}{18.925},ln\frac{E}{7.201},ln\frac{I}{3.086})}^T},$$

and

$$I_1\left(S\right)=\frac{31.8776}{S}{e}^{-\frac{{(lnS-2.9405)}^2}{0.0003134}},I_2\left(E\right)=\frac{17.4523}{E}{e}^{-\frac{{(lnE-1.9742)}^2}{0.0010456}},$$

$$I_3\left(I\right)=\frac{14.7632}{I}{e}^{-\frac{{(lnI-1.1269)}^2}{0.0014612}}.$$

The right-hand column in Figure 1 presents the frequency histograms and the marginal densities.

Numerical Test 2. The same data are presented in Ref. [4]: $\Lambda =10,\mu =0.3,$ $\sigma =0.2,$$\beta =0.2,\gamma =0.2,\delta =0.3,$${\theta }_1=0.9,$${\theta }_2=0.1$. This shows that the disease-free equilibrium $P^0$ of model (2) is globally asymptotically stable. Now, set the stronger white noise intensities ${\sigma }_{11}=0.5$, ${\sigma }_{12}=0.0001,{\sigma }_{21}=2.1,{\sigma }_{22}=0.9487,{\sigma }_{31}=0.01,{\sigma }_{32}=0.99$. Then, ${\mathfrak{R}}_0=0.5873,$ ${\displaystyle Q=1/16340590,{\int }_0^{\infty }\left|x-\frac{\Lambda }{\mu }\right|\pi \left(x\right)dx=16.696}$, ${\mathfrak{R}}_0^E=0.908<1$. According to Theorem 2, the disease will die out (see Figure 2).

6. Conclusions and Discussion

This paper investigates a second-order stochastic SEIR model with infectivity in the incubation period and homestead-isolation of the susceptible population. Firstly, we study the existence of a unique ergodic stationary distribution of model (3), which shows that the disease will last for a long time. Secondly, we obtain the corresponding probability density function with respect to the stationary solution to model (3). Finally, the sufficient conditions for disease’s extinction in model (3) are obtained.

However, because of the limitation of our mathematical approaches and the complexity of model (3), we have no ability to obtain a reasonable result concerning the threshold of disease eradication. In addition, we shall investigate the corresponding dynamics of the stochastic model with regime switching. We therefore leave these problems for our future work.