# Dynamic Analysis of a Stochastic SEQIR Model and Application in the COVID-19 Pandemic

Authors: Shijie Liu; Maoxing Liu
Journal: Discrete Dynamics in Nature and Society
 (2021)

---

## Abstract

In this study, a deterministic SEQIR model with standard incidence and the corresponding stochastic epidemic model are explored. In the deterministic model, the reproduction number is given, and the local asymptotic stability of the equilibria is proved. When the reproduction number is less than unity, the disease-free equilibrium is locally asymptotically stable, whereas the endemic equilibrium is locally asymptotically stable in the case of a reproduction number greater than unity. A stochastic expansion based on a deterministic model is studied to explore the uncertainty of the spread of infectious diseases. Using the Lyapunov function method, the existence and uniqueness of a global positive solution are considered. Then, the extinction conditions of the epidemic and its asymptotic property around the endemic equilibrium are obtained. To demonstrate the application of this model, a case study based on COVID-19 epidemic data from France, Italy, and the UK is presented, together with numerical simulations using given parameters.


---

## Body

## 1. Introduction

At the end of 2019, COVID-19 was reported in Wuhan, China [1], setting the beginning of a global epidemic. Now, there are more than 200 countries trapped in the epidemic disaster, which has seriously affected the development of the society. It was found that COVID-19 was caused by the SARS-CoV-2 coronavirus [2, 3]. This coronavirus mainly spreads in three ways: (i) direct transmission, which refers to the infection caused by patient’s sneezing, coughing, and talking droplets, and the direct inhalation of exhaled gas at a short distance; (ii) aerosol transmission, which refers to droplet mixing in the air to form an aerosol that can lead to infection after inhalation; (iii) contact transmission, which refers to the droplet deposited on the surface of an object, contact with contaminated hands, and then contact with the oral cavity, nasal cavity, eyes, and other mucous membranes, resulting in infection.In the absence of a specific therapeutic method, governments can only do their best to hinder the spread of the epidemic through measures such as population lockdown, reasonable allocation of medical resources, and calls for public health (mask wearing, avoidance of places with many people, meetings with reduced attendance, etc.) [4, 5]. These measures can effectively prevent epidemics and delay the arrival of epidemic peaks [6]. For now, through the aforementioned measures and vaccine production, the epidemic situation has been controlled in many countries; however, there are still some countries with severe epidemic situations.Mathematical modeling is a suitable approach to describe the spread of diseases and predict their epidemic trend [7]. Based on the spread characteristics of diseases, many epidemic models, e.g., SIR, SEIR, and SIQR [8–12], were proposed. Using mathematical models to describe infectious diseases can help to better understand the spread trends of diseases and propose reasonable control methods and suggestions. Novel coronavirus pneumonia models were reported by scholars who have been involved in this task since the outbreak of the epidemic. Tang et al. constructed an SEIR model in [13]; their results indicate that the number of infected people and the propagation of infectious diseases can be reduced when governments take certain measures. In [14], the author predicted that the spread of infectious disease would be costly. After that, the governments of the UK and USA established strict movement restrictions. Many models, related suggestions, and opinions on COVID-19 have been provided [15–21].Although there are many mathematical models, few considered the inherent uncertainties in these models. Manski and Molinari [22] showed that there is an uncertainty in the epidemic transmission process. Therefore, assuming that stochastic disturbance applied to a model is suitable, some stochastic epidemic models were proposed [18, 23–26]. In [18], the authors studied an SLIRD with the deterministic and stochastic models to analyze actual COVID-19-infected cases in Spain. In the numerical simulations, they fitted data from Spain and obtained the daily basic reproduction number. Thus, they predicted the epidemic situation for the next few months. In [23], by extending the SIR epidemiological model, the authors established a random propagation model for additional modeling of whether an individual is far from a crowded area. The results showed that the propagation of the epidemic in Japan would gradually slow down by reducing the time spent in crowded zones. Anwarud et al. analyzed a stochastic SIQ model and obtained sufficient conditions for disease extinction, disease existence, and stationary distribution in [24]. The numerical simulations were divided into two parts: fitting data from Khyber Pakhtunkhwa, Pakistan, and verification of the theoretical results for the given parameters. In [25], the authors studied an SEQIR model with Markovian switching, established a random threshold of disease extinction and persistence, and used the data from Indian states to confirm their conclusions. In [26], an epidemic system was considered through an additive fractional white noise. It was indicated that the epidemic under fractional random settings may be more suitable for modeling than deterministic modeling. Based on the above studies, an SEQIR model with standard incidence is analyzed in this study, together with deterministic and stochastic models. A case study of COVID-19 epidemic data from France, Italy, and the UK is presented, and numerical simulations using given parameters are reported.The remainder of this paper is organized as follows. In the next section, deterministic and stochastic SEQIR models with standard incidence are described. The dynamic analysis of the deterministic model is presented in Section3, and an analysis of the stochastic system is presented in Section 4. Numerical simulations and application to the COVID-19 epidemic are described in Section 5. The paper ends with a discussion in Section 6.


## 2. Model Formulation

The entire population is divided into five classes according to its states, namely, susceptibleSt, exposed Et, quarantined Qt, hospitalized infected It, and recovered Rt. Then, Nt=St+Et+Qt+It+Rt. In [19], Mandal et al. constructed an SEQIR model with optimized control and provided advice to Indian cities. Based on their work, Brahim et al. considered a model with Markovian switching and generalized incidence in [25]. Interestingly, they considered that people who are exposed and asymptomatic are infectious. The standard incidence rate is more suitable for infectious disease models with a large population. Based on these results, we consider an SEQIR model with standard incidence as follows:(1)dStdt=A−βStEtNt+b1Qt−dSt,dEtdt=βStEtNt−b2+α+dEt,dQtdt=b2Et−b1+c+dQt,dItdt=αEt+cQt−η+d+δIt,dRtdt=ηIt−dRt.The constant input of the susceptible population isA; β represents the transmission rate of the disease; α is the rate of the exposed class removed from the infected class; b2 is the portion of the exposed class moving to the quarantined class; b1 and c indicate the rate at which the quarantined people become susceptible and infected types, respectively; the recovery rate of hospitalized infected people is η; d is the natural mortality rate; and δ is the diseased death rate.A popular technique to incorporate stochasticity into a deterministic model is parametric perturbation. Following this method, the transmission rate,β, is replaced by β⟶β+σdBt. Therefore, we have a stochastic model defined as follows:(2)dSt=A−βStEtNt+b1Qt−dStdt−σStEtNtdBt,dEt=βStEtNt−b2+α+dEtdt+σStEtNtdBt,dQt=b2Et−b1+c+dQtdt,dIt=αEt+cQt−η+d+δItdt,dRt=ηIt−dRtdt.Note thatBt is a standard independent Brownian motion, and the intensity of the white noise is σ. An initial condition should be satisfied when analyzing systems (1) and (2):(3)S0≥0,E0≥0,Q0≥0,I0≥0,R0≥0.In the following analysis, we declare that the standard Brownian motion in the system is defined on a complete probability spaceΩ,ℱ,P with a filtration ℱtt∈R+ satisfying the usual conditions (i.e., it is right continuous and increasing while ℱ0 contains all P-null sets).


## 3. Analysis of the Deterministic Model

In this section, we analyze system (1) without stochastic perturbation. In this part, some of the basic dynamics are analyzed, including the basic reproduction number and stability of the equilibria.


### 3.1. The Basic Reproduction Number

An important parameter in infectious diseases is the basic reproduction number which can determine the outbreak or extinction of the disease. Therefore, we give it by the next-generation matrix method [27]. Because the infection classes involve E, Q, and I classes, we consider the three classes:(4)dEtdt=βStEtNt−b2+α+dEt,dQtdt=b2Et−b1+c+dQt,dItdt=αEt+cQt−η+d+δIt.Infection subsystem (4) has transmission matrix F and transition matrix V, where(5)F=β00000000,V=b2+α+d00−b2b1+c+d0−α−cη+d+δ.So, according to this method, we can get that the spectral radius ofFV−1 is R0=β/b2+α+d.


### 3.2. The Equilibria

System (1) has two equilibria, the disease-free equilibrium E0=A/d,0,0,0,0 and the endemic equilibrium E∗=S∗,E∗,Q∗,I∗,R∗, respectively, where(6)S∗=b2+α+dN∗β,E∗=−N∗dα+b2+d+Aβb1+c+dβα+db1+c+d+b2c+d,Q∗=Aβ−N∗db2+α+db2b1+b2+α+c+dd+b1+cα+b2cβ,I∗=αb1+c+d+b2cb2η+d+δQ∗,R∗=ηdI∗.Here,N∗=S∗+E∗+Q∗+I∗+R∗. It is worth noting that the endemic equilibrium E∗ exists when R0>1.


### 3.3. The Stability of Two Equilibria

In this part, we study the local asymptotic stability of two equilibria.Theorem 1.
The disease-free equilibriumE0 is locally asymptotically stable if R0<1.Proof.
The Jacobian matrix at the disease-free equilibrium of system (1) is given by(7)J|E0=−d−βb1000β−b2+α+d0000b2−b1+c+d000αc−η+d+δ0000η−d.
The characteristic equation of system (1) at its disease-free equilibrium is given by(8)λ+dλ+b2+α+d1−R0λ+b1+c+dλ+η+d+δλ+d=0.
It is easy to see that the characteristic values of the Jacobian matrix are negative if and only ifR0<1. Hence, the disease-free equilibrium of system (1) is locally asymptotically stable if R0<1.Theorem 2.
The endemic equilibriumE∗ is locally asymptotically stable for R0>1.Proof.
The Jacobian matrix at the endemic equilibrium of system (1) is given by(9)J|E∗=−βE∗N∗−d−βS∗N∗b100βE∗N∗βS∗N∗−b2+α+d0000b2−b1+c+d000αc−η+d+δ0000η−d.
The characteristic equation of system (1) at its endemic equilibrium is given by(10)λ+dλ+η+d+δλ3+a1λ2+a2λ+a3=0,where a1=b1+c+2d+βE∗/N∗, a2=βE∗/N∗+db1+c+d+β2E∗S∗/N∗2, and a3=βE∗/N∗b1α+d+c+db2+α+d. It is clear that the characteristic polynomial has two negative eigenvalues, and we study the last cubic polynomial. By calculation, we can verify a1,a2,a3>0 and a1a2−a3>0. Hence, according to the Routh–Hurwitz criterion, we can know that the endemic equilibrium of system (1) is locally asymptotically stable for R0>1.


## 3.1. The Basic Reproduction Number

An important parameter in infectious diseases is the basic reproduction number which can determine the outbreak or extinction of the disease. Therefore, we give it by the next-generation matrix method [27]. Because the infection classes involve E, Q, and I classes, we consider the three classes:(4)dEtdt=βStEtNt−b2+α+dEt,dQtdt=b2Et−b1+c+dQt,dItdt=αEt+cQt−η+d+δIt.Infection subsystem (4) has transmission matrix F and transition matrix V, where(5)F=β00000000,V=b2+α+d00−b2b1+c+d0−α−cη+d+δ.So, according to this method, we can get that the spectral radius ofFV−1 is R0=β/b2+α+d.


## 3.2. The Equilibria

System (1) has two equilibria, the disease-free equilibrium E0=A/d,0,0,0,0 and the endemic equilibrium E∗=S∗,E∗,Q∗,I∗,R∗, respectively, where(6)S∗=b2+α+dN∗β,E∗=−N∗dα+b2+d+Aβb1+c+dβα+db1+c+d+b2c+d,Q∗=Aβ−N∗db2+α+db2b1+b2+α+c+dd+b1+cα+b2cβ,I∗=αb1+c+d+b2cb2η+d+δQ∗,R∗=ηdI∗.Here,N∗=S∗+E∗+Q∗+I∗+R∗. It is worth noting that the endemic equilibrium E∗ exists when R0>1.


## 3.3. The Stability of Two Equilibria

In this part, we study the local asymptotic stability of two equilibria.Theorem 1.
The disease-free equilibriumE0 is locally asymptotically stable if R0<1.Proof.
The Jacobian matrix at the disease-free equilibrium of system (1) is given by(7)J|E0=−d−βb1000β−b2+α+d0000b2−b1+c+d000αc−η+d+δ0000η−d.
The characteristic equation of system (1) at its disease-free equilibrium is given by(8)λ+dλ+b2+α+d1−R0λ+b1+c+dλ+η+d+δλ+d=0.
It is easy to see that the characteristic values of the Jacobian matrix are negative if and only ifR0<1. Hence, the disease-free equilibrium of system (1) is locally asymptotically stable if R0<1.Theorem 2.
The endemic equilibriumE∗ is locally asymptotically stable for R0>1.Proof.
The Jacobian matrix at the endemic equilibrium of system (1) is given by(9)J|E∗=−βE∗N∗−d−βS∗N∗b100βE∗N∗βS∗N∗−b2+α+d0000b2−b1+c+d000αc−η+d+δ0000η−d.
The characteristic equation of system (1) at its endemic equilibrium is given by(10)λ+dλ+η+d+δλ3+a1λ2+a2λ+a3=0,where a1=b1+c+2d+βE∗/N∗, a2=βE∗/N∗+db1+c+d+β2E∗S∗/N∗2, and a3=βE∗/N∗b1α+d+c+db2+α+d. It is clear that the characteristic polynomial has two negative eigenvalues, and we study the last cubic polynomial. By calculation, we can verify a1,a2,a3>0 and a1a2−a3>0. Hence, according to the Routh–Hurwitz criterion, we can know that the endemic equilibrium of system (1) is locally asymptotically stable for R0>1.


## 4. Analysis of the Stochastic Model

Theorem 3.
For any initial valuesS0,E0,Q0,I0,R0∈R+5, there exists a unique global solution St,Et,Qt,It,Rt∈R+5 of system (2) for all t≥0, and the solution will remain in R+5 with probability 1, i.e., St,Et,Qt,It,Rt∈R+5 for all t≥0 almost surely.Proof.
Since the coefficients of system (2) are locally Lipschitz continuous on R+, for any initial value S0,E0,Q0,I0,R0∈R+5, there exists a unique local solution St,Et,Qt,It,Rt on 0,τe, where τe represents the explosion time [28]. If we can prove τe=∞ a.s, then the solution is global. To this end, let k0≥1 be sufficiently large such that S0,E0,Q0,I0, and R0 all lie in the interval 1/k0,k0. For each integer k≥k0, we can define the stopping time [28] through τk=inft∈0,τe:St∉1/k,korEt∉1/k,korQt∉1/k,korIt∉1/k,korRt∉1/k,k, where, throughout this paper, we set inf∅=∞ (as usual ∅ represents the empty set). Obviously, τk is increasing as k⟶∞. Let τ∞=limk⟶∞τk, whence τ∞≤τe a.s. If we can show τ∞=+∞, then τe=+∞ and St,Et,Qt,It,Rt∈R+5 a.s for all t≥0. That is to say, to complete the proof, all we need is to show that τ∞=+∞ a.s. If this assertion is false, then there exists a pair of constants T>0 and ε1∈0,1 such that Pτ∞≤T>ε1. There is an integer k1≥k0 such that Pτk≤T≥ε1 for all k≥k1. Define a Lyapunov function V:R+⟶R+ by VS,E,Q,I,R=S−1−lnS+E−1−lnE+Q−1−lnQ+I−1−lnI+R−1−lnR.
Sinceu−1−lnu≥0 for any u>0, the function V is nonnegative. Making use of Ito^’s formula (see [28]) to V, we get(11)LV=1−1SA−βSEN+b1Q−dS+σ2E22N2+1−1EβSEN−b2+α+dE+σ2S22N2+1−1Qb2E−b1+c+dQ+1−1IαE+cQ−η+d+δI+1−1RηI−dR=A−βSEN+b1Q−dS−AS+βEN−b1QS+d+σ2E22N2+βSEN−b2+α+dE−βSN+b2+α+d+σ2S22N2+b2E−b1+c+dQ−b2EQ+b1+c+d+αE+cQ−η+d+δI−αEI−cQI+η+d+δ+ηI−dR−ηIR+d≤A+β+α+b1+b2+η+δ+c+5d+σ2≜K,dV=LVdt+σE−SNdBt≤Kdt+σE−SNdBt,where K is a positive constant. Integrating from 0 to τk∧T=minτk,T and then taking the expectation on both sides, we obtain(12)EVSτk∧T,Eτk∧T,Qτk∧T,Iτk∧T,Rτk∧T≤VS0,E0,Q0,I0,R0+KEτk∧T.
Thus, we have(13)EVSτk∧T,Eτk∧T,Qτk∧T,Iτk∧T,Rτk∧T≤VS0,E0,Q0,I0,R0+KT.
SetΩk=τk≤T for k≥k1, and we get PΩk≥ε1. Note that, for every ω∈Ωk, there exists that Sτk,ω,Eτk,ω,Qτk,ω,Iτk,ω,Rτk,ω equals either kor1/k. Therefore, VSτk,ω,Eτk,ω,Qτk,ω,Iτk,ω,Rτk,ω is no less than either k−1−lnkor1/k−1−ln1/k=1/k−1+lnk. Therefore, we can see(14)VSτk,ω,Eτk,ω,Qτk,ω,Iτk,ω,Rτk,ω≥k−1−lnk∧1k−1+lnk.
Then, we get(15)VS0,E0,Q0,I0,R0+KT≥E1ΩkωVSτk,ω,Eτk,ω,Qτk,ω,Iτk,ω,Rτk,ω≥ε1k−1−lnk∧1k−1+lnk,where 1Ωkω is the indicator function of Ωk. Letting k⟶∞,(16)+∞>VS0,E0,Q0,I0,R0+KT=+∞.is a contradiction, and thus, we can get τ∞=+∞a.s., which implies that St,Et,Qt,It,Rt∈R+5 for all t∈R+. This completes the proof.Remark 1.
DefineΓ=St,Et,Qt,It,Rt∈R+5:A/d+δ≤St+Et+Qt+It+Rt≤A/d. By system (2), we can get d/dtS+E+Q+I+R=A−dS+E+Q+I+R−δI. Next, we have A−d+δS+E+Q+I+R≤A−dS+E+Q+I+R−δI≤A−dS+E+Q+I+R. Therefore, we can obtain the solution of system (2) in the region Γ. Then, we can give some general conclusions and results.Theorem 4.
LetSt,Et,Qt,It,Rt be the positive solution of system (2) with the initial value S0,E0,Q0,I0,R0∈R+5. If the assumption β2/2σ2<b2+α+d holds, then the disease will go to extinction exponentially with probability one, i.e.,(17)limt⟶∞Et=limt⟶∞Qt=limt⟶∞It=limt⟶∞Rt=0,a.s.Proof.
Applying Ito^’s formula, we have(18)dlnEt=βSN−b2+α+d−σ2S22N2dt+σSNdBt.
Integrating (18) from 0 to t, we get(19)lnEt−lnE0=∫0tβSuNu−b2+α+d−σ2S2u2N2udu+∫0tσSuNudBu=∫0t−σ2S2u2N2u+βSuNudu−b2+α+dt+Mt,where Mt=∫0tσSu/NudBu is a continuous local martingale [28] whose quadratic variation is(20)M,Mt=σ2∫0tS2uN2udu.
Applying Theorem 7.4 in [28], p44, we can get(21)Psup0≤t≤kMt−c2M,Mt>2clnk≤1k2,where 0<c<1 and k>1 is a random integer. Using the Borel–Cantelli lemma [28], we obtain that, for almost all ω∈Ω, there exists a random integer k0ω such that, for all k>k0,(22)sup0≤t≤kMt−c2M,Mt≤2clnk.
That is to say,(23)Mt≤c2σ2∫0tS2uN2udu+2clnk,for any0≤t≤k.
Taking it to (19), we obtain(24)lnEt−lnE0≤∫0t−121−cσ2S2uN2u+βSuNudu−b2+α+dt+2clnk.
Observing the integrand function, we can do the following:(25)−121−cσ2S2uN2u+βSuNu≤β221−cσ2.
Consequently, we get(26)lnEt−lnE0≤β221−cσ2t−b2+α+dt+2clnk.
Therefore, fork−1≤t≤k, one can see that(27)lnEtt≤lnE0t+β221−cσ2−b2+α+d+2clnkk−1.
Lettingk⟶+∞, t⟶+∞, lnE0/t⟶0, and lnk/k−1⟶0. We can obtain(28)limsupt⟶+∞lnEtt≤β221−cσ2−b2+α+da.s.
Lettingc⟶0 gives(29)limsupt⟶+∞lnEtt≤β22σ2−b2+α+da.s,which shows that limt⟶+∞Et=0a.s.
On the other hand, there exists a constantT such that Ew≤ε,t>T. Applying the third equation of system (2), we get(30)dQt=b2E−b1+c+dQdt≤b2ε−b1+c+dQdt.
Using the comparison theorem, we have(31)limsupt⟶+∞Qt≤b2εb1+c+da.s.
Lettingε⟶0, we get limt⟶+∞Qt=0a.s. Similarly, in the case that limt⟶+∞Et=0a.s. and limt⟶+∞Qt=0a.s., we can get(32)limt⟶+∞It=0a.s.,limt⟶+∞Rt=0a.s.
This completes the proof.Theorem 5.
Assume thatR0>1, m1=d−b1/2−σ2d+δb2+α+2dE∗/βA, m2=b2/2+α/2+d−b1/2, m3=c/2+d−b2/2, m4=η/2+d+δ−α/2−c/2, m5=d−η/2, M=minm1,m2,m3,m4,m5, and Λ=σ2d+δb2+α+2dE∗/βAS∗2+AS∗δb2+α+2d/dd+δ are all positive constants. Then, for any initial value S0,E0,Q0,I0,R0∈R+5, the solution of system (2) has the property(33)limsupt⟶∞1tE∫0tS−S∗2+E−E∗2+Q−Q∗2+I−I∗2+R−R∗2du≤ΛM.Proof.
System (1) has the endemic equilibrium, and we can get(34)A=βS∗E∗N∗−b1Q∗+dS∗,βS∗E∗N∗=b2+α+dE∗,b2E∗=b1+c+dQ∗,αE∗+cQ∗=η+d+δI∗,ηI∗=dR∗.
We construct a function(35)V=12S−S∗+E−E∗2+χE−E∗−E∗lnEE∗+12Q−Q∗2+12I−I∗2+12R−R∗2=V1+χV2+V3+V4+V5,where χ is a constant. By Ito^’s formula, we obtain(36)LV=LV1+χLV2+LV3+LV4+LV5.
And we get(37)LV1=S−S∗+E−E∗A−dS+b1Q−b2+α+dE=−dS−S∗2−b2+α+dE−E∗2−b2+α+2dS−S∗E−E∗+b1S−S∗Q−Q∗+b1E−E∗Q−Q∗≤−d−b12S−S∗2−b2+α+d−b12E−E∗2+b1Q−Q∗2−b2+α+2dS−S∗E−E∗,(38)LV2=1−E∗EβSEN−b2+α+dE+σ2S22N2E∗≤E−E∗βSd+δA−βS∗d+δA+βS∗δA+σ2S2d+δ22A2E∗≤βd+δAS−S∗E−E∗+σ2d+δ2S−S∗+S∗22A2E∗+βS∗δd≤βd+δAS−S∗E−E∗+σ2d+δ2S−S∗2A2E∗+σ2d+δ2S∗2A2E∗+βS∗δd,(39)LV3=Q−Q∗b2E−b1+c+dQ=b2E−E∗Q−Q∗−b1+c+dQ−Q∗2≤−b1+c+d−b22Q−Q∗2+b22E−E∗2,(40)LV4=I−I∗αE+cQ−η+d+δI=αE−E∗I−I∗+cQ−Q∗I−I∗−η+d+δI−I∗2≤−η+d+δ−α2−c2I−I∗2+α2E−E∗2+c2Q−Q∗2,(41)LV5=R−R∗ηI−dR=ηI−I∗R−R∗−dR−R∗2≤−d−η2R−R∗2+η2I−I∗2.
Taking (37)–(41) into (36), one can obtain that(42)LV≤−d−b12S−S∗2−b2+α+d−b12E−E∗2+b1Q−Q∗2−b2+α+2dS−S∗E−E∗+χβd+δAS−S∗E−E∗+σ2d+δ2S−S∗2A2E∗+σ2d+δ2S∗2A2E∗+βS∗δd+b22E−E∗2−b1+c+d−b22Q−Q∗2−η+d+δ−α2−c2I−I∗2+α2E−E∗2+c2Q−Q∗2−d−η2R−R∗2+η2I−I∗2≤−d−b12−σ2d+δb2+α+2dE∗βAS−S∗2−b22+α2+d−b12E−E∗2−c2+d−b22Q−Q∗2−η2+d+δ−α2−c2I−I∗2−d−η2R−R∗2+σ2d+δb2+α+2dE∗βAS∗2+AS∗δb2+α+2ddd+δ≜−m1S−S∗2−m2E−E∗2−m3Q−Q∗2−m4R−R∗2−m5R−R∗2+Λ,where χ=A/βd+δb2+α+2d, mii=1,2,3,4,5, and Λ are defined in Theorem 5.(43)dV=LVdt+σSE−E∗NdBt.
Integrating (39) from 0 to t and taking the expectation on both sides we have(44)0≤EV≤EVS0,E0,Q0,I0,R0+E∫0t−m1S−S∗2−m2E−E∗2−m3Q−Q∗2−m4I−I∗2−m5R−R∗2+Λdu.
Then, lettingM=minm1,m2,m3,m4,m5, dividing t on both sides of (44), and taking the superior limit, we get(45)limsupt⟶∞1tE∫0tS−S∗2+E−E∗2+Q−Q∗2+I−I∗2+R−R∗2du≤ΛM.


## 5. Application of the Model

### 5.1. Numerical Simulations

The numerical simulations are divided into two aspects. We first set values of the parameters to verify the theorem results obtained in this study. Specifically, we setA=50, β=0.8, b1=0.3, b2=0.5, α=0.3, c=0.2, η=0.18, d=0.2, and δ=0.25. Using different initial values for St,Et,Qt,It,Rt such as 1500,100,50,10,0, 1000,200,60,20,0, and 500,10,5,1,0, we numerically solve system (1). Figure 1 shows that when R0<1, the epidemic will go to extinction, and the solutions with different initial values will converge to the disease-free equilibrium expressed as E0=250,0,0,0,0, which coincides with the conclusion of Theorem 1. For β=1.2 and the other parameters set as before, it is verified that when R0>1, the epidemic will continue to spread, and the solutions will converge to the endemic equilibrium expressed as E∗=199.1,12.95,9.248,8.82,8.2 in Figure 1.Figure 1
The local asymptotic stability of the equilibria of system (1). The blue line is the case of R0>1, and the red dash line is the case of R0<1.If we setβ=1.1 and σ=0.8, the other parameters remain unchanged, and the initial value is given by 500,10,5,1,0, while these parameters satisfy β2/2σ2=0.945<b2+α+d=1, the condition of Theorem 4 is satisfied as shown in Figure 2, and in this case, the epidemic will die out. When β=1.1, σ=0.01, and σ=0.04, the other parameters and initial value are the same as before, and we obtain the asymptotic property of system (2) and the solutions with stochastic oscillation near the endemic equilibrium of system (1), as shown in Figure 3.Figure 2
The epidemic extinction with stochastic perturbation.Figure 3
The asymptotic solutions of the deterministic and stochastic system. The first row is the case ofσ=0.01, and the second row is the case of σ=0.04.


### 5.2. Case Study

In this section, we present some cases in which the proposed research model can be applied. After the WHO reported the first case on December 31, 2019, the COVID-19 pandemic has spread rapidly worldwide. To date, 220 countries have suffered from this pandemic. By May 30, 2021, there were 171,474,925 cases of global confirmed cases, and the cumulative mortality reached 3,565,330 cases. The infectious, noxious, and transmission speeds of the epidemic are shocking. Some countries are still trapped in the epidemic; here, we analyze data from France, Italy, and the UK. The data used in our study were extracted from Worldometers [29], an online freely available repository. In this study, we consider 95 days of time-series data of currently infected patients from three countries from February 24 to May 30, 2021. The source of the parameters is as follows. The population of Italy [30] is approximately 60,500,000, which we set it as the initial value of the total population. The value of A can be set on the basis of the number of births in Italy, 1427. The average life expectancy in Italy is 84.01 [31]. Given that 1/d is the average life expectancy, the value of d can be calculated as d=1/84.01×365≈3.2612×10−5. Through the same calculations, we derive some parameters for France and the UK. We also set the range of α, which represents the rate of exposure to infection according to [32]. Because the average recovery period in France is approximately 13.6 days [33], the parameter η is set to 1/13.6≈0.0735. According to [34], we take the average recovery period in Italy as 13.15 days; thus, in this case, we set the parameter η=1/13.15≈0.076. The recovery rate of the UK, β,b1,b2, and c are estimated by the least square method.We fit the parameters of the system by using the daily active COVID-19 cases in France, Italy, and the UK. The parameters resulting from fitting the actual active cases are shown in Table1; the results are shown in Figure 4. We consider hospitalized infected patients as active cases in this part. Note from Figure 4 that the simulations are correct, and the epidemic development of these three countries shows a decreasing trend. It is evident that the people infection will not stop shortly. The government should take measures to reduce the number of newly infected people and maintain a decreasing trend.Table 1
Description and values of parameters.
ParameterDefinitionFranceItalyThe UKSourceATotal input221814272356[30]βDisease transmission rate0.9260.5710.999EstimateddNatural death rate3.3616×10−53.2612×10−53.2929×10−5[31]ηThe recovery rate ofI0.073 50.0760.071 4[33, 34],estimatedb1The rate thatQ becomes S0.10.0120.1Estimatedb2The rate thatE becomes Q0.7180.335 40.773 4EstimatedαThe rate thatE becomes I0.19230.192 40.199[32]cThe rate thatQ becomes I0.630 80.80.1EstimatedδThe mortality rate ofI0.024 90.033 80.029 4[31]Figure 4
The fitting results of system (1). (a) France. (b) Italy. (c) The UK.
(a)(b)(c)According to the previous analysis, stochastic perturbation can change the solution of the system and the stage of disease development. The epidemic becomes extinct when the noise intensity is sufficiently large. However, if the noise intensity is small, the overall trend of disease development will not change. This situation can be seen in Figure5, which shows the solutions of the deterministic and stochastic models with σ=0.02, illustrating that the stochastic solution fluctuates around the deterministic solution. This leads to conclude that the government should take some measures to prevent new infections. Figure 6 shows the increase in b2, which represents the strengthening isolation, and α, which means that raising nucleic acid detection will effectively reduce the risk of disease transmission.Figure 5
The solutions of system (2) with σ=0.02. The first row is the case of France, the middle row is the case of Italy, and the last row is the case of the UK. All parameters are listed in Table 1.Figure 6
The relationship between parameters andR0. (a) β,b2, and R0. (b) β,α, and R0.
(a)(b)


## 5.1. Numerical Simulations

The numerical simulations are divided into two aspects. We first set values of the parameters to verify the theorem results obtained in this study. Specifically, we setA=50, β=0.8, b1=0.3, b2=0.5, α=0.3, c=0.2, η=0.18, d=0.2, and δ=0.25. Using different initial values for St,Et,Qt,It,Rt such as 1500,100,50,10,0, 1000,200,60,20,0, and 500,10,5,1,0, we numerically solve system (1). Figure 1 shows that when R0<1, the epidemic will go to extinction, and the solutions with different initial values will converge to the disease-free equilibrium expressed as E0=250,0,0,0,0, which coincides with the conclusion of Theorem 1. For β=1.2 and the other parameters set as before, it is verified that when R0>1, the epidemic will continue to spread, and the solutions will converge to the endemic equilibrium expressed as E∗=199.1,12.95,9.248,8.82,8.2 in Figure 1.Figure 1
The local asymptotic stability of the equilibria of system (1). The blue line is the case of R0>1, and the red dash line is the case of R0<1.If we setβ=1.1 and σ=0.8, the other parameters remain unchanged, and the initial value is given by 500,10,5,1,0, while these parameters satisfy β2/2σ2=0.945<b2+α+d=1, the condition of Theorem 4 is satisfied as shown in Figure 2, and in this case, the epidemic will die out. When β=1.1, σ=0.01, and σ=0.04, the other parameters and initial value are the same as before, and we obtain the asymptotic property of system (2) and the solutions with stochastic oscillation near the endemic equilibrium of system (1), as shown in Figure 3.Figure 2
The epidemic extinction with stochastic perturbation.Figure 3
The asymptotic solutions of the deterministic and stochastic system. The first row is the case ofσ=0.01, and the second row is the case of σ=0.04.


## 5.2. Case Study

In this section, we present some cases in which the proposed research model can be applied. After the WHO reported the first case on December 31, 2019, the COVID-19 pandemic has spread rapidly worldwide. To date, 220 countries have suffered from this pandemic. By May 30, 2021, there were 171,474,925 cases of global confirmed cases, and the cumulative mortality reached 3,565,330 cases. The infectious, noxious, and transmission speeds of the epidemic are shocking. Some countries are still trapped in the epidemic; here, we analyze data from France, Italy, and the UK. The data used in our study were extracted from Worldometers [29], an online freely available repository. In this study, we consider 95 days of time-series data of currently infected patients from three countries from February 24 to May 30, 2021. The source of the parameters is as follows. The population of Italy [30] is approximately 60,500,000, which we set it as the initial value of the total population. The value of A can be set on the basis of the number of births in Italy, 1427. The average life expectancy in Italy is 84.01 [31]. Given that 1/d is the average life expectancy, the value of d can be calculated as d=1/84.01×365≈3.2612×10−5. Through the same calculations, we derive some parameters for France and the UK. We also set the range of α, which represents the rate of exposure to infection according to [32]. Because the average recovery period in France is approximately 13.6 days [33], the parameter η is set to 1/13.6≈0.0735. According to [34], we take the average recovery period in Italy as 13.15 days; thus, in this case, we set the parameter η=1/13.15≈0.076. The recovery rate of the UK, β,b1,b2, and c are estimated by the least square method.We fit the parameters of the system by using the daily active COVID-19 cases in France, Italy, and the UK. The parameters resulting from fitting the actual active cases are shown in Table1; the results are shown in Figure 4. We consider hospitalized infected patients as active cases in this part. Note from Figure 4 that the simulations are correct, and the epidemic development of these three countries shows a decreasing trend. It is evident that the people infection will not stop shortly. The government should take measures to reduce the number of newly infected people and maintain a decreasing trend.Table 1
Description and values of parameters.
ParameterDefinitionFranceItalyThe UKSourceATotal input221814272356[30]βDisease transmission rate0.9260.5710.999EstimateddNatural death rate3.3616×10−53.2612×10−53.2929×10−5[31]ηThe recovery rate ofI0.073 50.0760.071 4[33, 34],estimatedb1The rate thatQ becomes S0.10.0120.1Estimatedb2The rate thatE becomes Q0.7180.335 40.773 4EstimatedαThe rate thatE becomes I0.19230.192 40.199[32]cThe rate thatQ becomes I0.630 80.80.1EstimatedδThe mortality rate ofI0.024 90.033 80.029 4[31]Figure 4
The fitting results of system (1). (a) France. (b) Italy. (c) The UK.
(a)(b)(c)According to the previous analysis, stochastic perturbation can change the solution of the system and the stage of disease development. The epidemic becomes extinct when the noise intensity is sufficiently large. However, if the noise intensity is small, the overall trend of disease development will not change. This situation can be seen in Figure5, which shows the solutions of the deterministic and stochastic models with σ=0.02, illustrating that the stochastic solution fluctuates around the deterministic solution. This leads to conclude that the government should take some measures to prevent new infections. Figure 6 shows the increase in b2, which represents the strengthening isolation, and α, which means that raising nucleic acid detection will effectively reduce the risk of disease transmission.Figure 5
The solutions of system (2) with σ=0.02. The first row is the case of France, the middle row is the case of Italy, and the last row is the case of the UK. All parameters are listed in Table 1.Figure 6
The relationship between parameters andR0. (a) β,b2, and R0. (b) β,α, and R0.
(a)(b)


## 6. Conclusions

In 2020, new coronavirus pneumonia seriously affected people’s daily life and hindered social development. More than 200 countries worldwide suffered from this epidemic. Prevention of this coronavirus pneumonia has become the most important issue worldwide. Infectious disease models constitute an effective tool for analyzing and understanding the trend of disease spread; therefore, we built a mathematical model with standard incidence based on the spread of COVID-19. We obtained the basic reproduction number, which affects the spread of the disease, and proved that whenR0<1, the disease-free equilibrium is locally asymptotically stable, and when R0>1, the endemic equilibrium is locally asymptotically stable. Then, we incorporated stochastic perturbation in the deterministic model because there exists stochastic influence in the process of disease spread. Specifically, we obtained the conditions under which the epidemic goes to extinction and the asymptotic property near the endemic equilibrium of the deterministic model.Finally, we verified the results in two aspects. First, we set some parameters to verify the proposed theorem. Second, we fitted the active data cases using the least square method. In Figure4, we effectively fitted the active cases in France, Italy, and the UK. The active cases of these three countries all decreased. In the first part of the numerical simulations, we took into account that adding a stochastic disturbance will prevent the spread of the disease. Therefore, we considered stochastic disturbance in the model for comparison with the deterministic model. The parameters in Table 1 do not satisfy the condition of disease extinction, so the epidemic will not go to extinction. It is well known that white noise with sufficiently large intensity can lead to the extinction of the disease. Hence, we can control the development of the disease by increasing the noise intensity.Random perturbance can be considered as a human-issue disturbance. Therefore, governments must take measures to decrease the number of infections, e.g., population lockdown, increase of media publicity, rational allocation of medical resources, and appeal to the public to pay attention to personal hygiene. Some suggestions, such as increasing the isolation and raising nucleic acid detection, should be made. With these suggestions, the number of infections will be reduced. Currently, novel coronavirus pneumonia vaccines have been developed. It is believed that a reasonable combination of government control and vaccination strategies will eliminate the epidemic. Next, we will also consider the model with immunity to analyze the impact of immune intensity on the spread of infectious diseases.



---

*Source: 6125064-2021-12-08.xml*