# A Stochastic SIR Epidemic System with a Nonlinear Relapse

Authors: Ali El Myr; Abdelaziz Assadouq; Lahcen Omari; Adel Settati; Aadil Lahrouz
Journal: Discrete Dynamics in Nature and Society
 (2018)

---

## Abstract

We investigate the conditions that control the extinction and the existence of a unique stationary distribution of a nonlinear mathematical spread model with stochastic perturbations in a population of varying size with relapse. Numerical simulations are carried out to illustrate the theoretical results.


---

## Body

## 1. Introduction

In medicine, relapse is the return of a disease or the signs and symptoms of a disease after a period of improvement. Relapse also refers to returning to the use of an addictive substance or behavior, such as cigarette smoking [1]. For example, for human tuberculosis, incomplete treatment can lead to relapse, but relapse can also occur in patients who take a full course of treatment and are declared cured. Recently, considerable attention has been paid to model the relapse phenomenon. In [2] Tuder developed one of the first epidemic models with relapse in a constant population with bilinear incidence rate. Moreira and Wang [3] included a nonlinear incidence rate in the model. Van Den Driessche and Zou in [4] formulated a SIRI epidemic model as an integro-differential system with the fraction P(t) of recovered individuals remaining in the recovered class t units after the recovery. In [5] the displacement of recovered individuals to the infective class due to relapse is given by ηRt. Inspired by the works cited above and the fact that relapse is due to contact with infected, it is more reasonable to consider a bilinear relapse rate σβRtIt. We consider the following SIR compartmental model in a population of varying size with a bilinear relapse rate.(1)dSt=Λ-μ1St-βStItdt,dIt=-μ2+αIt+βStIt+σβRtItdt,dRt=-μ3Rt+αIt-σβRtItdt.In this model each letter refers to a compartment in which an individual can reside. Let St denote the number of members of a population susceptible to the disease at time t, It the number of infective members, and Rt the number of members who have been removed from the possibility of infection with permanent or temporary immunity. The parameters that occur in the model have the following meaning. Λ is the rate at which new individuals enter the population. α is the rate at which the infective individuals become recovered. σ∈[0,1] is the parameter that measure the intensity of the relapse. The positive constants μ1, μ2, and μ3 satisfying(2)μ1≤min⁡μ2,μ3,represent the natural death rate of susceptible, infected, and recovered individuals, respectively. Another addition in the modeling of population dynamics of diseases is the introduction of stochasticity into epidemic models. Many scholars have studied the effect of stochasticity on epidemic models [6–10]. For instance, to include stochastic demographic variability, Allen [6] studied SDEs for simple SIS and SIR epidemic models with constant population size that was derived from a continuous time Markov chain model. In [7, 10], the situation of a white noise stochastic perturbations around the endemic equilibrium state was considered. Lahrouz et al. in [11] formulated a stochastic version of the classical SIS epidemic model with varying population size. The authors studied the long time behavior of the stochastic system. They also gave conditions for extinction and persistence of the disease in the population. According to the value of the threshold Rσ=βΛ/μ(μ+λ+(1/2)σ2), they showed that if Rσ<1, the disease will die out from the population with the probability one and the disease will persist if Rσ>1. In the case of persistence, they proved the existence of a stationary distribution.Let(Ω,F,{Ft}t≥0,P) be a complete probability space with a filtration {Ft}t≥0 satisfying the usual conditions (i.e., it is right continuous and F0 contains all Prob-null sets). In this paper, we assume that fluctuations in the environment will manifest themselves mainly as fluctuations in the death rates of S, I, and R. Specifically, μ1dt is replaced with μ1dt+σ1dB1; that is, the rate is perturbed by Gaussian white noise. The rates μ2 and μ3 are similarly perturbed by an independent Gaussian white noises. Therefore, the corresponding stochastic system to (1) can be described by the Itô equation:(3)dSt=Λ-μ1St-βStItdt-σ1StdBt1,dIt=-μ2+αIt+βStIt+σβRtItdt-σ2ItdBt2,dRt=-μ3Rt+αIt-σβRtItdt-σ3RtdBt3,where Bt1, Bt2, and Bt3 are independent Brownian motions; σ1, σ2, σ3 represent the intensities of the white noises. By using the same method as in [10, 12], the existence and uniqueness of positive solution for system (3) hold with probability 1, if we start from any positive initial value (S0,I0,R0). The main concern of the present paper is to establish a sufficient condition for the extinction and the persistence of solutions of the system (3).


## 2. Preliminaries

Throughout the rest of this paper, we denote(4)R+3=x1,x2,x3∣xi>0,i=1,2,3.In general, consider the n-dimensional stochastic differential equation: (5)dXt=fXt,tdt+gXt,tdBt,X0=x0∈Rn,wheref:Rn→Rn, g:Rn→Rn×d, and B(t) denotes a d-dimensional standard Brownian motion defined on the underlying probability space. If A is a vector or matrix, its transpose is denoted by AT. The n×n matrix (6)Σx≜gxTgxis called the diffusion matrix. For the convenience of a later presentation, we introduce the generator L associated with (5) as follows. For any twice continuously differentiable V(x)∈C2(Rn)(7)LVx=∇Vxfx+12traceΣx∇2Vx,where ∇V, ∇2V denote the gradient, Hessian of V, respectively.The following theorem gives a criterion for positive recurrence in terms of Lyapunov function (see [13] Theorem 3.13, p. 1164, Theorem 4.3, p. 1168, and Theorem 4.4, p. 1169).Theorem 1.
The system (5) is positive recurrent if there is a bounded open subset Δ of Rn with a regular boundary, and the following holds.
(i) There exist someδ∈0,1 such that, for all x∈Δ, (8)δξ2≤ξTΣxξ≤δ-1ξ2,for anyξ∈Rn,(ii) There exists a nonnegative function V:Δc→R such that V is twice continuously differentiable and that for some θ>0(9)LV≤-θ,for anyx∈ΔcMoreover, the positive Markov process X(t) has a unique ergodic stationary distribution π. That is, if h is a function integrable with respect to the measure π, then (10)Plimt→∞⁡1t∫0thXsds=∫Rnhxπdx=1.


## 3. Extinction of the Disease

In this section, we followed the methods of Lahrouz et al. [11] to establish sufficient condition for the extinction of the disease. Before this, let us prepare two useful lemmas. In the following lemma, we show that the positive solutions to (3) have finite moments.Lemma 2.
LetNt=St+It+Rt be the total population size in system (3) for an initial positive value (S0,I0,R0). Then for any q>1 such that μ1-(1/2)(q-1)σ¯2>0, where σ¯=max⁡(σ1,σ2,σ3), we have (11)ilimsupt→+∞⁡ENtq≤Λμ1-1/2q-1σ¯2q,(12)iisupn⁡Emaxt∈n,n+1Ntq<+∞.Proof.
( i ) Let q>1 such that μ1-(1/2)(q-1)σ¯2>0 and let θ>0 where μ1-(1/2)(q-1)σ¯2-θ/q>0.
Applying Itô’s formula leads to(13)deθtNtq=LeθtNtqdt-qeθtNtq-1σ1StdB1+σ2ItdB2+σ3RtdB3.where(14)LeθtNtq=θeθtNtq+eθtqNtq-1Λ-μ1St-μ2It-μ3Rt+qq-12eθtNtq-2σ12St2+σ22It2+σ32Rt2≤qeθtΛNtq-1-μ1-12q-1σ¯2-θqNtq.Note that the function Λxq-1-(μ1-(1/2)(q-1)σ¯2-θ/q)xq is bounded in (0,∞) with (15)maxx>0⁡Λxq-1μ1-12q-1σ¯2xq=q-1q-1Λqqqμ1-1/2q-1σ¯2-θ/qq-1.integrating (13); then taking expectation on both sides and using (14) and (15), we obtain (16)ENtq≤N0qe-θt+q-1q-1Λqqq-1θμ1-1/2q-1σ¯2-θ/qq-11-e-θt.Therefore, for all θ∈0,qμ1-1/2q-1σ¯2, we have (17)limsupt→+∞⁡ENtq≤q-1q-1Λqqq-1θμ1-1/2q-1σ¯2-θ/qq-1and for the fact that (18)min0<θ<qμ1-1/2q-1σ¯2⁡q-1q-1Λqqq-1θμ1-1/2q-1σ¯2-θ/kq-1=Λμ1-1/2q-1σ¯2q.we have (i).
(ii) From (13) and (14) we have, for all t∈[n,n+1], (19)Ntq≤en-tθNnq+∫ntes-tθΛNsq-1-μ1-12q-1σ¯2-θqNsqds-∫ntqeθs-tNsq-1σ1SsdB1s+σ2IsdB2s+σ3RsdB3sds≤Nnq+∫ntΛNsq-1-μ1-12q-1σ¯2-θqNsqds-∫ntqeθt-sNsq-1σ1SsdB1s+σ2IsdB2s+σ3RsdB3sds.so (20)maxt∈n,n+1⁡Ntq≤Nnq+∫nn+1ΛNsq-1-μ1-12q-1σ¯2-θqNsqds+maxt∈n,n+1Mt,where (21)Mt=-∫ntqeθt-sNsq-1σ1SsdB1s+σ2IsdB2s+σ3RsdB3sdsIn view of (11) and the continuity of Nt, there exists C>0 such that (22)supt≥0⁡ENtq≤C.So, there exists a positive constant C1 such that (23)Emaxt∈n,n+1⁡Ntq≤C1+Emaxt∈n,n+1⁡MtApplying Itô’s formula leads to (24)dNtq=qNtq-1Λ-μ1St-μ2It-μ3Rtdt+qq-12Ntq-2σ12St2+σ22It2+σ32Rt2dt-qNtq-1σ1StdB1t+σ2ItdB2t+σ3RtdB3t.From(25)-∫ntqNsq-1σ1SsdB1s+σ2IsdB2s+σ3RsdB3sds=Ntq-Nnq-∫ntqNsq-1Λ-μ1Ss-μ2Is-μ3Rsds-∫ntqq-12Nsq-2σ12Ss2+σ22Is2+σ32Rs2ds,we have(26)Mt≤∫ntqNsq-1σ1SsdB1s+σ2IsdB2s+σ3RsdB3sds≤Ntq+Nnq+∫nn+1qNsq-1Λ-μ1Ss-μ2Is-μ3Rsds+∫nn+1qq-12Nsq-2σ12Ss2+σ22Is2+σ32Rs2ds.So, there exists a positive constant C2 such that(27)Mt≤Ntq+Nnq+C2∫nn+1Nsq-1+Nsqds.Letting q′>1 we have (28)Emaxt∈n,n+1⁡Mtq′≤q′q′-1q′maxt∈n,n+1⁡EMtq′≤C3ENtqq′+Nnqq′+∫nn+1Nsq-1q′+Nsq-1q′ds,where C3 is a positive constant independent of n, and the first inequality is derived from the maximal inequality for martingales and the second by (27) and Jensen’s inequality.
Choosingq′=1+γ where 0<γ<2μ1/σ12 we have (29)q′>1,qq′>1andμ1-12qq′-1σ¯2>0.So, by (22) we get Emaxt∈[n,n+1]Mtq′<+∞, and by (23) we have (30)Emaxt∈n,n+1⁡Ntq≤C1+Emaxt∈n,n+1⁡Mtq′q′<+∞.Under the conditions of the preceding lemma we have the following lemma.Lemma 3.
Let(St,It,Rt) be the solution of (3) with the positive initial condition (S0,I0,R0).
Assume that(1/2)σ¯2<μ1. Then (31)ilimt→∞⁡Ntt=0.iilimt→∞⁡1t∫0tSudBu=limt→∞⁡1t∫0tIudBu=limt→∞⁡1t∫0tRudBu=0.Proof.
( i ) In view of (12) and Markov’s inequality there exists a positive constant C such that (32)Psupn≤t≤n+1⁡Ntt>n-q′=Psupn≤t≤n+1⁡Ntqtq>n-qq′≤Cnq-qq′,where q verifies the conditions of Lemma 3 and q′ is a positive constant such that q′<1-1/q.
An application of Borel-Cantelli lemma yields to for almostω∈Ω and there is a random integer n0(ω) such that for n≥n0(33)supn≤t≤n+1⁡Ntt≤n-q′,hence for all t∈n,n+1 and n>n0(34)Ntt≤n-q′.Let n→∞, then t→∞ so (35)lim∞⁡Ntt=0.(ii) We shall prove that limt→∞(1/t)∫0tSudB(u)=0, and the other limits can be obtained in the same way. Put Xt=∫0tSudB(u). In view of the assumption (1/2)σ¯2<μ1 and (11), one can derive that there exists C>0 such that (36)supt≥0ENt2≤C.Therefore (37)EXt2=E∫0tSu2du≤Ct.Beside, let δ1,δ2 such that δ1>2 and δ1+2≤δ2<2δ1. Thanks to Doob’s martingale inequality, and using (37), we get(38)Psupk-1δ1≤t≤kδ1⁡Xt2≥kδ2≤EXkδ12kδ2≤Ckδ2-δ1,Since Σ(1/kδ2-δ1)<∞, the Borel-cantelli lemma implies that for almost ω∈Ω, there is k0(ω) such that for k≥k0(39)supk-1δ1≤t≤kδ1⁡Xt2<kδ2.Hence, for all t∈[(k-1)δ1,kδ1] and k≥k0(40)X2tt2≤kδ2k-12δ1.Since 2δ1>δ2, we have kδ2/(k-1)2δ1~1/k2δ1-δ2 and limk→∞(1/k2δ1-δ2)=0 leads to (41)limt→∞⁡X2tt2=0.That is(42)limt→∞⁡1t∫0tSudBu=0Now, we are in the position to establish the threshold criterion for disease extinction in terms of the positive number: (43)RT≜Λβμ1α+μ2+σ22/2.Theorem 4.
Let(S(t),I(t),R(t)) be the solution of system (3) with any initial value (S(0),I(0),R(0))∈R+3. if (1/2)σ¯2<μ1 and RT<1, then we have the following property: (44)limt→∞⁡It=limt→∞⁡Rt=0.Proof.
Applying Itô’s formula leads to(45)dlog⁡It=-μ2+α+σ222+βSt+σβRtdt-σ2dB2≤-μ2+α+σ222+βSt+Rt-σ2dB2≤-μ2+α+σ222+βNt-σ2dB2,and (46)dNt=Λ-μ1St-μ2It-μ3Rtdt-σ1StdB1-σ2ItdB2-σ3RtdB3≤Λ-μ1Nt-σ1StdB1-σ2ItdB2-σ3RtdB3,so (47)dlog⁡It+βμ1Nt≤βΛμ1-μ2+α+σ222-σ2dB2-σ1StdB1-σ2ItdB2-σ3RtdB3.Integrating (47) and dividing by t>0 we obtain (48)log⁡Itt-log⁡I0t+βμ1Nt-N0t≤μ2+α+σ222RT-1-σ2B2-σ1t∫0tSudB1u-σ2t∫0tIudB2u-σ3t∫0tRudB3u,and using Lemma 3 yields (49)limsupt→∞⁡log⁡Itt≤μ2+α+σ222RT-1.Hence by condition RT<1 we have limt→∞It=0 a.s.
( i i ) Let (50)ψt=exp⁡μ3+σ322t+σ3B3.We know that limt→∞(B3/t)=0. Then (51)limt→∞⁡log⁡ψtt=μ3+σ322≜h′.From the fact that (52)limsupt→∞⁡log⁡Itt≤μ2+α+σ222RT-1≜h,we deduce that for any ϵ>0 and for almost all ω∈Ω there exists T(ω) such that for all t≥T(53)It≤K1eh+ϵtandK2eh′-ϵt≤ψt≤K3eh′+ϵtwhere Ki,i=1,2,3 are almost surely finite positive random variables. By (53) we have (54)dRt≤-μ3Rt+αK1eh+ϵtdt-σ3RtdB3.Let ϕ be the solution of stochastic equation: (55)dϕt=-μ3ϕt+αK1eh+ϵtdt-σ3ϕtdB3,and we can write (56)ϕt=ϕ0ψt-1+αK1ϕt-1∫0tψseh+ϵsds.So from (53) we have (57)ϕt≤ϕ0K2e-h′-ϵt+αK2e-h′-ϵt∫0Tψseh+ϵsds+αK1K3K2e-h′-ϵ∫Tteh′+ϵseh+ϵsds=ϕ0K2+αK2∫0Tψseh+ϵsds-αK1K3K2h+h′+2ϵeh+h′+2ϵTe-h′-ϵt+αK1K3K2h+h′+2ϵeh+h′+2ϵteh+3ϵt.By comparison theorem for stochastic differential equations we have (58)Rt≤ϕ0K2+αK2∫0Tψseh+ϵsds-αK1K3K2h+h′+2ϵeh+h′+2ϵTe-h′-ϵt+αK1K3K2h+h′+2ϵeh+h′+2ϵteh+3ϵt.Since ϵ is arbitrary, we get (59)Rt≤ϕ0K2+αK2∫0Tψsehsds-αK1K3K2h+h′eh+h′Te-h′t+αK1K3K2h+h′eh+h′teht.From h′>0, RT<1 and (59) we obtain limt→∞Rt=0.


## 4. Stationary Distribution and Positive Recurrence

In many papers, the existence of stationary distribution needs the construction of suitable Lyapunov functions that are based on the positive equilibrium state of the deterministic system, which gives strong sufficient conditions [14]. In the following theorem, to investigate the existence of an asymptotically invariant distribution for the solution of model (3), we did not use equilibrium state to construct Lyapunov functions.Theorem 5.
The solutions(St,It,Rt) are positive recurrent and admit a unique ergodic stationary distribution, provided that the following conditions hold: (60)RT>1and (61)μi>σi22,i=1,2,3Proof.
We construct a positive functionW such that (62)supS,I,R∈R+3∖Δρ⁡LWS,I,R≤-1for some ρ>0, where Δρ=(1/ρ,ρ)×(1/ρ,ρ)×(1/ρ,ρ).
First of all we have(63)L-log⁡S=-ΛS+βI+μ1+σ122.Second, we have (64)L12N2=NΛ-μ1S+μ2I+μ3R+σ122S2+σ222I2+σ322R2≤ΛN-μ1-σ122S2-μ1-σ222I2-μ1-σ322R2.By using the elementary inequality S≤ϵS2+1/4ϵ2, where ϵ is a positive constant to be choosen later, we have (65)L12N2≤3Λ4ϵ-μ1-σ122-ϵΛS2-μ2-σ222-ϵΛI2-μ3-σ322-ϵΛR2.Now we compute (66)LSlog⁡I=Λlog⁡I-μ1Slog⁡I-βSIlog⁡I-Sμ2+α+σβRS+βS2-σ222S.By using the elementary inequalities -1/e≤Ilog⁡I≤I2 and SR≤(1/2)(S2+R2) we have (67)LSlog⁡I≤Λlog⁡I-μ1Slog⁡I+βe-μ2-α-σ222S+σβ2R2+σβ2+βS2.Next, we compute (68)L12log⁡I2=log⁡I-μ2+α+βS+σβR+σ2221-log⁡I=-σ222+μ2+αlog⁡I+βSlog⁡I+σβRlog⁡I+σ222and (69)L-log⁡R=μ3-αIR+σβI+σ322and (70)LRlog⁡I=log⁡I-μ3R+αI-σβIR+R-μ2+α+βS+σβR-σ222R≤-μ3Rlog⁡I+αI2+σβe-σ222R+β2S2+β2+σβR2.Now let (71)W1S,I,R=-log⁡S-log⁡R+14log⁡I2+ω¯4N2and (72)W2S,I,R=ω¯4N2+βμ1Slog⁡I+σβμ3Rlog⁡I+14log⁡I2,where ω¯ is a positive constant. It is easy to check that W1(S,I,R) has a minimum point (S^,I^,R^) in R+3, and by choosing ω¯ sufficiently large we can assure the positivity of W2(S,I,R).
Finally consider the positive function:(73)WS,I,R=W1S,I,R-W1S^,I^,R^+W2.Combining (63), (65), (67), (68), (69), and (70) we get (74)LWS,I,R≤-ΛS+βI+μ1+σ122+μ3-αIR+σβI+σ322+ω¯3Λ4ϵ-ω¯μ1-σ122-ϵΛS2-ω¯μ2-σ222-ϵΛI2-ω¯μ3-σ322-ϵΛR2+βΛμ1log⁡I-βSlog⁡I+βμ1βe-μ2-α-σ222S+σβ22μ1R2+βμ1σβ2+βS2-σβRlog⁡I+σβαμ3I2+σ2β2eμ3-σβσ222μ3R+σβ22μ3S2+σβμ3β2+σβR2-σ222+μ2+αlog⁡I+βSlog⁡I+σβRlog⁡I+σ222.Simplification and reorganization of the above inequality give (75)LWS,I,R≤-ΛS-ω¯μ1-σ122-ϵΛ-βμ1σβ2+β-σβ22μ3S2+βμ1βe-μ2-α-σ222S+βΛμ1-σ222-μ2-αlog⁡I+β+σβI-ω¯μ2-σ222-ϵΛ-σβαμ3I2+σ2β2eμ3-σβσ222μ3R-ω¯μ3-σ322-ϵΛ-σβμ3β2+σβR2-αIR+ω¯3Λ4ϵ+μ1+μ3+σ122+σ222+σ322≜Υ1S+Υ2I+Υ3R-αIR+C≜ΥS,I,R,where (76)Υ1S=-ΛSt-ω¯μ1-σ122-ϵΛ-βμ1σβ2+β-σβ22μ3S2+βμ1βe-μ2-α-σ222S,Υ2I=βΛμ1-σ222-μ2-αlog⁡I+β+σβI-ω¯μ2-σ222-ϵΛ-σβαμ3I2,Υ3R=σ2β2eμ3-σβσ222μ3R-ω¯μ3-σ322-ϵΛ-σβμ3β2+σβR2andC=ω¯3Λ4ϵ+μ1+μ3+σ122+σ222+σ322Since (61), then we can choose ϵ>0 such that (77)μi-σi22-ϵΛ>0,i=1,2,3,so, for sufficiently large ω¯, we have (78)ω¯μ1-σ122-ϵΛ-βμ1σβ2+β-σβ22μ3>0,ω¯μ2-σ222-ϵΛ-σβαμ3>0
and(79)ω¯μ3-σ322-ϵΛ-σβμ3β2+σβ>0.From (75) we have Υ(S,I,R)≤Υ1(S)+Υ2(I)+Υ3(R)+C and by (79) we get (80)limS→∞⁡ΥS,I,R=limI→∞⁡ΥS,I,R=limR→∞⁡ΥS,I,R=-∞,and from condition (60) which is equivalent to βΛ/μ1-σ22/2-μ2-α>0 we have (81)limS→0⁡ΥS,I,R=limI→0⁡ΥS,I,R=-∞.In the following, we will prove (62).
Firstly, sincelimS→∞Υ1(S)=limS→0Υ1(S)=-∞ and continuity of Υ1(S), Υ1(S) has a maximum point in R+3. So, let (82)M1=supS∈R+3⁡Υ1S.Similarly let (83)M2=supI∈R+3⁡Υ2IandM3=supR∈R+3⁡Υ3R
Case 1 (1/ρ≤I≤ρ). From (75), we have Υ(S,I,R)≤M1+M2+M3+C-1/ρR, hence (84)limR→0⁡ΥS,I,R=-∞
Case 2 (I>ρ). From (75), we have Υ(S,I,R)≤M1+M2+C+Υ3(R)-ρ/R, hence (85)limR→0⁡ΥS,I,R=-∞
Case 3 (I<1/ρ). From (75), we have Υ(S,I,R)≤M1+M3+C+Υ2(I) and since limI→0Υ2(I)=-∞, then we can choose ρ sufficiently large such that Υ2(I)<-(M1+M3+C)-2. Hence for a sufficiently large ρ we have (86)ΥS,I,R≤-2.Thereby and from (80) and (81) we deduce that, for a sufficiently large ρ, (87)LWS,I,R≤-1∀S,I,R∈R+3∖Δρ,So the condition (ii) of Theorem 1 holds.
The diffusion matrix associated with the system (3) is given by(88)ΣS,I,R=σ12S2000σ22I2000σ32R2It is easy to see that Σ(S,I,R) satisfy the uniform ellipticity condition (i) of Theorem 1.


## 5. Conclusion and Numerical Simulation

In the current paper, a stochastic SIR epidemic model with nonlinear relapse rate is considered to model the influence of infected individuals on the recovered ones. The thresholdRT which determines the dynamical behavior of the system (3) is found. Precisely, if RT<1 the disease will die out from the population with the probability one, while RT>1 leads to the persistence of the disease with a unique positive stationary distribution. We remark that RT is independent of the relapse coefficient σ. Hence, relapse under the pressure of the infected individuals has no influence on the dynamics of the stochastic model (3). However, one can remark from numerical simulation that the nonlinear relapse phenomenon plays a crucial role in the speed of the extinction and the spread of the disease in the population. Indeed, in Figure 1, we show that disease dies out quickly as long as σ is small enough. Furthermore, in Figure 2, we observe that, for higher coefficient σ, the modes of the invariant stationary distributions of infected individuals become larger.Figure 1
The numerical simulation of one path of the solution(St,It,Rt) of the system (3) using the Milstein scheme with (σ=0.035 left), (σ=0.35 right) and the parameters values Λ=1, β=0.01, μ1=0.066, μ2=0.066, μ3=0.0666, α=0.16, σ1=0.2, σ2=0.25, and σ3=0.21. In this situation RT=0.59<1. Hence, the disease will die out.Figure 2
Computer simulation of the density function of the invariant stationary distribution of system (3) with different values of σ and the parameters values Λ=2, β=0.02, μ1=0.066, μ2=0.066, μ3=0.0666, α=0.16, σ1=0.021, σ2=0.02, and σ3=0.023. In this situation RT=2,66>1.



---

*Source: 5493270-2018-06-25.xml*