# Stationary Distribution and Extinction of a Stochastic SIQR Model with Saturated Incidence Rate

Authors: Qiuhua Zhang; Kai Zhou
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

In this paper, we consider a stochastic SIQR epidemic model with saturated incidence rate. By constructing a proper Lyapunov function, we obtain the existence and uniqueness of positive solution for this SIQR model. Furthermore, we study the dynamical properties of this stochastic SIQR model; that is, (i) we establish the sufficient condition for the existence of ergodic stationary distribution of the model; (ii) we obtain the extinction of the disease under some conditions. At last, numerical simulations are introduced to illustrate our theoretical results.


---

## Body

## 1. Introduction

Due to the importance of disease transmission and control, epidemic models are attracting more and more attention from scholars. Since the pioneer work of Kermack and McKendrick (see [1]), where they described a directly transmitted viral or bacterial agent in a closed population and divided the total population into three subparts as susceptibles (S), infectives (I), and recovers (R), mathematical models have proven their adequacy to describe and analyze the propagation and control of infectious diseases. From then on, various epidemic models have been developed in recent decades, such as SIR models, SIS models, and SEIR models; see [2–6].In modern society, due to the close interpersonal communication, the harm of disease transmission is more and more significant, and various measures have been taken to control the spread of disease. Among them, isolation measure was widely used by public health institutions and governments to control the spread of disease. In reality, it has been proved that the isolation measure is an effective method to suppress disease such as leprosy, plague, cholera, smallpox, tuberculosis, measles, rabies, foot and mouth, and ebola, which has ravaged West Africa in recent years [7]. Therefore, from the point of view of mathematic, the quarantine subpopulation (Q) was introduced into the epidemic models, and SIQ, SIQS, SIQR models have been proposed in the literature [8–10]. Afterwards Hethcote et al. [11] analyzed six SIQS and SIQR models with bilinear, standard, or quarantine-adjusted incidence and found that a Hopf bifurcation might occur in the SIQR model with quarantine-adjusted incidence rate. Recently, Nirwani et al. [12] considered a SIQR epidemic model with saturated incidence rate as follows:(1)dStdt=Λ-βStIt1+αIt-μSt,dItdt=βStIt1+αIt-δ+γ+μ+μ1It,dQtdt=δIt-κ+μ+μ2Qt,dRtdt=γIt+κQt-μRt,where they found that the equilibrium points are locally asymptotically stable under certain conditions.Recently, some researchers have considered the impact of environmental noise on the dynamics of epidemic models [13–16]. In reality, the environmental noise may involve the variations of factors such as climate, habitats, health habits, medical quality, and so on, which may affect the natural birth rate, death rate, and so on [17]. Especially, for human diseases, the nature of epidemic growth and spread is inherently stochastic, as the contacts of person-to-person are unpredictable [18]. Therefore, stochastic models may be a more accurate method to describe the spread of infectious diseases in humans (see, e.g., [19–26]). There are various approaches to include random perturbation in epidemic models, and among them, parameter perturbation method is a frequently adopted way. In this paper, we will adopt the method of Cai et al. [27]; that is, we replace -μdt in (1) with -μ~dt=-μdt+σidBi(t)(i=1,2,3,4). Therefore, we propose the following stochastic SIQR epidemic model with saturated incidence:(2)dSt=Λ-βStIt1+αIt-μStdt+σ1SdB1t,dIt=βStIt1+αIt-δ+γ+μ+μ1Itdt+σ2IdB2t,dQt=δIt-κ+μ+μ2Qtdt+σ3QdB3t,dRt=γIt+κQt-μRtdt+σ4RdB4t,where all parameters are positive. For the biological meaning of each coefficient, we briefly mention below: Λ is the recruitment rate of the population, β is the average number of adequate contacts of a person per unit time, α measures the appropriate preventions taken by infectious for epidemic control, and μ is the natural death rate of the population, while μ1,μ2 are the extra disease-related death rate in compartments I and Q, respectively. γ,κ denote the recover rates from group I,Q to R, respectively, and δ denotes the removal rate from I. Additionally, B1(t), B2(t), B3(t), and B4(t) are standard one-dimensional independent Brownian motions, and σi2>0 are the intensities of the white noises, i=1,2,3,4.This paper is organized as follows. In Section2, we give some preliminaries and notations which will be used in the following analysis. In Section 3, we prove that there is a unique global positive solution of system (2). In Section 4, we show the main result that there is an ergodic stationary distribution of system (2) under certain conditions. In Section 5, we establish sufficient conditions for extinction of the disease. Some numerical simulations are given to illustrate the analytical result in Section 6. Finally, we provide a brief summary of the main results and put forward some research prospects.


## 2. Preliminaries and Some Notations

Throughout this paper, letΩ,F,P be a complete probability space with a filtration Ftt≥0, which satisfies the usual conditions (i.e., it is increasing and right continuous while F0 contains all P-null stes) and Bi(t)(i=1,2,3,4) are defined on this complete probability space. Denote R+d=x∈Rd:xi>0,1<i<d and a∨b≔max⁡a,b, a∧b≔min⁡a,b.For thed-dimensional stochastic differential equation(3)dXt=fXt,tdt+gXt,tdBtfort≥0,with initial value X(0)=X0∈Rd, while B(t) is a d-dimensional standard Brownian motion defined on the complete probability space Ω,F,Ftt≥0,P, we define the differential operator L of (3) as follows:(4)L=∂∂t+∑i=1dfiX,t∂∂Xi+12∑i,j=1dgTX,tgX,tij∂2∂Xi∂Xj.For V(X,t)∈C2,1Rd×0,∞;R+, which represents the family of nonnegative functions defined on Rd×[0,∞) such that they are continuously twice differentiable in X and once in t, by acting L on this function V, we have(5)LVX,t=VtX,t+VXX,tfX,t+12tracegTX,tVXXX,tgX,t.By the Itô’s formula, if X(t)∈Rd, then(6)dVXt,t=LVXt,tdt+VXXt,tgXt,tdBt.At last, we give some theory about the stationary distribution. Firstly, we recall the definitions of homogeneous Markov process and stationary distribution.Definition 1 (see [28]).
The transition probability functionP(s,t,x,A) is said to be time-homogeneous (and the corresponding Markov process is called time-homogeneous) if the function P(s,t,x,A) is independent of s, where 0≤s≤t,x∈Rd, and A∈B and B denotes the σ-algebra of Borel sets in Rd.Definition 2 (see [29]).
DenotePγ the corresponding probability distribution of an initial distribution γ, which describes the initial state of (3) at t=0. Suppose that the distribution of X(t) with initial distribution γ converges in some sense to a distribution π=πγ, i.e.,(7)limt→∞⁡PγXt∈F=πFfor all measurable F, then we say that (3) has a stationary distribution π(·).LetX(t) be a homogeneous Markov process in Ed and satisfy the following stochastic differential equation:(8)dXt=bXdt+∑r=1kgrXdBrt.The diffusion matrix is denoted as A(x)=aijx, where aij(x)=∑r=1kgri(x)grj(x).Lemma 3 (see [28]).
The Markov processX(t) has a unique ergodic stationary distribution π(·) if there exists a bounded domain D⊂Ed with regular boundary Γ and
A 1: in the open domain D and some neighborhood thereof, the smallest eigenvalue of the diffusion matrix A(x) is bounded away from zero;
A 2: if x∈Ed∖D, the mean time τ which a path issuing from x reaches the set D is finite, and supx∈KExτ<∞ for every compact subset K∈Ed.
Then(9)PxlimT→∞⁡1T∫0TfXtdt=∫Edfxπdx=1for all x∈Ed, where f(·) is a function integrable with respect to the measure π.Remark 4.
As pointed in [24], to validate A1, we need to prove
A ~ 1: there is a positive number M such that ∑i,j=1daij(x)ξiξj≥Mξ2,x∈D,ξ∈Rd; and to validate A2, we need to prove
A ~ 2: there exists a nonnegative C2- function V such that LV is negative for any Ed∖D. Thus, we will verify the conditions A~1 and A~2 to obtain the existence of stationary distribution for Markov process X(t).


## 3. Existence and Uniqueness of the Positive Solution

In this section, we will give the existence and uniqueness of positive solutions for system (2), which is the basis of studying the long-time behavior of the model.Theorem 5.
For any initial valueX0=(S0,I0,Q0,R0)∈R+4, there is a unique positive solution X(t)=St,It,Qt,Rt of system (2) on t≥0 and the solution will maintain in R+4 with probability one.Proof.
The proof is motivated by the work Mao et al. [30]. Since the coefficients of system (2) satisfy the local Lipschitz condition, then for any initial value S0,I0,Q0,R0∈R+4, there exists a unique local solution St,It,Qt,Rt on [0,τe), where τe is the explosion time [31]. To verify that this solution is global, we only need to prove that τe=∞ a.s. For this purpose, let n0>0 be large enough such that every component of X0 lay within the interval [1/n0,n0]. For each integer n>n0, define the following stopping time:(10)τn=inf⁡t∈0,τe:min⁡St,It,Qt,Rt≤1normax⁡St,It,Qt,Rt≥n.For empty set Φ, we define infΦ=∞.
It is obvious thatτn is increasing with n→∞. Let τ∞=lim⁡supn→∞⁡τn, then τ∞≤τe a.s. In the following, we only need to prove that τ∞=∞ a.s. If this assertion is not valid, then there is a pair of constants T>0 and ϵ∈(0,1) such that Pτ∞≤T>ϵ. Hence, there exists an integer n1≥n0 such that(11)Pτn≤T≥ϵ,n≥n1. Define a C2- function V:R+4→R+ as follows:(12)VS,I,Q,R=S-1-ln⁡S+I-1-ln⁡I+Q-1-ln⁡Q+R-1-ln⁡R.Utilizing the Itô’s formula, we get that(13)dVS,I,Q,R=LVS,I,Q,Rdt+σ1S-1dB1t+σ2I-1dB2t+σ3Q-1dB3t+σ4R-1dB4t,where(14)LV=1-1SΛ-βSI1+αI-μS+1-1IβSI1+αI-δ+γ+μ+μ1I+1-1Q×δI-κ+μ+μ2Q+1-1RγI+κQ-μR+σ12+σ22+σ32+σ422=Λ+4μ+δ+γ+κ+μ1+μ2-μS+I+Q+R-μ1I-μ2Q+βI1+αI-βS1+αI-ΛS-δIQ-γIR-κQR+σ12+σ22+σ32+σ422<Λ+4μ+δ+γ+κ+βα+μ1+μ2+σ12+σ22+σ32+σ422≔G,where G is a positive constant which is independent of S,I,Q,R, and t. Therefore,(15)dVS,I,Q,R≤Gdt+σ1S-1dB1t+σ2I-1dB2t+σ3Q-1dB3t+σ4R-1dB4t.Integrating both sides of (15) from 0 to τn∧T, then (16)VSτn∧T,Iτn∧T,Qτn∧T,Rτn∧T-VS0,I0,Q0,R0≤Gτn∧T+∫0τn∧Tσ1S-1dB1t+σ2I-1dB2t+σ3Q-1dB3t+σ4R-1dB4tFurther, taking the expectation for the above inequality and due to the property of Brownian motion, we obtain(17)EVSτn∧T,Iτn∧T,Qτn∧T,Rτn∧T≤VS0,I0,Q0,R0+GT.Let Ωn=τn≤T for n≥n1. By (11) we have P(Ωn)≥ϵ. From the definition of the stopping time, we can deduce that, for every ω∈Ωn, S(τn,ω) or I(τn,ω) or Q(τn,ω) or R(τn,ω) equals either n or 1/n. Hence, it follows from (17) that (18)VS0,I0,Q0,R0+GT≥E1ΩnωVSτn,ω,Iτn,ω,Qτn,ω,Rτn,ω=PΩnEVSτn,ω,Iτn,ω,Qτn,ω,Rτn,ω≥ϵn-1-ln⁡n∧1n-1+ln⁡n,where 1Ωn represents the indicator function of Ωn. Letting n→∞ in above inequality yields a contradiction that(19)∞>VS0,I0,Q0,R0+GT=∞.Therefore, τ∞=∞ a.s. and this completes the proof.


## 4. Stationary Distribution and Ergodicity

In the field of infectious diseases, the persistence and extinction of diseases are an important research issue. In this section, we will show that system (2) has an ergodic stationary distribution under certain condition.Theorem 6.
Assume that(20)R0S≔βΛμ+σ12/2γ+δ+μ+μ1+σ22/2>1,then system (2) has a unique stationary distribution π(·) and it has the ergodic property.Proof.
By Theorem5, we have obtained that, for any initial value (S0,I0,Q0,R0)∈R+4, there is a unique global solution (S,I,Q,R)∈R+4.
The diffusion matrix of system (2) is given by(21)A=σ12S20000σ22I20000σ32Q20000σ42R2.Choose M=min(S,I,Q,R)∈D¯σ⊂Ed⁡σ12S2,σ22I2,σ32Q2,σ42R2, where Ed≕R+4 and Dσ is a bounded domain in Ed, which does not contain the origin. For ξ=(ξ1,ξ2,ξ3,ξ4)∈R4, we have that(22)∑i,j=14aijS,I,Q,Rξiξj=σ12S2ξ12+σ22I2ξ22+σ32Q2ξ32+σ42R2ξ42≥Mξ2,S,I,Q,R∈D¯σ.Then the condition A~1 holds.
Denote(23)η=μ+σ1221/2βΛγ+δ+μ+μ1+σ22/21/2-μ+σ1221/2.From the assumption R0S>1, we know that η is positive. Construct a function V^:R+4→R in the following form: (24)V^S,I,Q,R=M-ln⁡S-cln⁡I+1ρ+1S+I+Q+Rρ+1-ln⁡S-ln⁡Q-ln⁡R+S+I+Q+R≔MV1+V2+V3+V4+V5+V6,where(25)c=μ+σ12/2γ+δ+μ+μ1+σ22/2,0<ρ<2μσ12∨σ22∨σ32∨σ42and M is positive constant which satisfying the following condition:(26)-2Mη+Λ+B+3μ+κ+μ2+σ12+σ32+σ422≤-2, where B is determined later. It is evident that(27)lim⁡infk→∞,S,I,Q,R∈R+4∖UkV^S,I,Q,R=+∞,where Uk=(1/k,k)×(1/k,k)×(1/k,k)×(1/k,k). In addition, V^(S,I,Q,R) is a continuous function. Hence V^(S,I,Q,R) has a minimum point (S¯0,I¯0,Q¯0,R¯0) in the interior of R+4. Then we define a nonnegative C2-function V:R+4→R+ as follows:(28)VS,I,Q,R=V^S,I,Q,R-V^S¯0,I¯0,Q¯0,R¯0.Making use of Itô’s formula, we have(29)LV1=-ΛS+cβS1+αI+βI1+αI+μ+σ122+cγ+δ+μ+μ1+σ222.Utilizing the inequality a+b≥2ab, a,b>0, yields(30)LV1≤-2cβΛ1+αI1/2+βI1+αI+μ+σ122+cγ+δ+μ+μ1+σ222=-2βΛμ+σ12/2γ+δ+μ+μ1+σ22/21+αI1/2+2μ+σ122+βI1+αI=-2μ+σ1221/2βΛγ+δ+μ+μ1+σ22/21+αI1/2-μ+σ1221/2+βI1+αI.Meanwhile,(31)LV2=S+I+Q+RρΛ-μS-μ+μ1I-μ+μ2Q-μR+12ρS+I+Q+Rρ-1×σ12S2+σ22I2+σ32Q2+σ42R2≤S+I+Q+RρΛ-μS+I+Q+R+12ρS+I+Q+Rρ+1σ12∨σ22∨σ32∨σ42=ΛS+I+Q+Rρ-μ-12ρσ12∨σ22∨σ32∨σ42S+I+Q+Rρ+1≤B-12μ-12ρσ12∨σ22∨σ32∨σ42S+I+Q+Rρ+1≤B-12μ-12ρσ12∨σ22∨σ32∨σ42Sρ+1+Iρ+1+Qρ+1+Rρ+1,where(32)B=supS,I,Q,R∈R+4⁡ΛS+I+Q+Rρ-12μ-12ρσ12∨σ22∨σ32∨σ42S+I+Q+Rρ+1<∞.We can also obtain that(33)LV3=-ΛS+βI1+αI+μ+σ122,(34)LV4=-δIQ+κ+μ+μ2+σ322,(35)LV5=-γIR-κQR+μ+σ422and(36)LV6=Λ-μS+I+Q+R-μ1I-μ2Q≤Λ.Hence by (30), (31), (33), (34), (35), and (36), we can deduce that(37)LV≤-2Mμ+σ1221/2βΛγ+δ+μ+μ1+σ22/21+αI1/2-μ+σ1221/2+M+1βI1+αI-12μ-12ρσ12∨σ22∨σ32∨σ42Sρ+1+Iρ+1+Qρ+1+Rρ+1-ΛS-δIQ-γIR-κQR+Λ+B+3μ+κ+μ2+σ12+σ32+σ422. Next, we will construct a compact subset D such that the condition A2 in Lemma 3 holds. Thus, we define the following bounded subset:(38)D=ϵ1≤S≤1ϵ1,ϵ1≤I≤1ϵ1,ϵ2≤Q≤1ϵ2,ϵ2≤R≤1ϵ2,where ϵ1,ϵ2>0 are small enough and satisfy the following conditions: (39)-2Mη~+M+1βϵ1+Λ+B+3μ+κ+μ2+σ12+σ32+σ422≤-1,(40)-Λϵ1+E≤-1,(41)-δϵ1+E≤-1,(42)-γϵ1+E≤-1,(43)-12μ-12ρσ12∨σ22∨σ32∨σ421ϵ1ρ+1+E≤-1,(44)-12μ-12ρσ12∨σ22∨σ32∨σ421ϵ2ρ+1+E≤-1, where(45)η~≔μ+σ1221/2βΛγ+δ+μ+μ1+σ22/21+αϵ11/2-μ+σ1221/2.Note that, for sufficiently small ϵ1, (39) is valid due to inequality (26).
It can been easily seen thatR+4∖D is divided into eight domains, which can be denoted as (46)D1=S,I,Q,R∈R+4,0<S<ϵ1,D2=S,I,Q,R∈R+4,0<I<ϵ1,D3=S,I,Q,R∈R+4,S≥ϵ1,I≥ϵ1,0<Q<ϵ2,D4=S,I,Q,R∈R+4,I≥ϵ1,0<R<ϵ2,D5=S,I,Q,R∈R+4,S>1ϵ1,D6=S,I,Q,R∈R+4,I>1ϵ1,D7=S,I,Q,R∈R+4,Q>1ϵ2,D8=S,I,Q,R∈R+4,R>1ϵ2.Consequently, we consider the condition A2 from the following eight cases.
Case  1. If (S,I,Q,R)∈D1, from (37), we can get that(47)LV≤-ΛS+2Mμ+σ122+M+1βα+Λ+B+3μ+κ+μ2+σ12+σ32+σ422=-ΛS+E≤-Λϵ1+E,where(48)E=2Mμ+σ122+M+1βα+Λ+B+3μ+κ+μ2+σ12+σ32+σ422.Then it infers from (40) that(49)LV≤-1.
Case  2. If (S,I,Q,R)∈D2, we obtain that(50)LV≤-2Mμ+σ1221/2βΛγ+δ+μ+μ1+σ22/21+αI1/2-μ+σ1221/2+M+1βI+Λ+B+3μ+κ+μ2+σ12+σ32+σ422≤-2Mη~+M+1βϵ1+Λ+B+3μ+κ+μ2+σ12+σ32+σ422.In view of (39), we have that(51)LV≤-1.
Case  3. If (S,I,Q,R)∈D3, we have (52)LV≤-δIQ+2Mμ+σ122+M+1βα+Λ+B+3μ+κ+μ2+σ12+σ32+σ422≤-δϵ1ϵ2+E.Choosing(53)ϵ2=ϵ12,and combing with (41), we obtain that(54)LV≤-δϵ1+E≤-1.
Case  4. If (S,I,Q,R)∈D4, using the similar technique in Case 3 and (42), we have (55)LV≤-γIR+2Mμ+σ122+M+1βα+Λ+B+3μ+κ+μ2+σ12+σ32+σ422≤-γϵ1ϵ2+E=-γϵ1+E≤-1.
Case  5. If (S,I,Q,R)∈D5, we can derive that (56)LV≤-12μ-12ρσ12∨σ22∨σ32∨σ42Sρ+1+2Mμ+σ122+M+1βα+Λ+B+3μ+κ+μ2+σ12+σ32+σ422≤-12μ-12ρσ12∨σ22∨σ32∨σ421ϵ1ρ+1+E.By (43), we can conclude that LV≤-1 on D5.
Case  6. If (S,I,Q,R)∈D6, using the analogous estimation in Case 5 and (43), we have(57)LV≤-12μ-12ρσ12∨σ22∨σ32∨σ42Iρ+1+2Mμ+σ122+M+1βα+Λ+B+3μ+κ+μ2+σ12+σ32+σ422≤-12μ-12ρσ12∨σ22∨σ32∨σ421ϵ1ρ+1+E≤-1.
Case  7. If (S,I,Q,R)∈D7, it infers that(58)LV≤-12μ-12ρσ12∨σ22∨σ32∨σ42Qρ+1+2Mμ+σ122+M+1βα+Λ+B+3μ+κ+μ2+σ12+σ32+σ422≤-12μ-12ρσ12∨σ22∨σ32∨σ421ϵ2ρ+1+E.By (44), we can conclude that LV≤-1 on D7.
Case  8. At last, if (S,I,Q,R)∈D8, by (44), we have(59)LV≤-12μ-12ρσ12∨σ22∨σ32∨σ42Rρ+1+2Mμ+σ122+M+1βα+Λ+B+3μ+κ+μ2+σ12+σ32+σ422≤-12μ-12ρσ12∨σ22∨σ32∨σ421ϵ2ρ+1+E≤-1.According to the discussion of the above eight cases, we obtain that(60)LV≤-1∀S,I,Q,R∈R+4∖D.Therefore, condition A~2 is satisfied, and we can obtain that system (2) is ergodic and has a unique stationary distribution. The proof is complete.Remark 7.
Theorem6 gives the sufficient condition of the existence of stationary distribution for the solution of (2). In epidemiology, this means that the stochastic model (2) oscillates around the endemic equilibrium of the corresponding deterministic model and implies the persistence of the disease a.s. under some condition.


## 5. Extinction of the Disease

In this section, we will establish the sufficient conditions for the extinction of the disease in the system (2), which is another important issue in epidemiology. For simplicity, we define(61)xt=1t∫0txudu.Lemma 8.
LetSt,It,Qt,Rt be the solution of system (2) with any initial value (S0,I0,Q0,R0)∈R+4. Then(62)limt→∞⁡St+It+Qt+Rtt=0,a.s.Moreover, if μ>σ12∨σ22∨σ32∨σ42/2, then(63)limt→∞⁡1t∫0tSudB1u=0,limt→∞⁡1t∫0tIudB2u=0,a.s.,(64)limt→∞⁡1t∫0tQudB3u=0,limt→∞⁡1t∫0tRudB4u=0,a.s.Lemma8 can be proved by using the same methods as that in Lemmas 2.1 and 2.2 of Zhao and Jiang [32], so we omit it here.Define a constant as follows:(65)R^0S≔βΛ/μ-σ22/2δ+γ+μ+μ1=R0-σ222δ+γ+μ+μ1,where R0=βΛ/μ(δ+γ+μ+μ1) is the threshold value of the corresponding deterministic system for (2). Therefore, we obtain the following main result in this section.Theorem 9.
Assume thatμ>σ12∨σ22∨σ32∨σ42/2 and (S(t),I(t),Q(t),R(t)) is the solution of system (2) with any initial value (S0,I0,Q0,R0)∈R+4. If R^0S<1, then(66)limt→∞⁡ln⁡Itt≤δ+γ+μ+μ1R^0S-1a.s.,and(67)limt→∞⁡St=βΛμ,limt→∞⁡Qt=0,limt→∞⁡Rt=0a.s.Proof.
By system (2), we have(68)dSt+It+Qt+Rt=Λ-μSt+It+Qt+Rt-μ1It-μ2Qtdt+σ1StdB1t+σ2ItdB2t+σ3QtdB3t+σ4RtdB4t. Integrating the both sides of (68) yields (69)St-S0t+It-I0t+Qt-Q0t+Rt-R0t=Λ-μSt-μ+μ1It-μ+μ2Qt-μRt+σ1t∫0tSudB1u+σ2t∫0tIudB2u+σ3t∫0tQudB3u+σ4t∫0tRudB4u,then(70)St=Λμ-μ+μ1μIt-μ+μ2μQt-Rt+φt,where φ(t) is denoted as below: (71)φt=-1μSt-S0t+It-I0t+Qt-Q0t+Rt-R0t-σ1t∫0tSudB1u-σ2t∫0tIudB2u-σ3t∫0tQudB3u-σ4t∫0tRudB4u.Due to (62)-(64), we have(72)limt→∞⁡φt=0a.s.By utilizing Itô’s formula to the I component of system (2), we have(73)dln⁡I=βS1+αI-δ+γ+μ+μ1-12σ22dt+σ2dB2t≤βS-δ+γ+μ+μ1-12σ22dt+σ2dB2t.Integrating both sides of the above inequality from 0 to t and dividing t, we obtain(74)ln⁡It-ln⁡I0t≤βSt-δ+γ+μ+μ1-σ222+σ2B2tt.Substituting (70) into (74) leads to(75)ln⁡Itt≤βΛμ-δ+γ+μ+μ1-σ222+βφt+σ2B2tt+ln⁡I0t=δ+γ+μ+μ1R^0S-1+βφt+σ2B2tt+ln⁡I0t.Due to the fact that limt→∞(B2(t)/t)=0 a.s. and (72), by taking the superior limit on the both sides of (75), if R^0S<1, then(76)lim⁡supt→∞ln⁡Itt≤δ+γ+μ+μ1R^0S-1<0a.s.,which infers that(77)limt→∞⁡It=0a.s.Hence, from (77) and the last two equations of (2), we have(78)limt→∞⁡Qt=0,limt→∞⁡Rt=0a.s.At last, by (72), (77), and (78) and from (70), we obtain(79)limt→∞⁡St=βΛμa.s.This completes the proof.


## 6. Numerical Simulations

In this section, we give some numerical simulations in order to illustrate the theoretical results. For numerical simulation, by using Milstein’s method [33], we obtain the following discretization transformation of system (2):(80)Sj+1=Sj+Λ-βSjIj1+αIj-μSjΔt+σ1SjΔtε1,j+σ122Sjε1,j2-1Δt,Ij+1=Ij+βSjIj1+αIj-δ+γ+μ+μ1IjΔt+σ2IjΔtε2,j+σ222Ijε2,j2-1Δt,Qj+1=Qj+δIj-κ+μ+μ2QjΔt+σ3QjΔtε3,j+σ322Qjε3,j2-1Δt,,Rj+1=Rj+γIj+κQj-μRjΔt+σ4RjΔtε4,j+σ422Rjε4,j2-1Δt,where Δt>0, σi2>0(i=1,2,3,4) are the intensities of the white noise, εi,j(i=1,2,3,4,j=1,2,…,n) are independent Gaussian random variables which follow N(0,1). In order to verify the persistence and extinction results, we choose the parameters in system (80) by β=0.3, μ=μ1=0.1, μ2=γ=δ=0.2, κ=0.1, α=0.5, σ1=σ3=0.05, and σ4=0.1.Firstly, we start our numerical simulation withΛ=0.3, σ2=0.1 and starting from the initial point (2,0.5,0.5,1). Direct calculation leads to R0S≔βΛ/(μ+σ12/2)(γ+δ+μ+μ1+σ22/2)≈1.4692>1. Thus, from Theorem 6, we can conclude that system (2) has a unique ergodic stationary distribution which leads to the stochastically persistence of the disease. Figure 1 confirms this.Figure 1
The left column shows the solutions of stochastic system (2) under the noise intensities σ1=σ3=0.05,σ2=0.4,σ4=0.1. The right column displays the histogram of the probability density function of S(t),I(t),Q(t),R(t).Next, we increase the stochastic intensity ofI to σ2=0.4 and choose Λ=0.22. In this case, μ>σ12∨σ22∨σ32∨σ42/2 and R^0S≔R0-σ22/2(δ+γ+μ+μ1)≈0.9667<1<R0. According to Theorem 9, we obtain that system (2) has disease dies out with probability 1, and Figure 2 illustrates this.Figure 2
Numerical simulation for system (2) with σ1=σ3=0.05,σ2=σ4=0.1.


## 7. Conclusion and Discussion

In this paper, we consider a stochastic SIQR epidemic model with saturated incidence rate. By constructing proper Lyapunov function, we establish the existence and uniqueness of solutions of system (2). Furthermore, we obtain the existence of ergodic stationary distribution of the solution to (2) when R0S>1. Meanwhile, we establish sufficient conditions for the extinction of the disease as well, that is, μ>σ12∨σ22∨σ32∨σ42/2 and R^0S<1.It is worth further studying the following two topics. On the one hand, noting thatR0S≠R^0S, it is of great interest that whether can we promote the sufficient conditions for the extinction of the disease and get the threshold value for the extinction and persistence of the disease? On the other hand, we can propose SIQR epidemic model with some other types of random perturbations, such as telephone noise which can cause the population system to switch among different environmental regimes (see, e.g., [34–37]). In the literatures, it is common to use the continuous-time Markov chain with values in a finite state space to describe the random switching of environmental regimes. We will carry out these studies in the future.



---

*Source: 3575410-2019-06-24.xml*