Threshold Dynamics and the Density Function of the Stochastic Coronavirus Epidemic Model

Abstract

Since November 2019, each country in the world has been affected by COVID-19, which has claimed more than four million lives. As an infectious disease, COVID-19 has a stronger transmission power and faster propagation speed. In fact, environmental noise is an inevitable important factor in the real world. This paper mainly gives a new random infectious disease system under infection rate environmental noise. We give the existence and uniqueness of the solution of the system and discuss the ergodic stationary distribution and the extinction conditions of the system. The probability density function of the stochastic system is studied. Some digital simulations are used to demonstrate the probability density function and the extinction of the system.

1. Introduction

So far, infectious diseases have become one of the important factors endangering human health. Medical research shows that there are three outcomes of any infectious disease: the first is that the infectious virus is eliminated by human drugs; the second is that the virus exists only in a small area, such as Ebola, SARS (severe acute respiratory syndrome), and so on; the third is the long-term coexistence of viruses and humans, such as influenza, AIDS (acquired immunodeficiency syndrome), and so on [1].

In order to describe the dynamic behavior of the epidemic, mathematical modeling is considered to be an important tool [2]. According to the occurrence, transmission, and development law of infectious disease in the population, Mathematicians and ecologists have established several epidemic models to study and control various epidemics [3,4,5]. The authors gave the dynamics and stationary distribution of the hepatitis E model. Meanwhile, the authors obtained the optimal control analysis and the Atangana–Baleanu derivative for the dynamical analysis of the hepatitis E model in [6]. Through the qualitative analysis and numerical simulation of the dynamic behavior of the infectious disease model [7], the authors gave an exact expression of the probability density function of the stochastic model SVI (susceptible, vaccinated, infectious) around the unique endemic equilibrium of the deterministic system by solving the corresponding Fokker–Planck equation, which is guaranteed by a new critical value $R_0^s$ in [8], and other models, such as SIR (susceptible, infectious, recovered), SIRS (susceptible, infectious, recovered, susceptible), SEIR (susceptible, exposed, infectious, recovered), etc. [9,10].

Since November 2019, the world has been enveloped in COVID-19 (coronavirus disease 2019). As a contagious virus, COVID-19 is highly infectious [11]. Since April 2021, only half a year, the number of newly confirmed cases in the world has increased by 100 million [12]. At the same time, Griffin B.D., etc. [13], found that the new coronavirus could infect and spread among North American deer rats, which has increased the difficulty in controlling COVID-19. Many countries are working hard to prevent the spread of COVID-19. Mathematicians use simulated mathematical models to predict the future behavior of coronavirus disease transmission in 2019. So far, several typical 2019 coronavirus disease transmission models have been proposed and are being used in some decisions. Recently, people have grasped valuable inferences through mathematical modeling and obtained in-depth understanding of the novel coronavirus (COVID-19) [14]. A. Atangana confirmed the effect of lockdown as a possible adequate measure to help flatten the curve of deaths and infections with the epidemic model as follows: [15]:

$$\left\{\begin{array}{c}dS\left(t\right)=[\Lambda -\frac{\beta S\left(t\right)D\left(t\right)}{N}-(\delta +\mu )S\left(t\right)+\eta R\left(t\right)]dt,\hfill \\ dC\left(t\right)=[(\frac{\beta S\left(t\right)D\left(t\right)}{N}+\delta S\left(t\right)-(\beta +\mu +\pi )C\left(t\right)]dt,\hfill \\ dI\left(t\right)=\left[\delta \right(1-\theta \left)S\right(t)+\pi C(t)-(\tau +\mu +\sigma \left)I\right(t\left)\right]dt,\hfill \\ dR\left(t\right)=\left[\beta C\right(t)+\tau I(t)-(\mu +\eta \left)R\right(t\left)\right]dt,\hfill \\ dD\left(t\right)=\left[\sigma I\right(t\left)\right]dt,\hfill \end{array}\right.$$

(1)

where the parameters are in

Table 1. The definitions of the parameters.

Parameter	Definitions
$S\left(t\right)$	The susceptible class
$I\left(t\right)$	Infected people
$C\left(t\right)$	Carriers (dead corpse)
$R\left(t\right)$	Recovered persons
$D\left(t\right)$	Total number of deaths
$\mu$	Rate of natural death recruitment rate into $S\left(t\right)$
$\theta$	Probability of an $S\left(t\right)$ class to join $C\left(t\right)$ class
$\sigma$	Death rate induced by COVID-19
$\beta$	Recovery rate of $C\left(t\right)$ class
$\delta$	Force of infection of class $S\left(t\right)$
$\tau$	Recovery rate of $I\left(t\right)$ class
$\pi$	Rate at which an $C\left(t\right)$ class is recovered
$\eta$	Rate at which treated persons become $C\left(t\right)$ class

.

However, in the real world, due to the influence of various factors, such as the environment, a random model is constructed by random components with some distribution. Through the addition of some white noise, these distributions may reflect the uncertainty of the input content or random process [16]. Meanwhile, the quarantined measures play a very important role in fighting and preventing the increase in COVID-19. The authors found that the dynamic system with the external source was more reliable than the suspected people travelling, and that the rate of isolation is extremely important for controlling the increase in the cumulative confirmed people of COVID-2019 [17]. The authors in [18] put forward the stochastic coronavirus epidemic model with the parameter disturbance by the natural mortality rate by translating the quarantined factor as follows:

$$\left\{\begin{array}{c}dS\left(t\right)=[\Lambda -\frac{\beta S\left(t\right)I\left(t\right)}{N}-{\mu }_{0}S\left(t\right)]dt+{\eta }_{1}S\left(t\right)d{B}_1\left(t\right),\hfill \\ dI\left(t\right)=[(\frac{\beta S\left(t\right)I\left(t\right)}{N}-({\gamma }_1+{\mu }_1+{\mu }_0)I\left(t\right)+\sigma Q\left(t\right)]dt+{\eta }_{2}I\left(t\right)d{B}_2\left(t\right),\hfill \\ dQ\left(t\right)=[{\gamma }_{1}I\left(t\right)-({\mu }_0+\mu +\sigma )Q\left(t\right)]dt+{\eta }_{3}Q\left(t\right)d{B}_3\left(t\right),\hfill \end{array}\right.$$

(2)

where the definitions of the parameters are in

Table 2. The definitions of the parameters.

Parameter	Definitions
$S\left(t\right)$	The susceptible class
$I\left(t\right)$	Infected people
$Q\left(t\right)$	Quarantined people
N	Total population
$\Lambda$	Capita constant fecundity rate
$\beta$	Infection rate
${\mu }_0$	Infected natural mortality rate
${\mu }_1$	Quarantined natural mortality rate
$\mu$	Disease-related mortality rate
${\gamma }_1$	The constant rate of quarantining infected
$\sigma$	The quarantined rate from infected people
$B_i\left(t\right),i=1,2,3$	Brownian motion
${\eta }_i,i=1,2,3$	The intensity of $B_i\left(t\right)$

and $N=S\left(t\right)+I\left(t\right)+Q\left(t\right)$.

In the system (2), the important role of isolation in COVID-19 is pointed out, and the stable distribution of the model under extinction conditions is obtained. However, in the real world, according to the COVID-19 data in Pakistan [19], the disturbance of the infection rate coefficient plays a very important role in the spread of COVID-19. Meanwhile, vaccination and isolation measures can also affect the infection rate. In the present paper, we give the stochastic coronavirus epidemic model with the stochastic disturbance of the infection rate coefficient. The system is the following:

$$\left\{\begin{array}{c}dS\left(t\right)=[\Lambda -\frac{\beta S\left(t\right)I\left(t\right)}{N}-{\mu }_{0}S\left(t\right)]dt-\frac{\eta S\left(t\right)I\left(t\right)}{N}dB\left(t\right),\hfill \\ dI\left(t\right)=[(\frac{\beta S\left(t\right)I\left(t\right)}{N}-({\gamma }_1+{\mu }_1+{\mu }_0)I\left(t\right)+\sigma Q\left(t\right)]dt+\frac{\eta S\left(t\right)I\left(t\right)}{N}dB\left(t\right),\hfill \\ dQ\left(t\right)=[{\gamma }_{1}I\left(t\right)-({\mu }_0+\mu +\sigma )Q\left(t\right)]dt,\hfill \end{array}\right.$$

(3)

The main composition of the present paper is as follows. The second section gives the basic lemma and basic concepts of this paper. The existence and uniqueness of the global positive solution of the system (3) are obtained in the third section, and the fourth section gives the ergodic stationary distribution of the system (3). In order to better understand the degree of control of the virus, we consider the extinction condition of the system (3) in the fifth section. Meanwhile, in the sixth section, the probability density function of the system (3) is given to understand the trend of the coronavirus viruses in the system (3). In the last section, by numerical simulation, two examples give the extinction, long-term persistence, and the probability density function with the corresponding conditions.

2. Preliminaries

We give some basic conceptions as in [3,5,16,20]. Suppose $(\mathsf{\Omega },\mathcal{F},{\left\{{\mathcal{F}}_t\right\}}_{t\ge 0},\mathbb{P})$ is a complete probability space with a filtration ${\left\{{\mathcal{F}}_t\right\}}_{t\ge 0}$; we define ${\mathbb{R}}_{+}^3=\{x\in {\mathbb{R}}^3:x_i>0\mathrm{for} \mathrm{all}1\le i\le 3\}$ and ${\overline{\mathbb{R}}}_{+}^3=\{x\in {\mathbb{R}}^3:x_i\ge 0\mathrm{for} \mathrm{all}1\le i\le 3\}$. In addition, if $f\left(t\right)$ is an integral function on $t\in [0,\infty )$, we define $f^{\mu }=sup\{f\left(t\right)\mid t\ge 0\},$$f^l=inf\{f\left(t\right)\mid t\ge 0\}$. In the following, we give the It${\widehat{o}}^{\prime }$s formula.

Lemma 1 

([3]). Let $x\left(t\right)$ be an It${\widehat{o}}^{\prime }$ process with the stochastic differential

$$dx\left(t\right)=f\left(t\right)dt+g\left(t\right)d{B}_t,fort\ge t_0,$$

(4)

where $f\in L^1({\mathbb{R}}_{+},\mathbb{R})$ and $g\in L^1({\mathbb{R}}_{+},\mathbb{R})$. Let $V\in C^{2,1}(\mathbb{R}\times {\mathbb{R}}_{+},\mathbb{R})$. Then, $V\left(x\right(t),t)$ is again an It${\widehat{o}}^{\prime }$ process with the stochastic differential given by

$$dV(x\left(t\right),t)=[V_t(x\left(t\right),t)+V_x(x\left(t\right),t)f\left(t\right)+\frac{1}{2}{V}_{xx}(x\left(t\right),t)g^2\left(t\right)]dt+V_x(x\left(t\right),t)g\left(t\right)d{B}_{t}a.s.$$

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

$$dx\left(t\right)=f(x\left(t\right),t)dt+g(x\left(t\right),t)dB\left(t\right),fort\ge t_0$$

(5)

with initial value $x\left(t_0\right)=x_0\in {\mathbb{R}}^3$, where $B\left(t\right)$ denotes three-dimensional standard Brownian motion defined on the above probability space $(\mathsf{\Omega },\Gamma ,{\left\{{\Gamma }_t\right\}}_{t\ge 0},\mathbb{P})$. Define the differential operator L by Mao [3] as

$$L=\frac{\partial }{\partial t}+\Sigma f_i(x,t)\frac{\partial }{\partial x_i}+\frac{1}{2}\Sigma {\left[g^T(x,t)g(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}}^3\times {\overline{\mathbb{R}}}_{+};{\overline{\mathbb{R}}}_{+})$, where ${\overline{\mathbb{R}}}_{+}=\{x\in \mathbb{R}:x\ge 0\}$, then

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

where $V_t=\frac{\partial V}{\partial t},V_x=(\frac{\partial V}{\partial x_1},\cdots ,\frac{\partial V}{\partial x_3})$ and $V_{xx}={\left(\frac{{\partial }^{2}V}{\partial x_i\partial x_j}\right)}_{3\times 3}$.

By Lemma 1, we obtain

$$dV(x\left(t\right),t)=LV(x\left(t\right),t)dt+V_x(x\left(t\right),t)g(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}^3{g}_r^i\left(x\right)g_r^j\left(x\right),1\le i,j\le 3.$$

3. Existence and Uniqueness of the Global Positive Solution

The problem where the solution is global and positive is important in studying the dynamical behavior of the system (3). The coefficients of the system (3) are not the linear growth, and the solutions of the system (3) may explode at a fixed time. The main theorem is as follows.

Theorem 2. 

There is a unique positive solution $\left(S\right(t),I(t),Q(t\left)\right)$ of system (3) on $t\ge 0$ by the initial value $(S\left(0\right),I\left(0\right),Q\left(0\right))\in {\mathbb{R}}_{+}^3$, and the solution $(S\left(t\right),I\left(t\right),Q\left(t\right))\in {\mathbb{R}}_{+}^3$ for all $t\ge 0$ almost surely (a.s.).

Proof. 

Based on [5], we obtain the fact that there is a unique solution $\left(x\right(t),y(t),z(t),w(t\left)\right)$ on $[0,{\tau }_0)$ for the reason that the coefficients of the system (3) are the locally Lipschitz continuous, where ${\tau }_0$ is an explosion time. We can obtain the fact that the local solution is global when ${\tau }_0=\infty$ a.s. By the definitions in [3], we define a fundamental $C^2$-function $U:{\mathbb{R}}_{+}^3\to {\overline{\mathbb{R}}}_{+}$, which is

$$\begin{array}{c}U(S,I,Q)=\left(S\right(t)-1-lnS(t\left)\right)+a\left(I\right(t)-1-\beta lnI(t\left)\right)+b\left(Q\right(t)-1-lnQ(t\left)\right),\hfill \end{array}$$

(6)

where $a,b$ 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 [3], we obtain

$$\begin{array}{c}dU(S,I,Q)=LUdt+\frac{\eta \left(S\right(t)-1)I\left(t\right)}{N^2}dB\left(t\right)-\frac{\eta \left(I\right(t)-1)S\left(t\right)}{N^2}dB\left(t\right),\hfill \end{array}$$

(7)

where

$$\begin{array}{cc}LU\hfill & =(1-\frac{1}{S})(\Lambda -\frac{\beta SI}{N}-{\mu }_{0}S)-\frac{{\eta }^2}{2}{\left(\frac{I}{N}\right)}^2+(1-\frac{1}{I})[(\frac{\beta SI}{N}-({\gamma }_1+{\mu }_1+{\mu }_0)I+\sigma Q]+\frac{{\eta }^2}{2}{\left(\frac{S}{N}\right)}^2\hfill \\ & +(1-\frac{1}{Q})[{\gamma }_{1}I-({\mu }_0+\mu +\sigma )Q]\hfill \\ & =\Lambda -\frac{\beta SI}{N}-{\mu }_{0}S-\frac{\Lambda }{S}+\frac{\beta I}{N}+{\mu }_0-({\mu }_0+{\mu }_1+{\gamma }_1)I+\frac{\beta SI}{N}+\sigma Q-\frac{\beta S}{N}+({\mu }_0+{\mu }_1+{\gamma }_1)-\frac{\sigma Q}{I}+{\gamma }_{1}I\hfill \\ & -({\mu }_0+\mu +\sigma )Q-\frac{{\gamma }_{1}I}{Q}+({\mu }_0+\mu +\sigma )-\frac{{\eta }^2}{2}{\left(\frac{I}{N}\right)}^2+\frac{{\eta }^2}{2}{\left(\frac{S}{N}\right)}^2\hfill \\ & \le \Lambda -{\mu }_{0}S-\frac{\Lambda }{S}+\beta +{\mu }_0-({\mu }_0+{\mu }_1+{\gamma }_1)I+\sigma Q+({\mu }_0+{\mu }_1+{\gamma }_1)-({\mu }_0+\mu +\sigma )Q\hfill \\ & -\frac{\sigma Q}{I}+{\gamma }_{1}I-\frac{{\gamma }_{1}I}{Q}+({\mu }_0+\mu +\sigma )+\frac{{\eta }^2}{2}\hfill \\ & \le \Lambda -{\mu }_{0}S-\frac{\Lambda }{S}+\beta +{\mu }_0+({\mu }_0+{\mu }_1+{\gamma }_1)-\frac{\sigma Q}{I}-\frac{{\gamma }_{1}I}{Q}+({\mu }_0+\mu +\sigma )+\frac{{\eta }^2}{2}\hfill \\ & \le \Lambda +\beta +3{\mu }_0+{\mu }_1+{\gamma }_1+\mu +\sigma -2\sqrt{{\mu }_0\Lambda }-2\sqrt{\sigma {\gamma }_1}+\frac{{\eta }^2}{2}\hfill \\ & ={(\sqrt{\Lambda }-\sqrt{{\mu }_0})}^2+{(\sqrt{\sigma }-\sqrt{{\gamma }_1})}^2+\beta +2{\mu }_0+{\mu }_1+\mu +\frac{{\eta }^2}{2},\hfill \end{array}$$

(8)

Then, we can obtain

$$\begin{array}{ccc}LU\le {(\sqrt{\Lambda }-\sqrt{{\mu }_0})}^2+{(\sqrt{\sigma }-\sqrt{{\gamma }_1})}^2+\beta +2{\mu }_0+{\mu }_1+\mu +\frac{{\eta }^2}{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 [5]. Hence, we omit it here. □

4. Ergodic Stationary Distribution of the Stochastic Coronavirus Epidemic Model

In this section, the existence of ergodic stationary components of the system (3) is given. Firstly, we define $R_0^{*}$ as a stochastic reproductive ratio of the system (3), such as

$$R_0^{*}=\frac{{\mu }_0\beta \sigma {\gamma }_1}{({\mu }_0+\frac{1}{2}{\eta }^2){({\gamma }_1+{\mu }_1+{\mu }_0+\frac{1}{2}{\eta }^2)}^2({\mu }_0+\mu +\sigma )},$$

which is equal to $\frac{\beta \sigma {\gamma }_1}{{({\gamma }_1+{\mu }_1+{\mu }_0)}^2({\mu }_0+\mu +\sigma )}$ when $\eta =0$ [3]. The following is a known result.

Lemma 3 

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

• (B.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$;
• (B.2) there exists a nonnegative $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 μ.

Theorem 4. 

When $R_0^{*}>1$, for the solution $\left(S\right(t),I(t),Q(t\left)\right)$ of the system (3), there exists an ergodic unique stationary distribution.

Proof. 

We construct a $C^2$-function $\tilde{V}:{\mathbb{R}}_{+}^3\to \mathbb{R}$ as follows:

$$\tilde{V}=N\left(t\right)-c_{1}lnS\left(t\right)-c_{2}lnI\left(t\right)-c_{3}lnQ\left(t\right).$$

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

$$\begin{array}{cc}L\tilde{V}\hfill & =(\Lambda -{\mu }_{0}N-{\mu }_{1}I-\mu Q)+c_1[-\frac{\Lambda }{S}+\frac{\beta I}{N}+{\mu }_0+\frac{1}{2}{\left(\eta IN\right)}^2]+c_2[-\frac{\beta S}{N}+({\gamma }_1+{\mu }_1+{\mu }_0)-\sigma \frac{Q}{I}+\frac{1}{2}{\left(\frac{\eta S}{N}\right)}^2]\hfill \\ & +c_3[-{\gamma }_1\frac{I}{Q}+({\mu }_0+\mu +\sigma )]\hfill \\ & =c_1\frac{\beta I}{N}+[-{\mu }_{0}N-c_1\frac{\Lambda }{S}-c_2\frac{\beta S}{N}-c_2\sigma \frac{Q}{I}-c_3{\gamma }_1\frac{I}{Q}]+c_1({\mu }_0+\frac{1}{2}{\eta }^2)+c_2({\gamma }_1+{\mu }_1+{\mu }_0+\frac{1}{2}{\eta }^2)\hfill \\ & +c_3({\mu }_0+\mu +\sigma )+\Lambda -{\mu }_{1}I-\mu Q\hfill \\ & \le c_1\frac{\beta I}{N}-5{({\mu }_0{c}_1{c}_2^2{c}_3\beta \sigma \Lambda {\gamma }_1)}^{\frac{1}{5}}+c_1({\mu }_0+\frac{1}{2}{\eta }^2)+c_2({\gamma }_1+{\mu }_1+{\mu }_0+\frac{1}{2}{\eta }^2)+c_3({\mu }_0+\mu +\sigma )+2\Lambda \hfill \\ & \le c_1\frac{\beta I}{N}-5\Lambda [{\left(R_0^{*}\right)}^{\frac{1}{5}}-1],\hfill \end{array}$$

(10)

and $R_0^{*}=\frac{{\mu }_0\beta \sigma {\gamma }_1}{({\mu }_0+\frac{1}{2}{\eta }^2){({\gamma }_1+{\mu }_1+{\mu }_0+\frac{1}{2}{\eta }^2)}^2({\mu }_0+\mu +\sigma )}$. We choose $c_1=\frac{\Lambda }{{\mu }_0+\frac{1}{2}{\eta }^2}$, $c_2=\frac{\Lambda }{{\gamma }_1+{\mu }_1+{\mu }_0+\frac{1}{2}{\eta }^2}$, $c_3=\frac{\Lambda }{{\mu }_0+\mu +\sigma }.$

When $R_0^{*}>1$, we suppose

$$\overline{V}=M\tilde{V}-lnS\left(t\right)-lnI\left(t\right)-lnQ\left(t\right)+N\left(t\right),$$

and $V:=\overline{V}(S,I,Q)-\overline{V}(S_0,I_0,Q_0).$

Applying It${\widehat{o}}^{\prime }$s formula to V, we obtain

$$\begin{array}{cc}LV\hfill & =ML\tilde{V}-LlnS\left(t\right)-LlnI\left(t\right)-LlnQ\left(t\right)+LN\left(t\right)\hfill \\ & =-5\Lambda M[{\left(R_0^{*}\right)}^{\frac{1}{5}}-1]+c_{1}M\frac{\beta I}{N}+(-\frac{\Lambda }{S}+\frac{\beta I}{N}+{\mu }_0+\frac{1}{2}{\eta }^2+[-\frac{\beta S}{N}+({\gamma }_1+{\mu }_1+{\mu }_0)-\sigma \frac{Q}{I}-\frac{1}{2}{\eta }^2]\hfill \\ & +[-{\gamma }_1\frac{I}{Q}+({\mu }_0+\mu +\sigma )]+\Lambda -{\mu }_{0}N-{\mu }_{1}I-\mu Q\hfill \\ & \le -5\Lambda M[{\left(R_0^{*}\right)}^{\frac{1}{5}}-1]+(c_{1}M+1)\beta -\frac{\Lambda }{S}-\sigma \frac{Q}{I}-{\gamma }_1\frac{I}{Q}-\frac{\beta S}{N}+3{\mu }_0+{\gamma }_1+{\mu }_1+\mu +\sigma +\Lambda \hfill \\ & -{\mu }_{0}N-{\mu }_{1}I-\mu Q\hfill \\ & \le (c_{1}M+1)\beta +3{\mu }_0+{\gamma }_1+{\mu }_1+\mu +\sigma +\Lambda -\frac{\Lambda }{S}-\sigma \frac{Q}{I}-{\gamma }_1\frac{I}{Q}-{\mu }_{0}S-({\mu }_1+{\mu }_0)I-(\mu +{\mu }_0)Q.\hfill \end{array}$$

(11)

Define

$$f_1\left(S\right)=(c_{1}M+1)\beta +3{\mu }_0+{\gamma }_1+{\mu }_1+\mu +\sigma +\Lambda -\frac{\Lambda }{S}-{\mu }_{0}S,$$

$$f_2\left(I\right)=-\sigma \frac{Q}{I}-({\mu }_1+{\mu }_0)I,$$

$$f_3\left(Q\right)=-{\gamma }_1\frac{I}{Q}-(\mu +{\mu }_0)Q.$$

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

$$D_1=\{(S,I,Q)\in {\mathbb{R}}_{+}^3:0<S<\epsilon \};D_2=\{(S,I,Q)\in {\mathbb{R}}_{+}^3:S>\frac{1}{\epsilon }\};$$

$$D_3=\{(S,I,Q)\in {\mathbb{R}}_{+}^3:0<I<\epsilon \};D_4=\{(S,I,Q)\in {\mathbb{R}}_{+}^3:I>\frac{1}{\epsilon }\};$$

$$D_5=\{(S,I,Q)\in {\mathbb{R}}_{+}^3:Q<{\epsilon }^2,I>\epsilon \};D_6=\{(S,I,Q)\in {\mathbb{R}}_{+}^3:Q>\frac{1}{\epsilon }\};$$

Clearly, $D_{\epsilon }={\bigcup }_{j=1}^6{D}_j$. In the following text, we will show that $LV(S,I,Q)\le -1$ on ${\mathbb{R}}_{+}^3\backslash D_{\epsilon }$.

Case 1. If $(S,I,Q)\in D_1\bigcup D_2$, one can choose

$$M<-\frac{{\mu }_0+\frac{1}{2}{\eta }^2}{\Lambda },$$

and

$$LV(S,I,Q)\le 3{\mu }_0+{\gamma }_1+{\mu }_1+\mu +\sigma +\Lambda -\frac{\Lambda }{S}-{\mu }_{0}S\le -2;$$

Case 2. If $(S,I,Q)\in D_3\bigcup D_4$,

$$LV(S,I,Q)\le (c_{1}M+1)\beta +3{\mu }_0+{\gamma }_1+{\mu }_1+\mu +\sigma +\Lambda +f_2\left(I\right)\le -2;$$

Case 3. If $(S,I,Q)\in D_5\bigcup D_6$,

$$LV(S,I,Q)\le (c_{1}M+1)\beta +3{\mu }_0+{\gamma }_1+{\mu }_1+\mu +\sigma +\Lambda +f_3\left(Q\right)\le -2;$$

therefore, for all $(S,I,Q))\in {\mathbb{R}}_{+}^3\backslash D_{\epsilon }$, $LV(S,I,Q))\le -1,$ which indicates that assumption $(B.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 Stochastic Coronavirus Epidemic Model

It is a very important topic to consider the dynamic behavior of the epidemic virus to obtain the conditions for the virus to be eliminated in a long time. We mainly discuss the extinction conditions of the system (3). According to the results in [10], we can obtain the following lemma.

Lemma 5. 

For any initial value, the solution of 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.$$

(12)

$$\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..$$

(13)

Moreover,

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

(14)

Theorem 6. 

Let $\left(S\right(t),I(t),Q(t\left)\right)$ be the solution of system (3) with any initial value $\left(S\right(0),I(0),$$Q\left(0\right))\in {\mathbb{R}}_{+}^3$. If $R_0^S<1$, then the solution $\left(S\right(t),I(t),Q(t\left)\right)$ of system (3) satisfies

$$\underset{t\to \infty }{lim\; sup}\frac{lnI\left(t\right)}{t}\le \frac{1}{{\mu }_0+{\gamma }_1+{\mu }_1+\mu +\sigma }(R_0^S-1)<0a.s.,$$

where $R_0^S=\frac{\beta }{2{\eta }^2}({\mu }_0+{\gamma }_1+{\mu }_1+\mu +\sigma )$. Namely, the disease will be eradicated in the long term.

Proof. 

Applying It${\widehat{o}}^{\prime }$s formula to $lnI\left(t\right)$, we obtain

$$\begin{array}{cc}dlnI\left(t\right)\hfill & =\frac{dI\left(t\right)}{I\left(t\right)}=(\frac{\beta S}{N}-({\gamma }_1+{\mu }_1+{\mu }_0)+\sigma \frac{Q}{I})dt-\frac{1}{2}{\left(\frac{\eta S}{N}\right)}^{2}dt-\frac{\eta S}{N}dB\left(t\right)\hfill \\ & \le [\frac{\beta S}{N}-\frac{1}{2}{\left(\frac{\eta S}{N}\right)}^2-({\gamma }_1+{\mu }_1+{\mu }_0)+\sigma ]dt-\frac{\eta S}{N}dB\left(t\right)\hfill \\ & \le \{-\frac{{\eta }^2}{2}[{\left(\frac{S}{N}\right)}^2-\frac{2\beta }{{\eta }^2}\frac{S}{N}+{\left(\frac{\beta }{{\eta }^2}\right)}^2-{\left(\frac{\beta }{{\eta }^2}\right)}^2]-({\gamma }_1+{\mu }_1+{\mu }_0)+\sigma \}dt-\frac{\eta S}{N}dB\left(t\right)\hfill \\ & \le [\frac{{\beta }^2}{2{\eta }^2}-({\gamma }_1+{\mu }_1+{\mu }_0-\sigma )]dt-\frac{\eta S}{N}dB\left(t\right).\hfill \end{array}$$

(15)

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

$$lnI\left(t\right)-lnI\left(0\right)\le {\int }_0^t[\frac{{\beta }^2}{2{\eta }^2}-({\gamma }_1+{\mu }_1+{\mu }_0-\sigma )]ds-{\int }_0^t\frac{\eta S}{N}dB\left(t\right).$$

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

$$\underset{t\to 0}{lim}\frac{1}{t}{\int }_0^t\frac{\eta S}{N}dB\left(t\right)=0a.s..$$

and we can obtain

$$\begin{array}{cc}\underset{t\to \infty }{lim\; sup}\frac{lnI\left(t\right)}{t}\hfill & \le {\int }_0^t[\frac{{\beta }^2}{{\eta }^4}-({\gamma }_1+{\mu }_1+{\mu }_0-\sigma )]ds\hfill \\ & \le \frac{{\beta }^2}{{\eta }^4}-({\gamma }_1+{\mu }_1+{\mu }_0-\sigma )\hfill \\ & =({\gamma }_1+{\mu }_1+{\mu }_0-\sigma )(\frac{{\beta }^2}{{\eta }^4({\gamma }_1+{\mu }_1+{\mu }_0-\sigma )}-1)\hfill \\ & <0a.s..\hfill \end{array}$$

(16)

We choose $R_0^S=\frac{\beta }{2{\eta }^2}({\mu }_0+{\gamma }_1+{\mu }_1+\mu +\sigma )<1,$ which is equal to $\eta >{\left(\frac{{\eta }^2}{{\mu }_0+{\gamma }_1+{\mu }_1+\mu +\sigma }\right)}^{\frac{1}{4}}.$

Therefore, the above indicates that

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

and we can obtain the fact that the viral will be eradicated, which completes the proof. □

6. The Probability Density Function of the Stochastic Coronavirus Epidemic Model

Let $N=S+I+Q$; we can transfer system (3) into the following system:

$$\left\{\begin{array}{c}dS\left(t\right)=[\Lambda -\frac{\beta S\left(t\right)\left(N\right(t)-S(t)-Q(t\left)\right)}{N}-{\mu }_{0}S\left(t\right)]dt-\frac{\eta S\left(t\right)\left(N\right(t)-S(t)-Q(t\left)\right)}{N}dB\left(t\right)\hfill \\ dQ\left(t\right)=({\gamma }_{1}N\left(t\right)-{\gamma }_{1}S\left(t\right)-({\gamma }_1+{\mu }_0+\mu +\sigma )Q\left(t\right))dt\hfill \\ dN\left(t\right)=[\Lambda -({\mu }_1+{\mu }_0)N\left(t\right)+{\mu }_{1}S-({\mu }_1-\mu )Q\left(t\right)]dt.\hfill \end{array}\right.$$

(17)

Hence, we can consider the probability density function of system (17) in place of system (3).

Theorem 7. 

We consider the condition that $R_0^{*}>1$, for any initial value $\left(S\right(0),Q(0),N(0\left)\right)$$\in {\mathbb{R}}_{+}^3$; then, the solution $\left(S\right(t),Q(t),N(t\left)\right)$ of system (3) with a weak kernel will have a normal probability density function $\Phi \left(S\right(t),Q(t),N(t\left)\right)$ around $(S^{*},Q^{*},N^{*})$, which is given by

$$\Phi (S\left(t\right),Q\left(t\right),N\left(t\right))={\left(2\pi \right)}^{-\frac{3}{2}}\mid \Sigma {\mid }^{-\frac{1}{2}(S\left(t\right),Q\left(t\right),N\left(t\right)){\Sigma }^{-1}{(S\left(t\right),Q\left(t\right),N\left(t\right))}^T},$$

where Σ is a positive definite matrix and satisfies

$$\Sigma =\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}{2a_3(a_1{a}_2-a_3)}\end{array}\right),$$

where $a_1={\gamma }_1+2{\mu }_0+\mu +\sigma +\frac{\beta (N^{*}-2S^{*}-Q^{*})}{N^{*}}$, $a_2=({\mu }_0-{\mu }_1)({\gamma }_1+2{\mu }_0+\mu +\sigma )+$$\frac{\beta (N^{*}-2S^{*}-Q^{*})}{N^{*}}(2{\mu }_0+\mu +\sigma )+\frac{{\gamma }_1\beta (N^{*}-S^{*})}{N^{*}},a_3=\frac{\beta S^{*}(N^{*}-S^{*}-Q^{*})}{{\left(N^{*}\right)}^2}+({\gamma }_1\mu +({\mu }_0-{\mu }_1)({\gamma }_1+{\mu }_0+\mu +\sigma ))({\mu }_0+\frac{\beta (N^{*}-2S^{*}-Q^{*})}{N^{*}})+{\mu }_1(2{\gamma }_1+{\mu }_0+\mu +\sigma )$.

Proof. 

Firstly, we can obtain the linear system of system (17) at point $(y_1,y_2,y_3)=(S^{*},Q^{*},N^{*}).$

$$\left\{\begin{array}{c}d{y}_1=(b_{11}{y}_1\left(t\right)+b_{12}{y}_2\left(t\right)+b_{13}{y}_3\left(t\right))dt-\frac{\eta S^{*}(N^{*}-S^{*}-Q^{*})}{N^{*}}dB\left(t\right)\hfill \\ d{y}_2=(b_{21}{y}_1\left(t\right)+b_{22}{y}_2\left(t\right)+b_{23}{y}_3\left(t\right))dt\hfill \\ d{y}_3=(b_{31}{y}_1\left(t\right)+b_{32}{y}_2\left(t\right)+b_{33}{y}_3\left(t\right))dt,\hfill \end{array}\right.$$

(18)

where

$$A=\left(\begin{array}{ccc}{b}_{11}& b_{12}& b_{13}\\ b_{21}& b_{22}& b_{23}\\ b_{31}& b_{32}& b_{33}\end{array}\right)=\left(\begin{array}{ccc}-({\mu }_0+\frac{\beta (N^{*}-2S^{*}-Q^{*})}{N^{*}})& \frac{\beta S^{*}}{N^{*}}& -\frac{\beta S^{*}(S^{*}+Q^{*})}{{\left(N^{*}\right)}^2}\\ -{\gamma }_1& -({\gamma }_1+{\mu }_0+\mu +\sigma )& {\gamma }_1\\ {\mu }_1& \mu -{\mu }_1& -({\mu }_1+{\mu }_0)\end{array}\right),$$

$$B=\left(\begin{array}{ccc}\frac{\eta S^{*}(N^{*}-S^{*}-Q^{*})}{N^{*}}& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right),$$

Let $Y={(y_1,y_2,y_3)}^T$, $G=diag(\frac{\eta S^{*}(N^{*}-S^{*}-Q^{*})}{N^{*}},0,0)$; then,

$$dY=AYdt+GdB\left(t\right).$$

By the Roozen [21], we can obtain the fact that the density function $\Phi \left(Y\right)=\Phi (y_1,y_2,y_3)$ of system (18) nearby the origin point can approximate to the Fokker–Plank equation as follows:

$$\begin{array}{cc}-\frac{{\eta }^2}{2}\frac{{\partial }^2\Phi \left(t\right)}{\partial y_1^2}\hfill & +\frac{\partial }{\partial y_1}[(b_{11}{y}_1\left(t\right)+b_{12}{y}_2\left(t\right)+b_{13}{y}_3\left(t\right))\Phi \left(t\right)]\hfill \\ & +\frac{\partial }{\partial y_2}[(b_{21}{y}_1\left(t\right)+b_{22}{y}_2\left(t\right)+b_{23}{y}_3\left(t\right))\Phi \left(t\right)]\hfill \\ & +\frac{\partial }{\partial y_3}[(b_{31}{y}_1\left(t\right)+b_{32}{y}_2\left(t\right)+b_{33}{y}_3\left(t\right))\Phi \left(t\right)]\hfill \\ & =0.\hfill \end{array}$$

(19)

By Gaussian distribution,

$$\Phi \left(Y\right)=Cexp\{-\frac{1}{2}(Y-Y^{*})P{(Y-Y^{*})}^T\},$$

(20)

where P is a real symmetric matrix that satisfies

$$P{G}^{2}P+A^{T}P+PA=0.$$

(21)

Let $P^{-1}=\Sigma$; then, we have

$$G^2+A\Sigma +\Sigma A^T=0.$$

(22)

We know there exists a matrix

$$\tilde{M}=\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& \frac{{\mu }_1}{{\gamma }_1}& 1\end{array}\right),$$

satisfying

$$\tilde{M}A{\tilde{M}}^{-1}=\left(\begin{array}{ccc}-a_{11}& a_{12}& a_{13}\\ a_{21}& -a_{22}& a_{23}\\ a_{31}& a_{32}& -a_{33}\end{array}\right)=\left(\begin{array}{ccc}-({\mu }_0+\frac{\beta (N^{*}-2S^{*}-Q^{*})}{N^{*}})& \frac{\beta S^{*}}{N^{*}}& -\frac{\beta S^{*}(S^{*}+Q^{*})}{{\left(N^{*}\right)}^2}\\ -{\gamma }_1& -({\gamma }_1+{\mu }_0+\mu +\sigma )& {\gamma }_1\\ 0& \frac{-{\mu }_1(2{\gamma }_1+{\mu }_0+\mu +\sigma )+{\gamma }_1\mu }{{\gamma }_1}& -{\mu }_0\end{array}\right).$$

Hence, we have the characteristic polynomials of A as

$${\phi }_A\left(\lambda \right)={\lambda }^3+a_1{\lambda }^2+a_2\lambda +a_3.$$

Denote

$$dY=d\left(\begin{array}{c}{y}_1\\ y_2\\ y_3\end{array}\right)=\left(\begin{array}{ccc}-a_1& -a_2& -a_3\\ 1& 0& 0\\ 0& 1& 0\end{array}\right)\left(\begin{array}{c}{y}_1\\ y_2\\ y_3\end{array}\right)dt,$$

where $a_1=a_{11}+a_{22}+a_{33}={\gamma }_1+2{\mu }_0+\mu +\sigma +\frac{\beta (N^{*}-2S^{*}-Q^{*})}{N^{*}}>0,a_2=a_{11}{a}_{22}+a_{11}{a}_{33}+a_{22}{a}_{33}+a_{23}{a}_{32}+a_{12}{a}_{21}=({\mu }_0-{\mu }_1)({\gamma }_1+2{\mu }_0+\mu +\sigma )+\frac{\beta (N^{*}-2S^{*}-Q^{*})}{N^{*}}(2{\mu }_0+\mu +\sigma )+\frac{{\gamma }_1\beta (N^{*}-S^{*})}{N^{*}}>0,a_3=a_{11}{a}_{22}{a}_{33}+a_{11}{a}_{23}{a}_{32}+a_{12}{a}_{21}{a}_{33}+a_{13}{a}_{21}{a}_{32}=\frac{\beta S^{*}(N^{*}-S^{*}-Q^{*})}{{\left(N^{*}\right)}^2}+({\gamma }_1\mu +({\mu }_0-{\mu }_1)({\gamma }_1+{\mu }_0+\mu +\sigma ))({\mu }_0+\frac{\beta (N^{*}-2S^{*}-Q^{*})}{N^{*}})+{\mu }_1(2{\gamma }_1+{\mu }_0+\mu +\sigma )>0$. We can easily obtain $a_1{a}_2-a_3>0.$ Therefore, by some transformation, the standard $R_1$ matrix of A is unique. By the same method of the Lemma 3 in [22], we can obtain a positive matrix.

$$\Sigma =\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}{2a_3(a_1{a}_2-a_3)}\end{array}\right).$$

Hence, $\Sigma$ is a positive definite, and we complete the proof. □

7. Examples and Numerical Simulations

In this section, we give the numerical simulation of system (3) by using the discrete equation with the same method as [10]. The equations are as follows:

$$\left\{\begin{array}{c}S(k+1)=S\left(k\right)+[\Lambda -\frac{\beta S\left(t\right)I\left(t\right)}{N\left(t\right)}-{\mu }_{0}S\left(t\right)]▵t-\eta \frac{S\left(t\right)I\left(t\right)}{N\left(t\right)}\sqrt{▵t}{\xi }_k-\frac{{\eta }^2}{2}\frac{S\left(t\right)I\left(t\right)}{N\left(t\right)}▵t({\xi }_k^2-1),\hfill \\ I(k+1)=I\left(k\right)+[(\frac{\beta S\left(t\right)I\left(t\right)}{N}-({\gamma }_1+{\mu }_1+{\mu }_0)I\left(t\right)+\sigma Q\left(t\right)]▵t+\eta \frac{S\left(t\right)I\left(t\right)}{N\left(t\right)}\sqrt{▵t}{\xi }_k+\frac{{\eta }^2}{2}\frac{S\left(t\right)I\left(t\right)}{N\left(t\right)}▵t({\xi }_k^2-1),\hfill \\ Q(k+1)=Q\left(k\right)+[{\gamma }_{1}I\left(t\right)-({\mu }_0+\mu +\sigma )Q\left(t\right)]▵t,\hfill \end{array}\right.$$

(23)

where the time increment $▵t>0$, and ${\xi }_k$ is a the Gaussian random variable $(k=0,1,2,\cdots n).$

Example 1. 

Here, in system (3), we use the environmental noise parameter as $\eta =0.1$. In addition, following the biological feasibility result, the values of the parameters are as shown in

Table 3. Parameters value.

Notation	Value	References
$\Lambda$	0.028	[18]
$\beta$	0.2	Estimated
${\mu }_0$	0.011	[18]
${\mu }_1$	0.2	Estimated
$\gamma$	0.06	[18]
$\sigma$	0.3	Estimated
$\mu$	0.5	[18]

.

In addition, we can choose the following real data $S\left(0\right)=\mathrm{355,250},I\left(0\right)=1453,Q\left(0\right)=\mathrm{51,343},$ in Pakistan on 7 October 2021 [19]. Then, $R_0^{*}=3.208>1,$ where $R_0^{*}$ is defined in Section 4. By the results of Theorem 4, we can find that system (3) will persist for a long time by a distribution $\mu (·)$. The numerical simulations (Figure 1) confirm this.

Example 2. 

Here, in system (3), we use the environmental noise parameter $\eta =0.1$. In addition, following the biological feasibility result, the values of the parameters are as shown in

Table 4. Parameters value.

Notation	Value	References
$\Lambda$	0.5	[15]
$\beta$	0.6	Estimated
${\mu }_0$	0.2	[15]
${\mu }_1$	0.2	Estimated
$\gamma$	0.3	[15]
$\sigma$	0.1	Estimated
$\mu$	0.2	[15]

.

In addition, we can choose the following real data $S\left(0\right)=\mathrm{355,250},I\left(0\right)=1453,Q\left(0\right)=\mathrm{51,343},$ in Pakistan on 7 October 2021 [19]. Then, $R_0^{*}=0.0769<1,$ where $R_0^{*}$ is defined in Section 4. We can find that system (3) will be extinct in a long time. The numerical simulations (Figure 2) confirm this.

In fact, using the statistical data of Pakistan for September to December 2021, it can be seen from the figure (Figure 3, Figure 4 and Figure 5) that the control of the isolation number will affect the disturbance of the infection rate and control the increase in the infection number. At the same time, when the infection rate is disturbed by other factors, such as vaccine injection, the number of deaths decreases with the decrease in the number of infections. This is basically consistent with the research results of system (3) in this paper.

8. Discussion

So far, the COVID-19 coronavirus disease is still one of the most serious diseases in the world. Until today, there is no appropriate treatment. At the same time, due to the strong transmission of the virus, with the existence of many uncertain factors (human activities, animal activities, express delivery, etc.), it also contains a lot of randomness. With the help of stochastic theory, we developed a model for the new 2019 coronavirus disease, and considered studying the transmission characteristics of the disease and understanding its transmission dynamics in the change in population and environment. The important role of isolation measures in controlling transmission is introduced. By disturbing the infection coefficient, the existence and positivity of Lyapunov function theory are studied. In this paper, in order to further discuss the extinction and stable distribution, we gave a new random infectious disease system under infection rate environmental noise. We give the existence and uniqueness of the solution of the system and discuss the ergodic stationary distribution and the extinction conditions of the system. The probability density function of the stochastic system is studied. Some digital simulations are used to demonstrate the probability density function and the extinction of system (3). Through numerical simulation, we analyzed the above results and drew a conclusion with the support of graphics. This work shows that random analysis is a better method used to study the dynamics of infectious diseases, especially the new 2019 coronavirus disease.

We know that the control of infectious diseases needs to consider a variety of random disturbance factors, and we will consider system (3) in more general random phenomena and their persistence and extinction properties in our future research.