A Stochastic Mathematical Model for Understanding the COVID-19 Infection Using Real Data

Abstract

Natural symmetry exists in several phenomena in physics, chemistry, and biology. Incorporating these symmetries in the differential equations used to characterize these processes is thus a valid modeling assumption. The present study investigates COVID-19 infection through the stochastic model. We consider the real infection data of COVID-19 in Saudi Arabia and present its detailed mathematical results. We first present the existence and uniqueness of the deterministic model and later study the dynamical properties of the deterministic model and determine the global asymptotic stability of the system for ${\mathcal{R}}_0\le 1$. We then study the dynamic properties of the stochastic model and present its global unique solution for the model. We further study the extinction of the stochastic model. Further, we use the nonlinear least-square fitting technique to fit the data to the model for the deterministic and stochastic case and the estimated basic reproduction number is ${\mathcal{R}}_0\approx 1.1367$. We show that the stochastic model provides a good fitting to the real data. We use the numerical approach to solve the stochastic system by presenting the results graphically. The sensitive parameters that significantly impact the model dynamics and reduce the number of infected cases in the future are shown graphically.

1. Introduction

The study of differential equations with symmetry (also known as equivariant differential equations) entails studying the structure of differential equation solutions that obey the limitations imposed by the symmetry group. Symmetries can be spatial or temporal in nature. Aside from their application in modeling, differential equations with symmetry have fascinating dynamical properties that are not seen in non-symmetric systems, such as durable heteroclinic cycles. In ordinary differential equations (including Hamiltonian systems) and delay-differential equations, I have looked at symmetry’s importance in bifurcations from equilibria and periodic orbits. Mathematical models to study infectious diseases, in particular COVID-19, are gaining attention from researchers day by day, see [1,2,3]. The COVID-19 infection which brings negative effects by increasing infection and deaths is still active in many countries of the world. As per the follow-up of the World Health Organization (WHO) instructions, many countries are overcoming this infection while some of them are still facing it. Saudi Arabia also suffered from this infection but their policies and restrictions on their citizens to follow the WHO suggestions have succeeded in curtailing the COVID-19 infection. At present, the total number of coronavirus cases since the disease onset in Saudi Arabia is 819,083, with 9383 deaths. The number of individuals recovered so far is 805,670. The percentage of recovery is 99% while the death percentage is 1%. The proper management from the government and the response of citizens in following the guidelines of the World Health Organization can minimize future cases.

Mathematical models that studied COVID-19 infections are enormous in the literature. We highlight some of them in the present study. For example, the authors utilized the infected cases of coronavirus in UAE and presented a deterministic stochastic model in [4]. A mathematical study that investigated the coronavirus infection using asymptomatic and symptomatic classes is presented in [5]. In [6], the authors considered a mathematical model to emphasize the dynamics of healthcare workers and the community. A numerical study to study the coronavirus infection with environmental contact rates was investigated in [7]. The authors in [8] studied the COVID-19 infection with treatment. The COVID-19 infection with vaccination strategies using reinforcement learning is utilized in [9]. The authors in [10] considered an algorithm to investigate coronavirus infections for undetected individuals. Using the reported cases of coronavirus in Bangkok Thailand through the mathematical study was analyzed in [11]. A mathematical study was given on the COVID-19 disease with fractional derivative modeling for the coronavirus in [12]. A mathematical study to determine the upper bound for new COVID-19 infection in Germany is explored in [13]. Some recent literature related to COVID-19 infection can be seen in [14,15,16,17,18,19]. For example, the discovery of a new variant of the coronavirus called omicron has been modeled by the authors mathematically and the results were established [14]. The stability and the analysis of the second wave of the coronavirus infection were investigated in [15]. A non-integer system to study the COVID-19 infection in India was considered in [16]. The control theory was used in [17] by the authors to show some strategies for the elimination of the virus. The COVID-19 pandemic with the implementation of the control strategies was investigated in [18]. In [19], a fractional order was considered to investigate the COVID-19 infection, and the results regarding disease elimination were presented.

The concept of stochastic differential equations (SDEs) and their applications has been reported in the literature for various physical and biological problems. Some of the recent literature related to COVID-19 infection is presented here, where the concept of stochastic differential equations is utilized. For instance, the authors in [20] used stochastic modeling to understand the coronavirus infection dynamics. The authors in [21] consider the stochastic modeling with time delay and obtained results for the COVID-19 infection. Stochastic differential equations and their application to the study of COVID-19 real data have been obtained in [22]. Analysis of the COVID-19 infection using stochastic modeling with Jump was analyzed in [23]. The coronavirus infection and its detection in the workplace using the concept of SDEs was considered in [24]. The important theory about stochastic modeling, formulation, analysis, and their numerical simulations have been discussed in detail for stochastic differential equations [25]. We followed this work [25] with care and present our model results.

This paper focuses on the transmission dynamics of COVID-19 disease in stochastic differential equations. We study the basic properties related to the deterministic model and then present in detail the analysis for the stochastic case. We utilize the real infection data of COVID-19 in Saudi Arabia and fit well the data into the model. We consider the data fitting using the deterministic and stochastic differential equations model and show that the stochastic case provides better fitting. Some important graphical results for the model are presented for the sensitive parameters that can be useful for disease elimination. The rest of the work in this paper is organized as follows: The formulation of the new COVID-19 model in deterministic and stochastic differential equations is discussed in Section 2. In Section 3, we discuss the analysis of the deterministic model. In Section 4, we carry out the mathematical results for the stochastic model. Section 5 discusses the parameter estimations for the model. Section 6 and Section 7, respectively, present the numerical results and their discussion, and the conclusions of the present work.

2. Model Formulation

The COVID-19 infection that caused many infections and deaths throughout the world still needs proper attention from researchers around the world to create some useful controls, treatments, and effective vaccinations for its elimination. Mathematical modeling in this regard provides interesting information about the disease’s current propagation and its future possible controls. In this section, we consider a mathematical model for understanding the coronavirus infection using the stochastic ordinary nonlinear differential equation model. We consider the total humans denoted by $N\left(t\right)$ into five components: the healthy or the susceptible population $S\left(t\right)$ (the individuals in this compartment are not yet infected or immune, but capable of attracting infection while exposed to symptomatic, asymptomatic or exposed individuals); the healthy individuals that have close contact with infected or become exposed to the disease are shown by $E\left(t\right)$; individuals in the exposed class who complete their incubation and show clear disease symptoms are symptomatically infected and are given by $I\left(t\right)$; while those who do not show clear disease symptoms are known as asymptomatic infected and are described by $A\left(t\right)$; the people who are recovered either from the symptomatic, asymptomatic or exposed stage are recovered and join the recovered class, $R\left(t\right)$, so that $N\left(t\right)=S\left(t\right)+E\left(t\right)+I\left(t\right)+A\left(t\right)+R\left(t\right)$. It should be noted that the exposed individuals have the ability to infect other healthy individuals, and this route of transmission has been included in this paper [26]. It is obvious that asymptomatic people that do not have obvious disease symptoms play a significant role in disease transmission, and it is very difficult to control the disease until we can increase the number of testing for COVID-19 people. The spread due to exposed individuals has been documented in several clinical studies [27,28,29,30]. The above discussion can be shown mathematically in the form of a nonlinear model governed by the evolutionary differential equations given by:

$$\begin{array}{ccc}\frac{dS}{dt}\hfill & =\hfill & \Lambda -({\beta }_{1}E+{\beta }_{2}I+{\beta }_{3}A)\frac{S}{N}-dS,\hfill \\ \frac{dE}{dt}\hfill & =\hfill & ({\beta }_{1}E+{\beta }_{2}I+{\beta }_{3}A)\frac{S}{N}-(\delta +d)E,\hfill \\ \frac{dI}{dt}\hfill & =\hfill & (1-\tau )\delta E-(d+d_1+{\gamma }_1)I,\hfill \\ \frac{dA}{dt}\hfill & =\hfill & \tau \delta E-(d+{\gamma }_2)A,\hfill \\ \frac{dR}{dt}\hfill & =\hfill & {\gamma }_{1}I+{\gamma }_{2}A-dR,\hfill \end{array}$$

(1)

where the non-negative initial values of the model variables are

$$S\left(0\right)=S_0\ge 0,E\left(0\right)=E_0\ge 0,I\left(0\right)=I_0\ge 0,A\left(0\right)=A_0\ge 0,R\left(0\right)=R_0\ge 0.$$

(2)

The population of healthy people is increased through the birth rate $\Lambda$ while it is decreased through the natural mortality rate d. In each compartment of the model, the natural death rate is given by d. The healthy individuals after interacting with exposed individuals can become infective after completing their incubation period through the effective contact rate ${\beta }_1$, and the route of the transmission for this is given by ${\beta }_{1}SE/N$. The healthy individuals that interact with infected people (that show clinical symptoms) and become infected through the contact rate are given by ${\beta }_2$, and the route of the transmission for this rate is given by ${\beta }_{2}SI/N$. Asymptomatic individuals with no clinical symptoms are considered more dangerous than those showing clinical symptoms who come in contact with healthy people and become infected with a rate given by ${\beta }_3$, so the route of the transmission for this is given by ${\beta }_{3}SA/N$. The parameter $\delta$ defines the successful incubation period of the individuals after coming in contact with exposed, symptomatic, or asymptomatic individuals and after completing the incubation period become infected and are either asymptomatic or symptomatic. At the rate $(1-\tau )\delta$, the exposed individuals become symptomatically infected while $\tau \delta$ are the individuals that join the asymptomatic individuals class. The recovery from the infection in the symptomatic or asymptomatic classes is shown, respectively, by ${\gamma }_1$ and ${\gamma }_2$. The natural death due to infection of the individuals in the symptomatic class is given by $d_1$.

Using the concept of stochastic differential equations, one can represent the ordinary differential equations system, given by (1), in the form given by:

$$\begin{array}{ccc}dS\hfill & =\hfill & \left(\Lambda -({\beta }_{1}E+{\beta }_{2}I+{\beta }_{3}A)\frac{S}{N}-dS\right)dt+{\sigma }_{1}Sd{B}_1\left(t\right),\hfill \\ dE\hfill & =\hfill & \left(({\beta }_{1}E+{\beta }_{2}I+{\beta }_{3}A)\frac{S}{N}-(\delta +d)E\right)dt+{\sigma }_{2}Ed{B}_2\left(t\right),\hfill \\ dI\hfill & =\hfill & \left((1-\tau )\delta E-(d+d_1+{\gamma }_1)I\right)dt+{\sigma }_{3}Id{B}_3\left(t\right),\hfill \\ dA\hfill & =\hfill & \left(\tau \delta E-(d+{\gamma }_2)A\right)dt+{\sigma }_{4}Ad{B}_4\left(t\right),\hfill \\ dR\hfill & =\hfill & \left({\gamma }_{1}I+{\gamma }_{2}A-dR\right)dt+{\sigma }_{5}Rd{B}_5\left(t\right),\hfill \end{array}$$

(3)

where ${\sigma }_i$ for $i=1,\cdots 5$ denote the real constant that describes the intensity of the stochastic differential equations while $B_i\left(t\right)$ for $i=1,\cdots ,5$ are refereed to be the stochastic Brownian motion. Due to the biological meaning of the classes $\left(S\right(t),E(t),I(t),A(t),R(t\left)\right)$, our focus is to study the model in the first quadrant:

$${\mathbb{R}}_{+}^5=\left\{(S,E,I,A,R)\in {\mathbb{R}}^5:S\ge 0,E\ge 0,I\ge 0,A\ge 0,R\ge 0\right\}.$$

In the following section, we first investigate the dynamics of the model (3) in the absence of stochastic noises.

3. Dynamics of the Deterministic Model

This section studies the dynamics of the deterministic system when ${\sigma }_1={\sigma }_2={\sigma }_3={\sigma }_4={\sigma }_5=0$, so the model (1) can be easily obtained. We now prove the existence of the solution of the system (1).

3.1. Existence and Uniqueness

We present here the existence and uniqueness of the system (1). We assume that for every $t\in [0,T]$ (T is refereed “final time”), the functions $f_S,f_E,f_I,f_A$ and $f_R$ are bounded, i.e., ${\parallel S\parallel }_{\infty }<M_S$, ${\parallel E\parallel }_{\infty }<M_E$, ${\parallel I\parallel }_{\infty }<M_I$, ${\parallel A\parallel }_{\infty }<M_A$, and ${\parallel R\parallel }_{\infty }<M_R$. We write the model (1) in the form given by:

$$\begin{array}{ccc}{f}_S(t,\Psi )\hfill & =\hfill & \Lambda -({\beta }_{1}E+{\beta }_{2}I+{\beta }_{3}A)\frac{S}{N}-dS,\hfill \\ f_E(t,\Psi )\hfill & =\hfill & ({\beta }_{1}E+{\beta }_{2}I+{\beta }_{3}A)\frac{S}{N}-(\delta +d)E,\hfill \\ f_I(t,\Psi )\hfill & =\hfill & (1-\tau )\delta E-(d+d_1+{\gamma }_1)I,\hfill \\ f_A(t,\Psi )\hfill & =\hfill & \tau \delta E-(d+{\gamma }_2)A,\hfill \\ f_R(t,\Psi )\hfill & =\hfill & {\gamma }_{1}I+{\gamma }_{2}A-dR,\hfill \end{array}$$

where $\Psi =(S,E,I,A,R)$. The parameters of the system (1) positive because it is dealing with the human population. We now show that these functions $f_S$,…$f_R$ hold the linear growth property and the Lipschitz condition. Let us start with linear growth property,

$$\begin{array}{ccc}|f_S(t,\Psi )|\hfill & =\hfill & |\Lambda -({\beta }_{1}E+{\beta }_{2}I+{\beta }_{3}A)\frac{S}{N}-dS|,\hfill \\ \hfill & <\hfill & |\Lambda -({\beta }_{1}E+{\beta }_{2}I+{\beta }_{3}A)S-dS|,\hfill \\ \hfill & \le \hfill & \Lambda +\left({\beta }_1\right|E|+{\beta }_2|I|+{\beta }_3|A\left|\right)\left|S\right|+d\left|S\right|,\hfill \\ \hfill & \le \hfill & \Lambda +({\beta }_1{\displaystyle \underset{t\in D_E}{sup}}\left|E\left(t\right)\right|+{\beta }_2{\displaystyle \underset{t\in D_I}{sup}}\left|I\left(t\right)\right|+{\beta }_3{\displaystyle \underset{t\in D_A}{sup}}\left|A\left(t\right)\right|){\displaystyle \underset{t\in D_S}{sup}}\left|S\left(t\right)\right|\hfill \\ & \hfill & +d{\displaystyle \underset{t\in D_S}{sup}}\left|S\left(t\right)\right|,\hfill \\ \hfill & \le \hfill & \Lambda +({\beta }_1{M}_E+{\beta }_2{M}_I+{\beta }_3{M}_A)M_S+d{M}_S,\hfill \\ \hfill & =\hfill & {M_S}_S,\hfill \\ \hfill & <\hfill & \infty ,\hfill \end{array}$$

$$\begin{array}{ccc}|f_E(t,\Psi )|\hfill & \le \hfill & |({\beta }_{1}E+{\beta }_{2}I+{\beta }_{3}A)S-(\delta +d)E|,\hfill \\ \hfill & \le \hfill & \left({\beta }_1\right|E|+{\beta }_2|I|+{\beta }_3|A\left|\right)\left|S\right|+(\delta +d)\left|E\right|,\hfill \\ \hfill & \le \hfill & ({\beta }_1{\displaystyle \underset{t\in D_E}{sup}}\left|E\left(t\right)\right|+{\beta }_2{\displaystyle \underset{t\in D_I}{sup}}\left|I\left(t\right)\right|+{\beta }_3{\displaystyle \underset{t\in D_A}{sup}}\left|A\left(t\right)\right|){\displaystyle \underset{t\in D_S}{sup}}\left|S\left(t\right)\right|\hfill \\ \hfill & \hfill & +(\delta +d){\displaystyle \underset{t\in D_E}{sup}}\left|E\left(t\right)\right|,\hfill \\ \hfill & \le \hfill & ({\beta }_1{M}_E+{\beta }_2{M}_I+{\beta }_3{M}_A)M_S+(\delta +d)M_E,\hfill \\ \hfill & =\hfill & {M_E}_E,\hfill \\ \hfill & <\hfill & \infty ,\hfill \end{array}$$

$$\begin{array}{ccc}|f_I(t,\Psi )|\hfill & =\hfill & |(1-\tau )\delta E-(d+d_1+{\gamma }_1)I|,\hfill \\ \hfill & \le \hfill & (1-\tau )\delta \left|E\right|+(d+d_1+{\gamma }_1)\left|I\right|,\hfill \\ \hfill & \le \hfill & (1-\tau )\delta {\displaystyle \underset{t\in D_E}{sup}}\left|E\left(t\right)\right|+(d+d_1+{\gamma }_1){\displaystyle \underset{t\in D_I}{sup}}\left|I\left(t\right)\right|,\hfill \\ \hfill & \le \hfill & (1-\tau )\delta M_E+(d+d_1+{\gamma }_1)M_I,\hfill \\ \hfill & =\hfill & {M_I}_I,\hfill \\ \hfill & <\hfill & \infty .\hfill \end{array}$$

$$\begin{array}{ccc}|f_A(t,\Psi )|\hfill & =\hfill & |\tau \delta E-(d+{\gamma }_2\left)A\right|,\hfill \\ \hfill & \le \hfill & \tau \delta \left|E\right|+(d+{\gamma }_2)\left|A\right|,\hfill \\ \hfill & \le \hfill & \tau \delta {\displaystyle \underset{t\in D_E}{sup}}\left|E\left(t\right)\right|+(d+{\gamma }_2){\displaystyle \underset{t\in D_A}{sup}}\left|A\left(t\right)\right|,\hfill \\ \hfill & \le \hfill & \tau \delta M_E+(d+{\gamma }_2)M_A,\hfill \\ \hfill & =\hfill & {M_A}_A,\hfill \\ \hfill & <\hfill & \infty .\hfill \end{array}$$

$$\begin{array}{ccc}|f_R(t,\Psi )|\hfill & =\hfill & |{\gamma }_{1}I+{\gamma }_{2}A-dR|,\hfill \\ \hfill & \le \hfill & {\gamma }_1\left|I\right|+{\gamma }_2\left|A\right|+d\left|R\right|,\hfill \\ \hfill & \le \hfill & {\gamma }_1{\displaystyle \underset{t\in D_I}{sup}}\left|I\left(t\right)\right|+{\gamma }_2{\displaystyle \underset{t\in D_A}{sup}}\left|A\left(t\right)\right|+d{\displaystyle \underset{t\in D_R}{sup}}\left|R\left(t\right)\right|,\hfill \\ & \le \hfill & {\gamma }_1{M}_I+{\gamma }_2{M}_A+d{M}_R,\hfill \\ \hfill & =\hfill & {M_R}_R,\hfill \\ & <\hfill & \infty .\hfill \end{array}$$

Now, we prove the Lipschitz conditions. Let us start with the first equation of the system (1):

$$\begin{array}{ccc}|f_S(t,{\Psi }_{S_1})-f_S(t,{\Psi }_{S_2})|\hfill & =\hfill & |({\beta }_{1}E+{\beta }_{2}I+{\beta }_{3}A)(S_2-S_1)+d(S_2-S_1)|,\hfill \\ & \le \hfill & |({\beta }_{1}E+{\beta }_{2}I+{\beta }_{3}A)(S_2-S_1)|+d|S_2-S_1|,\hfill \\ & \le \hfill & \left|\right({\beta }_1\left|E\right|+{\beta }_2\left|I\right|+{\beta }_3\left|A\right|\left)\right||S_1-S_2|+d|S_1-S_2|,\hfill \\ & \le \hfill & ({\beta }_1{M}_E+{\beta }_2{M}_I+{\beta }_3{M}_A+d)|S_1-S_2|,\hfill \\ & \le \hfill & K_S|S_1-S_2|,\hfill \end{array}$$

where $K_S=({\beta }_1{M}_E+{\beta }_2{M}_I+{\beta }_3{M}_A+d)$.

$$\begin{array}{ccc}|f_E(t,{\Psi }_{E_1})-f_E(t,{\Psi }_{E_2})|\hfill & =\hfill & |{\beta }_{1}S(E_1-E_2)+(\delta +d)(E_2-E_1)|,\hfill \\ & \le \hfill & {\beta }_1\left|S\right||E_1-E_2|+(\delta +d)|E_1-E_2|,\hfill \\ & \le \hfill & ({\beta }_1{\displaystyle \underset{t\in D_S}{sup}}\left|S\left(t\right)\right|+\delta +d\left)\right|E_1-E_2|,\hfill \\ & \le \hfill & ({\beta }_1{M}_S+\delta +d)|E_1-E_2|,\hfill \\ \hfill & \le \hfill & K_E|E_1-E_2|,\hfill \end{array}$$

where $K_E=({\beta }_1{M}_S+\delta +d)$.

$$\begin{array}{ccc}|f_I(t,{\Psi }_{I_1})-f_I(t,{\Psi }_{I_2})|\hfill & \le \hfill & |(d+d_1+{\gamma }_1)|I_1-I_2|,\hfill \\ & \le \hfill & K_I|I_1-I_2|,\hfill \end{array}$$

where $K_I=({\beta }_1{M}_S+\delta +d)$.

$$\begin{array}{ccc}|f_A(t,{\Psi }_{A_1})-f_I(t,{\Psi }_{A_2})|\hfill & \le \hfill & (d+{\gamma }_2)|A_1-A_2|,\hfill \\ & \le \hfill & K_A|A_1-A_2|,\hfill \end{array}$$

where $K_A=(d+{\gamma }_2)$.

$$\begin{array}{ccc}|f_R(t,{\Psi }_{R_1})-f_R(t,{\Psi }_{R_2})|\hfill & \le \hfill & d|R_1-R_2|,\hfill \\ \hfill & \le \hfill & K_R|R_1-R_2|,\hfill \end{array}$$

where $K_R=d$. Therefore, the given system (1) satisfies the linear growth as well as the Lipschitz conditions. Hence, the model (1) admits a unique system of solutions.

3.2. Equilibrium Points Analysis

The dynamics of the model (1) can be studied while investigating their possible equilibrium points. Usually, the disease epidemic models associated with human populations, often involve the disease free and the endemic equilibrium. We find the disease-free equilibrium of the model (1) denoted by $E_0=(S^0,0,0,0,0)$, and it is given by

$$E_0=(\frac{\Lambda }{d},0,0,0,0).$$

The local asymptotic stability of an epidemic model at its equilibrium points is determined with the available threshold number. For disease elimination, the value of the threshold should be less than 1, while for the disease to exist in the population permanently, it should be greater than unity. To study the asymptotic stability analysis of the model, we first investigate the basic reproduction number. The basic reproduction number can be stated as “the number of people infected by a single infected person when introduced into a healthy population by generating secondary infection”. Usually, the basic reproduction number is denoted by ${\mathcal{R}}_0$, and for our proposed model (1), we use the next-generation approach given in [31] to obtain the basic reproduction number. According to this approach [31], we obtain the following matrices,

$$F=\left(\begin{array}{ccc}{\beta }_1& {\beta }_2& {\beta }_3\\ 0& 0& 0\\ 0& 0& 0\end{array}\right),V=\left(\begin{array}{ccc}d+\delta & 0& 0\\ \delta (\tau -1)& d+d_1+{\gamma }_1& 0\\ -\delta \tau & 0& d+{\gamma }_2\end{array}\right).$$

The spectral radius that gives the basic reproduction number ${\mathcal{R}}_0$ is given by $\rho \left(F{V}^{-1}\right)$:

$${\mathcal{R}}_0=\underset{{\mathcal{R}}_1}{\underbrace{\frac{{\beta }_1}{d+\delta }}}+\underset{{\mathcal{R}}_2}{\underbrace{\frac{{\beta }_2\delta (1-\tau )}{\left({\gamma }_1+d+d_1\right)(d+\delta )}}}+\underset{{\mathcal{R}}_3}{\underbrace{\frac{{\beta }_3\delta \tau }{\left({\gamma }_2+d\right)(d+\delta )}}}.$$

Here, ${\mathcal{R}}_1$ accounts for the infection generated through the exposed, ${\mathcal{R}}_2$ through the symptomatic, and ${\mathcal{R}}_3$ through the asymptomatic individuals.

3.3. Endemic Equilibria

We find the expression for the endemic equilibria and their existence. In the endemic case, the disease persists in the population for ${\mathcal{R}}_0>1$. For the existence of unique positive endemic equilibria of the model proposed, the requirement is that the basic reproduction number ${\mathcal{R}}_0>1$. Here, we find the model (1) by denoting it $E_1$ and obtain the following,

$$E_1=\left(S^{*},E^{*},I^{*},A^{*},R^{*}\right)$$

where

$$\left\{\begin{array}{c}{S}^{*}=\frac{\Lambda }{d+{\lambda }^{*}}\hfill \\ E^{*}=\frac{\lambda S^{*}}{d+\delta }\hfill \\ I^{*}=\frac{\delta (1-\tau )E^{*}}{{\gamma }_1+d+d_1}\hfill \\ A^{*}=\frac{\delta \tau E^{*}}{{\gamma }_2+d}\hfill \\ R^{*}=\frac{{\gamma }_2{A}^{*}+{\gamma }_1{I}^{*}}{d}.\hfill \end{array}\right.$$

where

$${\lambda }^{*}=\frac{{\beta }_{1}E+{\beta }_{2}I+{\beta }_{3}A}{N}.$$

Using the above into

$${\lambda }^{*}=\frac{{\beta }_{1}E+{\beta }_{2}I+{\beta }_{3}A}{N},$$

we obtain the following,

$$b_0{\lambda }^{*}+b_1=0,$$

where

$$\begin{array}{ccc}{b}_0\hfill & =\hfill & \left({\gamma }_2+d\right)\left(\left({\gamma }_1+d\right)(d+\delta )+d_1(\delta \tau +d)\right),\hfill \\ b_1\hfill & =\hfill & d\left({\gamma }_1+d+d_1\right)\left({\gamma }_2+d\right)(d+\delta )(1-{\mathcal{R}}_0).\hfill \end{array}$$

The endemic positive equilibrium exists only if ${\mathcal{R}}_0>1$. The disease will spread in the community if ${\mathcal{R}}_0>1$ otherwise it will vanish if ${\mathcal{R}}_0\le 1$.

Lemma 1.

The plane given by $S+E+I+A+R=\frac{\Lambda }{d}$ is a manifold of model (1), and it is attracting in the first octant.

Proof. 

Adding all five equations in system (1), we obtain

$$\frac{dN}{dt}=\Lambda -dN-d_{1}I\le \Lambda -dN.$$

(4)

It is obvious that $N\left(t\right)=\Lambda /d$ is the solution of (4) and for any $N\left(t\right)\ge 0$, we may obtain the general solution for (4) given by

$$N\left(t\right)=\frac{\Lambda }{d}+\left(N\left(0\right)-\frac{\Lambda }{d}\right)e^{-dt}.$$

(5)

So,

$$\underset{t\to \infty }{lim}N\left(t\right)=\frac{\Lambda }{d},$$

(6)

that summarizes the result. □

3.4. Global Stability Disease Free Case

The following theorem shows the stability of the model (1) in the disease-free case. The global stability of epidemic models is important in the sense that there are no possibilities of backward bifurcation if the model is globally asymptotically stable.

Theorem 1.

The deterministic model (1) is globally asymptotically stable when ${\mathcal{R}}_0\le 1$.

Proof. 

We define the Lyapunov function given by

$$L=({\gamma }_1+d+d_1)E+{\beta }_{2}I+\frac{({\gamma }_1+d+d_1)\left({\beta }_3\right)}{({\gamma }_2+d)}A.$$

Time differentiation along the solution of the model (1) gives the following:

$$\begin{array}{ccc}\dot{L}\hfill & =\hfill & ({\gamma }_1+d+d_1)\dot{E}+{\beta }_2\dot{I}+\frac{({\gamma }_1+d+d_1)\left({\beta }_3\right)}{({\gamma }_2+d)}\dot{A},\hfill \\ \hfill & =\hfill & ({\gamma }_1+d+d_1)[({\beta }_{1}E+{\beta }_{2}I+{\beta }_{3}A)\frac{S}{N}-(\delta +d)E]\hfill \\ & \hfill & +{\beta }_2[(1-\tau )\delta E-(d+d_1+{\gamma }_1)I]\hfill \\ & \hfill & +\frac{({\gamma }_1+d+d_1)\left({\beta }_3\right)}{({\gamma }_2+d)}[\tau \delta E-(d+{\gamma }_2)A],\hfill \\ \hfill & \le \hfill & \left(\frac{{\beta }_1}{\delta +d}+\frac{{\beta }_2\delta (1-\tau )}{(\delta +d)({\gamma }_1+d+d_1)}+\frac{{\beta }_3\tau \delta }{({\gamma }_2+d)(\delta +d)}-1\right)({\gamma }_1+d+d_1)(\delta +d)E,\hfill \\ \hfill & =\hfill & ({\mathcal{R}}_0-1)({\gamma }_1+d+d_1)(\delta +d)E.\hfill \end{array}$$

So,

$$\begin{array}{ccc}\dot{L}\hfill & \le \hfill & ({\mathcal{R}}_0-1)({\gamma }_1+d+d_1)(\delta +d)E.\hfill \end{array}$$

Here, $\dot{L}\le 0$ if and only if ${\mathcal{R}}_0\le 1$. Thus, the deterministic model given by (1) is globally asymptotically stable. □

4. Dynamics of the Stochastic Model

To study the dynamic behavior of stochastic models (3), we first discuss whether the solution of the proposed model is non-negative and whether it concerns global existence. We need to follow [32]. Now let us take a closer look at the dynamics of the probabilistic model (3). First, we present some important results that will help us in later results.

4.1. Preliminaries

This subsection provides a preliminary concept regarding stochastic calculus and some other useful theorem that is based on the onward results in the present paper. We assumed that $(\Xi ,\mathcal{Q},{\left\{{\mathcal{Q}}_t\right\}}_{t\ge 0},\mathbb{P})$ defined to be a complete probability space with the filtration ${\left\{{\mathcal{Q}}_t\right\}}_{t\ge 0}$. Then, we introduce the following notations:

$$\begin{array}{cc}\hfill {\mathbb{R}}_{+}^m& =\{u=(u_1,\cdots ,u_m)\in {\mathbb{R}}^m;u_k>0,1\le k\le m\},\hfill \\ \hfill {\overline{\mathbb{R}}}_{+}^m& =\{u=(u_1,\cdots ,u_m)\in {\mathbb{R}}^m;u_k\ge 0,1\le k\le m\}.\hfill \end{array}$$

The d-dimensional stochastic ordinary differential equation [32] can be written as:

$$du\left(t\right)=\Phi \left(u\left(t\right)\right)dt+v\left(u\left(t\right)\right)dG\left(t\right)fort\ge t_0,$$

(7)

with the initial condition $u\left(0\right)=u_0\in {\mathbb{R}}^m$ where, $G\left(t\right)$ is the d-dimensional standard Brownian motion that is defined on the given probability space. Suppose $C^2({\mathbb{R}}^m;{\overline{\mathbb{R}}}_{+})$ describes the family of all functions that are non-negative $\mathcal{V}$ and defined on ${\mathbb{R}}^m$, which has the property to be continuously differentiable twice in u. Follows [32], we can present for model (7), the differential operator $\mathcal{L}$ shown by:

$$\mathcal{L}=\sum _{k=1}^m{\Phi }_k(u,t)\frac{\partial }{\partial u_k}+\frac{1}{2}\sum _{k,l=1}^m{\left[v^T(u,t)v(u,t)\right]}_{k,l}\frac{{\partial }^2}{\partial u_k\partial u_l}.$$

When $\mathcal{L}$ acts on a function $\mathcal{V}\in C^2({\mathbb{R}}^m;{\overline{\mathbb{R}}}_{+})$, then

$$\mathcal{L}\mathcal{V}=V_u\left(u\right)\Phi \left(u\right)+\frac{1}{2}trace\left[v^T\left(u\right)V_{uu}\left(u\right)g\left(u\right)\right],$$

where, $V_u=\left(\frac{\partial V}{\partial u_1},\cdots ,\frac{\partial V}{\partial u_m}\right),V_{uu}={\left(\frac{{\partial }^{2}V}{\partial u_k\partial u_l}\right)}_{m\times m}$.

By Ito’s lemma [32], if $u\in {\mathbb{R}}^m$, we have

$$d\mathcal{V}\left(u\left(t\right)\right)=\mathcal{L}\mathcal{V}\left(u\left(t\right)\right)dt+V_u\left(u\left(t\right)\right)v\left(u\left(t\right)\right)dG\left(t\right).$$

Consider that $B\left(t\right)$ defines the Brownian motion while $Z\left(t\right)$ be an Ito drift-diffusion process that satisfies the following stochastic differential equation:

$$dZ\left(t\right)=\mu \left(Z\right(t),t)dt+\sigma \left(Z\right(t),t)dG\left(t\right).$$

If $f(z,t)\in C^2\left({\mathbb{R}}^2,\mathbb{R}\right)$ then $\Phi \left(Z\right(t),t)$ is also an Ito drift-diffusion process, with their differential, given by:

$$d(\Phi (Z\left(t\right),t))=\frac{\partial \Phi }{\partial t}(Z\left(t\right),t)dt+{\Phi }^{\prime }(Z\left(t\right),t)dZ+\frac{1}{2}{\Phi }^{″}(Z\left(t\right),t)dZ{\left(t\right)}^2,$$

with $dZ{\left(t\right)}^2$ given by: $d{t}^2=0,dtdG\left(t\right)=0$ and $dG{\left(t\right)}^2=dt$.

Now, we give the following result:

4.2. Existence of the Positive Unique Global Solution

We provide the results for the unique positive global solution existence of the system (3) which are shown in the following:

Theorem 2.

For the initial values of the model variables, $X\left(0\right)=(S\left(0\right),E\left(0\right),I\left(0\right),A\left(0\right),R\left(0\right))\in R_{+}^5$, a solution (non-negative) exists, say, $X\left(t\right)=\left(S\right(t),E(t),I(t),A(t),R(t\left)\right)$ for the system (3) on $t\ge 0$ and it is remains in $R_{+}^5$ with the unit probability.

Proof. 

The coefficients given in (3) of all the equations are locally continuous in Lipschitz sense for the initial data $(S\left(0\right),E\left(0\right),I\left(0\right),A\left(0\right),R\left(0\right))\in {\mathbb{R}}_{+}^5$, so there should be a unique solution or, in other words, a local solution, say $\left(S\right(0),E(0),I(0),A(0),R(0\left)\right)$ on $t\in [0,{\chi }_m)$, where ${\chi }_m$ defines to be the explosion time, the readers can see more details in the work [22,23,24,33,34]. To obtain a global solution, one has to show ${\chi }_m=\infty$. We need to choose a positive real number, say, $j_0$, to be large enough so that all the initial values associated with the states should lie in $[\frac{1}{j_0},j_0]$. Further, the stopping time can be defined as,

$$\begin{array}{ccc}{\chi }_j\hfill & =\hfill & \{t\in [0,{\chi }_m):\frac{1}{j}\ge min\{S\left(t\right),E\left(t\right),I\left(t\right),A\left(t\right),R\left(t\right)\}\hfill \\ & \hfill & \mathrm{or}max\left\{S\right(t),E(t),I(t),A(t),R(t\left)\right\}\ge j\}\hfill \end{array}$$

(8)

for every non-negative integer j ($j\ge j_0$). Considering $inf\varphi =\infty$ for the case $\varphi$ when it shows the null set. The definition of stopping time ensures that ${\chi }_j$ has the monotonically increasing behavior for $j\to \infty$. Consider, ${lim}_{j\to \infty }={\chi }_{\infty }$ with ${\chi }_m\ge {\chi }_{\infty }$ a.s.

For the non-negative value of t, we consider the case and claim that ${\chi }_{\infty }=\infty$ a.s., then it can be said that ${\chi }_m=\infty$ and a.s. $(S\left(t\right),E\left(t\right),I\left(t\right),A\left(t\right),R\left(t\right))\in R_{+}^5$. So, we have to show that ${\chi }_m=\infty$ a.s. if such a result is wrong then there may exists two constants such as $0<T$ and $ϵ\in (0,1)$ such that

$$P\{T\ge {\chi }_{\infty }\}>ϵ.$$

(9)

Define the following Lyapunov function

$$H(S,E,I,A,R)=S+I+E+A+R-5-(\log S+\log E+\log I+\log A+\log R).$$

(10)

Applying $It\widehat{o}$ formula on the Equation (10), we obtain,

$$\begin{array}{cc}\hfill dH(S,E,I,A,R)& =\left(1-\frac{1}{S}\right)dS+{\sigma }_1(S-1)d{B}_1\left(t\right)+\left(1-\frac{1}{E}\right)dE+{\sigma }_2(E-1)d{B}_2\left(t\right)\hfill \\ \hfill & +\left(1-\frac{1}{I}\right)I+{\sigma }_3(I-1)d{B}_3\left(t\right)+\left(1-\frac{1}{A}\right)A+{\sigma }_3(A-1)d{B}_4\left(t\right)\hfill \\ \hfill & +\left(1-\frac{1}{R}\right)dR+{\sigma }_5(R-1)d{B}_5\left(t\right),\hfill \\ \hfill & =LH(S,E,I,A,R)dt+{\sigma }_1(S-1)d{B}_1\left(t\right)+{\sigma }_2(E-1)d{B}_2\left(t\right)\hfill \\ \hfill & +{\sigma }_3(I-1)d{B}_3\left(t\right)+{\sigma }_4(A-1)d{B}_4\left(t\right)+{\sigma }_5(R-1)d{B}_5\left(t\right).\hfill \end{array}$$

(11)

In Equation (11), the relation $LH:R_{+}^5\to R_{+}$ can be shown to the form given by,

$$\begin{array}{cc}\hfill LH(S,E,I,A,R)& =\left(1-\frac{1}{S}\right)\left(\Lambda -\left({\beta }_{1}E+{\beta }_{2}I+{\beta }_{3}A\right)\frac{S}{N}-dS\right)\hfill \\ \hfill & +\left(1-\frac{1}{E}\right)\left(\left({\beta }_{1}E+{\beta }_{2}I+{\beta }_{3}A\right)\frac{S}{N}-(\delta +d)E\right)\hfill \\ \hfill & +\left(1-\frac{1}{I}\right)\left((1-\tau )\delta E-\left(d+d_1+{\gamma }_1\right)I\right)\hfill \\ \hfill & +\left(1-\frac{1}{A}\right)\left(\tau \delta E-\left(d+{\gamma }_2\right)A\right)+\left(1-\frac{1}{R}\right)\left({\gamma }_{1}I+{\gamma }_{2}A-dR\right)\hfill \\ \hfill & +\frac{{\sigma }_1^2+{\sigma }_2^2+{\sigma }_3^2+{\sigma }_4^2+{\sigma }_5^2}{2},\hfill \\ \hfill & =\Lambda -\left({\beta }_{1}E+{\beta }_{2}I+{\beta }_{3}A\right)\frac{S}{N}-dS-\frac{\Lambda }{S}+\left({\beta }_{1}E+{\beta }_{2}I+{\beta }_{3}A\right)\frac{1}{N}+d\hfill \\ \hfill & +(\left({\beta }_{1}E+{\beta }_{2}I+{\beta }_{3}A\right)\frac{S}{N}-(\delta +d)E-(\left({\beta }_{1}E+{\beta }_{2}I+{\beta }_{3}A\right)\frac{S}{NE}+(\delta +d)\hfill \\ \hfill & +(1-\tau )\delta E-(d+d_1+{\gamma }_1)I-\frac{(1-\tau )\delta E}{I}+(d+d_1+{\gamma }_1)\hfill \\ \hfill & +\tau \delta E-\left(d+{\gamma }_2\right)A-\frac{\tau \delta E}{A}+\left(d+{\gamma }_2\right)+{\gamma }_{1}I+{\gamma }_{2}A-dR-\frac{{\gamma }_{1}I}{R}+\frac{{\gamma }_{2}A}{R}+d\hfill \\ \hfill & +\frac{{\sigma }_1^2+{\sigma }_2^2+{\sigma }_3^2+{\sigma }_4^2+{\sigma }_5^2}{2},\hfill \\ \hfill & \le \Lambda +{\beta }_1+{\beta }_2+{\beta }_3+d+(\delta +d)+(d+d_1+{\gamma }_1)+(d+{\gamma }_2)+d\hfill \\ \hfill & +\frac{{\sigma }_1^2+{\sigma }_2^2+{\sigma }_3^2+{\sigma }_4^2+{\sigma }_5^2}{2}\hfill \\ \hfill & \le \Lambda +{\beta }_1+{\beta }_2+{\beta }_3+5d+\delta +d_1+{\gamma }_1+{\gamma }_2\hfill \\ \hfill & +\frac{{\sigma }_1^2+{\sigma }_2^2+{\sigma }_3^2+{\sigma }_4^2+{\sigma }_5^2}{2}:=K.\hfill \end{array}$$

(12)

The constant K is obviously positive, which is free from independent state variables. Therefore,

$$\begin{array}{cc}\hfill dH(S,E,I,A,R)& \le Kdt+{\sigma }_1(S-1)d{B}_1\left(t\right)+{\sigma }_{2}Ed{B}_2\left(t\right)\hfill \\ \hfill & +{\sigma }_3(I-1)d{B}_3\left(t\right)+{\sigma }_4(A-1)d{B}_3\left(t\right)+{\eta }_5(R-1)d{B}_5\left(t\right).\hfill \end{array}$$

(13)

Integrating both sides of Equation (13) from 0 to ${\chi }_j\wedge T$

$$\begin{array}{cc}\hfill & E\left[H(S({\chi }_j\wedge T),E({\chi }_j\wedge T),I({\chi }_j\wedge T),A({\chi }_j\wedge T),R({\chi }_j\wedge T))\right],\hfill \\ \hfill & \le H(S\left(0\right),E\left(0\right),I\left(0\right),A\left(0\right),R\left(0\right))+E\left[{\int }_0^{{\chi }_j\wedge T}K\right],\hfill \\ \hfill & \le H\left(S\right(0),E(0),I(0),A(0),R(0\left)\right)+TK.\hfill \end{array}$$

(14)

Setting ${\Omega }_j=\{T\ge {\chi }_j\}$ for $j_1\le j$ and thus by Equation (9), $P\left({\Omega }_j\right)\ge ϵ$. It should be kept in mind that for every $\omega$ in ${\Omega }_j$, there should be at least one $S({\chi }_j,\omega )$, $E({\chi }_j,\omega )$, $I({\chi }_j,\omega )$,$A({\chi }_j,\omega )$,$R({\chi }_j,\omega )$, that equals $\frac{1}{j}$ or j. Hence, $(S\left({\chi }_j\right),E\left({\chi }_j\right),I\left({\chi }_j\right),I\left({\chi }_j\right),R\left({\chi }_j\right))$ is not less than $j-logj-1$ or $logj-1+\frac{1}{j}$. As a result

$$\left(logj-1+\frac{1}{j}\right)\wedge E\left(j-logj-1\right)\le H(S\left({\chi }_j\right),E\left({\chi }_j\right),I\left({\chi }_j\right),A\left({\chi }_j\right),R\left({\chi }_j\right)).$$

(15)

Referring to Equations (9) and (14), we can write

$$\begin{array}{cc}\hfill H\left(S\right(0),E(0),I(0),A(0),R(0\left)\right)+TK& \ge E[1_{\Omega \left(\omega \right)}H\left(S\left({\chi }_j\right),E\left({\chi }_j\right),I\left({\chi }_j\right),A\left({\chi }_j\right),R\left({\chi }_j\right))\right],\hfill \\ \hfill & \ge ϵ\left[(-1+j-logj)\wedge \left(-1+\frac{1}{j}+logj\right)\right].\hfill \end{array}$$

(16)

Here, the notion $1_{\Omega \left(\omega \right)}$ describes the indicator function of $\Omega$. Considering $j\to \infty$ leads to dichotomy $\infty >H\left(S\left(0\right),E\left(0\right),I\left(0\right),A\left(0\right),R\left(0\right))\right)+TK=\infty$, that implies ${\chi }_{\infty }=\infty$ a.s. □

4.3. Extinction for the Proposed Model

Here, we find the conditions for the extinction of the system (3). Some definitions and notations are given here, which shall be utilized later in the proof of the result. Consider,

$$〈x\left(t\right)〉=\frac{1}{t}{\int }_0^{t}x\left(r\right)dr.$$

For the stochastic model (3), we give the following threshold ${\mathbf{R}}_{\mathbf{0}}$:

$${\mathbf{R}}_{\mathbf{0}}=\frac{\beta }{(3d+d_1+{\gamma }_1+{\gamma }_2+\frac{{\sigma }_2^2}{2}+\frac{{\sigma }_3^2}{2}+\frac{{\sigma }_4^2}{2})}.$$

(17)

Theorem 3.

Let us assume that $(S,E,I,A,R)$ shows the solution of the system (3) with the given initial conditions $(S\left(0\right),E\left(0\right),I\left(0\right),A\left(0\right),R\left(0\right))\in R_{+}^5$, then

$$\begin{array}{c}\hfill \underset{t\to \infty }{lim}\frac{A\left(t\right)+E\left(t\right)+I\left(t\right)+A\left(t\right)+R\left(t\right)}{t}=0,\mathrm{a}.\mathrm{s}.\end{array}$$

(18)

then

$$\begin{array}{cc}\hfill & \underset{t\to \infty }{lim}\frac{1}{t}{\int }_0^{t}S\left(u\right)d{B}_1\left(u\right)=0,\hfill \\ \hfill & \underset{t\to \infty }{lim}\frac{1}{t}{\int }_0^{t}Ed{B}_2\left(u\right)=0,\hfill \\ \hfill & \underset{t\to \infty }{lim}\frac{1}{t}{\int }_0^{t}I\left(u\right)d{B}_3\left(u\right)=0,\hfill \\ \hfill & \underset{t\to \infty }{lim}\frac{1}{t}{\int }_0^{t}A\left(u\right)d{B}_4\left(u\right)=0,\hfill \\ \hfill & \underset{t\to \infty }{lim}\frac{1}{t}{\int }_0^{t}R\left(u\right)d{B}_5\left(u\right)=0.\hfill \end{array}$$

(19)

To obtain the proof of Theorem 3, the readers can see the work of Zhao and Jiang [35], which follows the same procedure to obtain the proof of Lemmas $\left(2.1\right)$ and $\left(2.2\right)$. Moreover, one can see in the work [36] the same steps to obtain the result for the Lemma 3. So, the proof of this is omitted.

We give the following result for disease extinction.

Theorem 4.

Let $(S,E,I,A,R)$ denote the solution to the system (3) with the initial conditions given by

$(S\left(0\right),E\left(0\right),I\left(0\right),A\left(0\right),R\left(0\right))\in R^5$. If ${\mathbf{R}}_{\mathbf{0}}<1$, then

$$\underset{t\to \infty }{lim}\frac{〈logE\left(t\right)〉}{t}<0,\underset{t\to \infty }{lim}\frac{〈logI\left(t\right)〉}{t}<0and\underset{t\to \infty }{lim}\frac{〈logA\left(t\right)〉}{t}<0,\mathrm{a}.\mathrm{s}.$$

Which means that disease will die out with probability one. In addition

$$\begin{array}{cc}\hfill \underset{t\to \infty }{lim}〈S\left(t\right)〉& =\frac{\Lambda }{d},\hfill \\ \hfill \underset{t\to \infty }{lim}〈E\left(t\right)〉& =0,\hfill \\ \hfill \underset{t\to \infty }{lim}〈I\left(t\right)〉& =0,\hfill \\ \hfill \underset{t\to \infty }{lim}〈A\left(t\right)〉& =0,\hfill \\ \hfill \underset{t\to \infty }{lim}〈R\left(t\right)〉& =0,\hfill \end{array}$$

(20)

Proof. 

We directly integrate the model (3), and obtain the following equations:

$$\begin{array}{cc}\hfill & \frac{S\left(t\right)-S\left(0\right)}{t}=\Lambda -\left({\beta }_1〈E〉+{\beta }_2〈I〉+{\beta }_3〈A〉\right)\frac{〈S〉}{〈N〉}-d〈S〉+\frac{{\sigma }_1}{t}{\int }_0^{t}S\left(r\right)d{B}_1\left(r\right),\hfill \\ \hfill & \frac{E\left(t\right)-E\left(0\right)}{t}=\left({\beta }_1〈E〉+{\beta }_2〈I〉+{\beta }_3〈A〉\right)\frac{〈S〉}{〈N〉}-(\delta +d)〈E〉+\frac{{\sigma }_2}{t}{\int }_0^{t}E\left(r\right)d{B}_2\left(r\right),\hfill \\ \hfill & \frac{I\left(t\right)-I\left(0\right)}{t}=(1-\tau )\delta 〈E〉-\left(d+d_1+{\gamma }_1\right)〈I〉+\frac{{\sigma }_3}{t}{\int }_0^{t}I\left(r\right)d{B}_3\left(r\right),\hfill \\ \hfill & \frac{A\left(t\right)-I_A\left(0\right)}{t}=\tau \delta 〈E〉-\left(d+{\gamma }_2\right)〈A〉+\frac{{\sigma }_4}{t}{\int }_0^{t}A\left(r\right)d{B}_4\left(r\right),\hfill \\ \hfill & \frac{R\left(t\right)-R\left(0\right)}{t}={\gamma }_1〈I〉+{\gamma }_2〈A〉-d〈R〉+\frac{{\sigma }_5}{t}{\int }_0^{t}R\left(r\right)d{B}_5\left(r\right).\hfill \end{array}$$

(21)

Define $\Phi =E+I+A$, and using the $It\widehat{o}$ formula to the second equation of model (3), we have

$$\begin{array}{cc}\hfill d\log \Phi \left(t\right)& =\frac{1}{\Phi }\left[({\beta }_{1}E+{\beta }_{2}I+{\beta }_{3}A)\frac{S}{N}-(d+d_1+{\gamma }_1)I-(d+{\gamma }_2)A\right]dt\hfill \\ \hfill & -\frac{{\sigma }_2^2{E}^{2}dt}{2{\Phi }^2}-\frac{{\sigma }_3^2{I}^{2}dt}{2{\Phi }^2}-\frac{{\sigma }_4^2{A}^{2}dt}{2{\Phi }^2}+\frac{{\sigma }_{2}Ed{B}_2\left(t\right)}{\Phi }+\frac{{\sigma }_{3}Id{B}_3\left(t\right)}{\Phi }\hfill \\ \hfill & +\frac{{\sigma }_{4}Ad{B}_4\left(t\right)}{\Phi }.\hfill \end{array}$$

(22)

We integrate the Equation (22) in the interval given by $[0,t]$ and then divide with t, which gives

$$\begin{array}{ccc}\frac{\log \Phi \left(t\right)-\log \Phi \left(0\right)}{t}\hfill & =\hfill & \frac{1}{\Phi }[({\beta }_{1}E+{\beta }_{2}I+{\beta }_{3}A)\frac{S}{N}-dE-(d+d_1+{\gamma }_1)I\hfill \\ & \hfill & -(d+{\gamma }_2)A-\frac{{\sigma }_2^2{E}^2}{2{\Phi }^2}-\frac{{\sigma }_3^2{I}^2}{2{\Phi }^2}-\frac{{\sigma }_4^2{A}^2}{2{\Phi }^2}]\hfill \\ & \hfill & +\frac{{\sigma }_2}{t}{\int }_0^t\frac{Ed{B}_2\left(t\right)}{\Phi }+\frac{{\sigma }_3}{t}{\int }_0^t\frac{Id{B}_3\left(t\right)}{\Phi }+\frac{{\sigma }_4}{t}{\int }_0^t\frac{Ad{B}_4\left(t\right)}{\Phi }.\hfill \end{array}$$

(23)

Moreover, $M\left(t\right)=\frac{{\sigma }_2}{t}{\int }_0^t\frac{Ed{B}_2\left(t\right)}{\Phi }+\frac{{\sigma }_3}{t}{\int }_0^t\frac{Id{B}_3\left(t\right)}{\Phi }+\frac{{\sigma }_4}{t}{\int }_0^t\frac{Ad{B}_4\left(t\right)}{\Phi }$, which is a local continuous martingale and $M\left(0\right)=0$. Applying the result given in (3), we have

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

(24)

Then

$$\begin{array}{cc}\hfill \frac{\log \Phi \left(t\right)}{t}& \le (\beta -(d+(d+d_1+{\gamma }_1)+(d+{\gamma }_2)+\frac{{\sigma }_2^2}{2}+\frac{{\sigma }_3^2}{2}+\frac{{\sigma }_4^2}{2})+\frac{\log \Phi \left(0\right)}{t},\hfill \\ \hfill & \le (\beta -(3d+d_1+{\gamma }_1+{\gamma }_2+\frac{{\sigma }_2^2}{2}+\frac{{\sigma }_3^2}{2}+\frac{{\sigma }_4^2}{2})+\frac{\log \Phi \left(0\right)}{t}.\hfill \end{array}$$

(25)

where $\beta =Max\{{\beta }_1,{\beta }_2,{\beta }_3\}$, and we know that $S\le N$, and $E+I+A\le ({\beta }_{1}E+{\beta }_{2}I+{\beta }_{3}A)$, and $\frac{\log \Phi \left(0\right)}{t}=0$.

If ${\mathbf{R}}_{\mathbf{0}}<1$ is satisfied, then Equation (25) becomes

$$\underset{t\to \infty }{lim}\sup \frac{\log \Phi \left(t\right)}{t}\le \left(3d+d_1+{\gamma }_1+{\gamma }_2+\frac{{\sigma }_2^2}{2}+\frac{{\sigma }_3^2}{2}+\frac{{\sigma }_4^2}{2}\right)\left({\mathbf{R}}_{\mathbf{0}}-1\right)<0,\mathrm{a}.\mathrm{s}.$$

(26)

The above Equation (26) implies that

$$\underset{t\to \infty }{lim}〈\Phi \left(t\right)〉=0,\mathrm{a}.\mathrm{s}.$$

(27)

Then

$$\underset{t\to \infty }{lim}〈\Phi \left(t\right)〉=\underset{t\to \infty }{lim}(E+I+A)=0,\mathrm{a}.\mathrm{s}.$$

(28)

Therefore, as a result,

$$\begin{array}{cc}\hfill \underset{t\to \infty }{lim}E& =0,\hfill \\ \hfill \underset{t\to \infty }{lim}I& =0,\hfill \\ \hfill \underset{t\to \infty }{lim}A& =0.\hfill \end{array}$$

(29)

Furthermore, considering the first equation of the model (21) and integrating with the limit ranging from 0 to t, after that dividing by t, and using Equation (29), then, we have

$$\begin{array}{cc}\hfill & \frac{S\left(t\right)-S\left(0\right)}{t}=\Lambda -\left({\beta }_1〈E〉+{\beta }_2〈I〉+{\beta }_3〈A〉\right)\frac{〈S〉}{〈N〉}-d〈S〉+\frac{{\sigma }_1}{t}{\int }_0^{t}S\left(r\right)d{B}_1\left(r\right),\hfill \\ \hfill & 〈S〉=\frac{1}{d}\left[\Lambda -\left({\beta }_1〈E〉+{\beta }_2〈I〉+{\beta }_3〈A〉\right)\frac{〈S〉}{〈N〉}+\frac{S\left(0\right)-S\left(t\right)}{t}+\frac{{\sigma }_1}{t}{\int }_0^{t}S\left(r\right)d{B}_1\left(r\right)\right].\hfill \end{array}$$

(30)

Thus implying that

$$\underset{t\to \infty }{lim}〈S\left(t\right)〉=\frac{\Lambda }{d},\mathrm{a}.\mathrm{s}.$$

(31)

Similarly, the fourth equation of system (21) by integrating within the range from 0 to t and later dividing it by t, and using Equation (29), then we have

$$\begin{array}{cc}\hfill & \frac{R\left(t\right)-R\left(0\right)}{t}={\gamma }_1〈I〉+{\gamma }_2〈A〉-d〈R〉+\frac{{\sigma }_5}{t}{\int }_0^{t}R\left(r\right)d{B}_5\left(r\right),\hfill \\ \hfill & 〈R\left(t\right)〉=\frac{1}{(d+\sigma +{\gamma }_3)}\left[{\gamma }_1〈I〉+{\gamma }_2〈A〉+\frac{R\left(0\right)-R\left(t\right)}{t}+\frac{{\sigma }_5}{t}{\int }_0^{t}R\left(r\right)d{B}_5\left(r\right)\right].\hfill \end{array}$$

(32)

Thus implying that

$$\underset{t\to \infty }{lim}〈R\left(t\right)〉=0,\mathrm{a}.\mathrm{s}.$$

(33)

The prove the result. □

5. Parameters Estimation

This section studies the estimation of the parameters for the model using the nonlinear least square curve fitting method. The data is considered for the infected cases of COVID-19 in Saudi Arabia. The cases are available in worldometers [37]. We consider the cases per day, thus, we consider the time unit given by per-day. Among these model parameters, we consider the demographic parameters of birth and the natural death rate that are obtained using estimations while the rest of the parameter values have been obtained through model fitting to the data. The stochastic model (3) is used to obtain the parameter estimations. We consider the infected cases of the Kingdom of Saudi Arabia of coronavirus infection for the period starting March–July 2020 [38] is considered. The total population is considered to be $N\left(0\right)=34,813,871$ in Saudi Arabia in 2020 [39]. The average lifespan in the Kingdom of Saudi Arabia is given by $1/74.87$ [40]. The birth rate is determined from the equation $\lambda =d\ast N\left(0\right)$ which is $\Lambda =1273.94$ per day. The initial conditions are $S\left(0\right)=34,811,870$, $E\left(0\right)=2000$, $I\left(0\right)=1$, $A\left(0\right)=R\left(0\right)=0$. The basic reproduction number obtained through the data fitting is approximate ${\mathcal{R}}_0=1.1367$. The graphics for the data fitting are given in Figure 1 and Figure 2, while the parameters definitions and their realistic values are given in

Table 1. Parameters used in data fitting.

Parameter	Description	Value	Source
$\Lambda$	Birth rate	1273.94	Estimated
d	Natural death rate	$\frac{1}{74.87\times 365}$	Estimated
${\beta }_1$	Contact rate among E and S	0.1321	Fitted
${\beta }_2$	Contact rate among I and S	0.5926	Fitted
${\beta }_3$	Contact rate among A and S	0.5465	Fitted
$\delta$	Incubation period	0.7526	Fitted
$\tau$	Incubation period	0.0746	Fitted
$d_1$	Natural mortality due to disease	0.3691	Fitted
${\gamma }_1$	Recovery from I	0.2588	Fitted
${\gamma }_2$	Recovery from A	0.4639	Fitted

. Moreover, in Figure 1 and Figure 2, the legend “Data” referred to the infected cases while the “Model simulation” is the solution of the model”. Figure 1 is the model versus data fitting in the absence of ${\sigma }_i=0$, for $i=1,2,\dots ,5$ while Figure 2 is the data fitting versus model using ${\sigma }_1={\sigma }_2={\sigma }_4={\sigma }_5=0.002$ and ${\sigma }_5=0.004$. Figure 3 describes the comparison of stochastic and deterministic models versus infected data. It is clear that the stochastic case provides good results for data compared to the deterministic case. It should be noted that in Figure 1 and Figure 2, the legend “Data” referred to the infected cases while the “Model simulation” is the solution of the model. The fitting with the stochastic case is better than the deterministic case.

6. Numerical Results

This section is divided into three subsections where we study the numerical scheme, the results, and finally in the last subsection, the discussion on the results achieved for Model (3). We first present a numerical algorithm for solving the model, then obtain graphical results with a detailed discussion.

6.1. Numerical Scheme

We give the following numerical scheme for the model (3) using the scheme adopted from [41] called the Milstein method to obtain the numerical results. It follows from [41] by writing our model (3) in the following form:

$$\begin{array}{cc}\hfill S_{i+1}& =S_i+\left[\Lambda -\left({\beta }_1{E}_i+{\beta }_2{I}_i+{\beta }_3{A}_i\right)\frac{S_i}{N}-d{S}_i\right]△t+{\sigma }_1{S}_i\sqrt{△t}{\varsigma }_{1,i}+\frac{{\sigma }_1^2}{2}{S}_i({\varsigma }_{1,i}^2-1)△t,\hfill \\ \hfill E_{i+1}& =E_i+\left[\left({\beta }_1{E}_i+{\beta }_2{I}_i+{\beta }_3{A}_i\right)\frac{S_i}{N}-(\delta +d)E_i\right]△t+{\sigma }_2{E}_i\sqrt{△t}{\varsigma }_{2,i}+\frac{{\sigma }_2^2}{2}{E}_i({\varsigma }_{2,i}^2-1)△t,\hfill \\ \hfill I_{i+1}& =I_i+\left[(1-\tau )\delta E_i-\left(d+d_1+{\gamma }_1\right)I_i\right]△t+{\xi }_3{I}_i\sqrt{△t}{\varsigma }_{3,i}+\frac{{\sigma }_3^2}{2}{I}_i({\varsigma }_{3,i}^2-1)△t,\hfill \\ \hfill A_{i+1}& =A_i+\left[\tau \delta E_i-\left(d+{\gamma }_2\right)A_i\right]△t+{\sigma }_4{A}_i\sqrt{△t}{\varsigma }_{4,i}+\frac{{\sigma }_4^2}{2}{A}_i({\varsigma }_{4,i}^2-1)△t,\hfill \\ \hfill R_{i+1}& =R_i+\left[{\gamma }_1{I}_i+{\gamma }_2{A}_i-d{R}_i\right]△t+{\sigma }_5{R}_i\sqrt{△t}{\varsigma }_{5,i}+\frac{{\sigma }_5^2}{2}{R}_i({\varsigma }_{5,i}^2-1)△t.\hfill \end{array}$$

(34)

In the above scheme, ${\varsigma }_{i,j}$, when $i=1,\dots ,5$, displays the independent Gaussian random variables with the given distribution $N(0,1),$ and where the step size is shown by $\Delta t$, and ${\sigma }_i>0$, for $i=1,\dots ,5$ represents the numerical values of the noise.

6.2. Results

We solve the stochastic model (3) using the values of the parameters given in Table 1 and obtain the required graphical results, see Figure 4, Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9. We use in these graphical results the values of the stochastic noises are ${\sigma }_1=0.01$, ${\sigma }_2=0.02$, ${\sigma }_3=0.04$, ${\sigma }_4=0.02$, ${\sigma }_5=0.02$. The subgraphs in Figure 4, Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9 represent the population of susceptible, exposed, infected with disease symptoms, infected with no disease symptoms, and the recovery individuals, respectively, in (a–e). Figure 4 represents the simulation of the stochastic model versus the deterministic model.

Figure 5 represents the simulation of the model for different values of ${\beta }_1$. Reducing the contact between exposed and susceptible decreases the number of infective cases, also, one can see a better reduction in future cases if we restrict healthy individuals by following WHO recommendations.

Figure 6 gives the simulation of the stochastic model for many values of the parameter ${\beta }_2$. We can obtain a better decrease if the contact among healthy and symptomatic individuals is decreased, so ${\beta }_2$ can decrease best the number of infected cases caused by the symptomatic infected individuals.

Similarly, we can see from Figure 7 the decrease in the parameter ${\beta }_3$ decreases the number of infected people due to the asymptomatic infection.

The impact of the parameters $\delta$ and $\tau$ is shown in Figure 8 and Figure 9.

6.3. Discussion

The graphical results obtained in this paper can be helpful regarding disease elimination in the country. The comparison of the stochastic and the deterministic model is shown in Figure 4. In particular, in Figure 4e, it can be seen that the number of recovered cases increases in the stochastic case as compared to the deterministic case. Increasing the testing of people and identifications of the symptomatic and asymptomatic individuals by quarantine and educating them about the disease can best minimize the future infected cases in the population. It is useful if individuals in a certain country increase COVID-19 testing upon identifying the asymptomatic and symptomatic individuals. The asymptomatic individuals that do not show visible disease symptoms should be isolated in order to decrease the future risk of infected cases. The symptomatic individuals should be quarantined, possibly at home, and can also be restricted from visiting other places which can possibly increase the population of infected cases in society. It can be seen from the graphical results that the model behaves well and leads to the equilibrium point.

7. Conclusions

In this work, we presented the analysis of the COVID-19 model using stochastic differential equations. We first formulated the model in ordinary differential equations and then we extended it to stochastic differential equations. We studied with care the stochastic COVID-19 mathematical model and presented their mathematical as well as numerical results in detail. Initially, we formulated the model by taking into account the assumptions of the transmission coefficients, such as exposed, symptomatic, and asymptomatic with their interaction with healthy people. It is well-known that these interactions are possible methods of increasing the cases of COVID-19 in a population. With these assumptions, we extended the model into stochastic differential equations. We studied the analysis of the model for the deterministic case and presented the related mathematical results for it. We proved the existence and uniqueness of the deterministic model and found the existence and uniqueness of the model. We studied the existence of the endemic equilibria and found that the model has a unique endemic equilibrium. We found that the deterministic model is globally asymptotically stable when the basic reproduction number ${\mathcal{R}}_0\le 1$.

Further, we studied the stochastic model and presented their unique global positive solution. We carried out the results of extinction for the stochastic model. The extinction results for the stochastic model have been provided.

We used the infected cases of the Kingdom of Saudi Arabia for the period March–July 2020 and parameterized the model. Using the nonlinear square curve fitting method, we obtained the data fitting to the model and presented it graphically. The results of the stochastic case for different ${\sigma }_i$ values have good fitting as compared to the deterministic case.

The obtained parameters from the least square curve fit have been used to obtain the numerical results. The basic reproduction number computed for the given set of parameters is approximately, ${\mathcal{R}}_0\approx 1.1367$.

We also compared the data fitting to the model using deterministic and stochastic and showed that stochastic fitting is better than the deterministic case for this work, and present a comparative graph of both solutions. We gave a comparison of graphical results for the case of deterministic and stochastic solutions. Some important parameters that can possibly best decrease future cases have been shown graphically. Among these parameters, the contact between susceptible and exposed, susceptible and symptomatic, and susceptible and asymptomatic, can best decrease future cases, if the recommendations of the World Health Organization (WHO) are properly followed.

According to the results of our simulations, one can see that our results are in line with the WHO recommendations and can be useful for disease elimination in the country of Saudi Arabia, see for more details [42,43].