# Extinction and Ergodic Property of Stochastic SIS Epidemic Model with Nonlinear Incidence Rate

Authors: Qixing Han; Daqing Jiang; Chengjun Yuan
Journal: Abstract and Applied Analysis
 (2013)

---

## Abstract

We investigate a stochastic SIS model with nonlinear incidence rate. We show that there exists a unique nonnegative solution to the system, and condition for the infectious individualsI(t) to be extinct is given. Moreover, we prove that the system has ergodic property. Finally, computer simulations are carried out to verify our results.


---

## Body

## 1. Introduction

More attention has been paid to the epidemics models in order to monitor and curb the spread of some human diseases. A classical model is proposed by Kermack and McKendrick in 1927 [1]. They divided the population into three classes denoted by S(t), I(t), and R(t), which expressed the number of susceptible individuals, infective individuals, and removed individuals at time t, respectively. The model is called susceptible-infected-removed (SIR) model, and SIR models were investigated by many researchers [2–4].However some diseases, such as some sexually transmitted and bacterial diseases, do not have permanent immunity. In [5], they introduced a SIS model to describe the spread of the disease, which takes the following form:
(1)dS(t)=[μN-βS(t)I(t)+γI(t)-μS(t)]dt,dI(t)=[βS(t)I(t)-(μ+γ)I(t)]dt
with initial values S0+I0=N. N is the total size of the population. Where S(t) and I(t) express the number of susceptible individuals and infective individuals at time t, respectively, μ is the per capita birth (and death) rate, γ is the rate at which infected individuals become cured, and β is the per capita contact rate. In model (1), they assumed that a rate of contacts by an infective individual with a susceptible individual is proportional to population size, and model (1) has been well studied [6].In fact, population dynamics is inevitably affected by environmental white noise, which is always present. Since the parameters in the deterministic models are constant, they have some limitations when we describe the epidemics systems. Some researchers have paid their attention to the stochastic epidemics model [7–9]. Especially, in [10], Gray et al. consider that the parameter β in (1) is perturbed with
(2)β⟶β+σB˙(t),
where B(t) is Brownian motions and σ represents the intensities of the white noise. Corresponding to the deterministic model system (1), the stochastic system takes the following form:
(3)dS(t)=[μN-βS(t)I(t)+γI(t)-μS(t)]dt-σS(t)I(t)dB(t),dI(t)=[βS(t)I(t)-(μ+γ)I(t)]dt+σS(t)I(t)dB(t).
Since S(t)+I(t)=N, then (3) is reduced to
(4)dI(t)=[β(N-I(t))I(t)-(μ+γ)I(t)]dt+σ(N-I(t))I(t)dB(t).
For model (4), they pointed out that(i)
if(βN/(μ+γ)-σ2N2/2(μ+γ))<1 and σ2≤β/N, the disease I(t) will die out with probability one;(ii)
if(βN/(μ+γ)-σ2N2/2(μ+γ))>1, then model (4) has a unique stationary distribution.The incidence rate in (4) is bilinear, and several authors pointed out that the disease transmission process may have a nonlinear incidence rate [11, 12]. In [13], Xiao and Ruan propose an incidence rate
(5)g(I)S=βIS1+αI2,
where α is the parameter that measures the psychological or inhibitory effect, and 1/(1+αI2) describes the psychological or inhibitory effect from the behavioral change of the susceptible individuals when the number of infective individuals is very large. It can be used to explain some phenomena; for example, the outbreak of severe acute respiratory syndrome (SARS) had such psychological effects on the general public [14]: for a very large number of infective individuals, the infection force may decrease as the number of infective individuals increases. Some control measures and policies, such as border screening, mask wearing, quarantine, isolation, and so forth, can decrease the infection rate although the number of infective individuals was getting relatively larger. Equation (1) with nonlinear incidence rate (5) and the disturbed parameter β (β→β+σB˙(t)) can be written as follows:
(6)dS(t)=[μN-βS(t)I(t)1+αI2(t)+γI(t)-μS(t)]dt-σS(t)I(t)1+αI2(t)dB(t),dI(t)=[βS(t)I(t)1+αI2(t)-(μ+γ)I(t)]dt+σS(t)I(t)1+αI2(t)dB(t).
The parameters appearing in (6) have the same meaning as those above. Given that S(t)+I(t)=N, it is sufficient to study the SDE for I(t),
(7)dI(t)=[β(N-I(t))I(t)1+αI2(t)-(μ+γ)I(t)]dt+σ(N-I(t))I(t)1+αI2(t)dB(t),
with initial value I(0)=I0∈(0,N). Notice that when α=0, system (7) becomes (4). In this paper, we will analyze the dynamical behaviors of (7).The organization of this paper is as follows. In the next section, we show that there exists a unique positive solution to (7). In Section 3, we carry out a qualitative analysis of the model (7) and extinction conditions for I(t) is derived. We prove that the system has ergodic property under some condition in Section 4. In Section 5, we present some numerical simulations to illustrate our mathematical findings. A brief conclusion is given in Section 6.Throughout this paper, let(Ω,ℱ,P) be a complete probability space with a filtration {ℱt}t≥0 satisfying the usual conditions (i.e., it is increasing and right continuous while ℱ0 contains all P-null sets) and let B(t) be a scalar Brownian motion defined on the probability space.


## 2. Existence and Uniqueness of the Global Positive Solution

In order for the model to make sense, we need to show the solution is global and nonnegative. However, theorem of existence and uniqueness (cf. Arnold [15] and Mao [16]) is not satisfied in (7). By using tools established by Mao et al. [17], we will show existence and uniqueness of the global positive solution of (7).Theorem 1.
For any given initial dataI(0)=I0∈(0,N), there exists a unique solution I(t)∈(0,N) for all t≥0 with probability 1.Proof.
It is obvious that the coefficients of the SDE (7) are locally continuous. For any given initial data I(0)=I0∈(0,N), there exists a unique maximal local solution I(t) on t∈[0,τe), where τe is the explosion time. In order to show that the solution is global, it is sufficient to show τe=∞ a.s. Let m0>0 be sufficiently large so that I0 lies within the interval [1/m0,N-1/m0]. For each integer m≥m0, we define the stopping time
(8)τm=inf{t∈[0,τe):I(t)∉(1m,N-1m)},
where, throughout this paper, we set inf∅=∞ (as usual ∅ denotes the empty set). It is clear that τm is increasing as m→∞. Let τ∞=limm→∞τm, whence τ∞≤τe. It is easy to show that τ∞=∞ a.s. implies τe=∞ a.s. and I(t)∈(0,N) a.s. for all t≥0. Therefore, to complete this proof, it is enough to show that τ∞=∞ a.s.
Define a functionV:(0,N)→R+ as follows:
(9)V(x)=1x+1N-x.
By Itô’s formula, we get
(10)dV(x)={σ2x2(N-x)2(1+αx2)2(1x3+1(N-x)3)x(-1x2+1(N-x)2)[β(N-x)1+αx2-μ-γ]+σ2x2(N-x)2(1+αx2)2(1x3+1(N-x)3)[β(N-x)1+αx2-μ-γ]}dt+{[-1x2+1(N-x)2]σx(N-x)(1+αx2)}dB(t):=LV(x)dt+{[-1x2+1(N-x)2]σx(N-x)(1+αx2)}dB(t),
where
(11)LV(x)=x(-1x2+1(N-x)2)[β(N-x)1+αx2-μ-γ]+σ2x2(N-x)2(1+αx2)2(1x3+1(N-x)3)≤μ+γx+βNN-x+σ2N2(1x+1N-x)≤[(μ+γ)∨(βN)+σ2N2]V(x).
By almost the same method in the proof of [10], the desired result will be obtained.


## 3. Extinction

In this section, we will point out the condition forI(t) to be extinct. We firstly do some preparation work.Consider the following stochastic equation:(12)dX(t)=b(X(t))dt+σ(X(t))dB(t),
and assume that the coefficients σ:J→R, b:J→R satisfy(1)
σ2(x)>0, ∀x∈J,(2)
∀x∈I, ∃ϵ>0 such that ∫x-ϵx+ϵ‍((1+|b(y)|)/σ2(y))dy<∞,whereJ=(l,r); -∞≤l<r≤∞.Lemma 2 (see [18]).
Assume that (1) and (2) hold, and letX(t) be a weak solution of (12) in J, with nonrandom initial condition X0=x∈J. Let p be given by
(13)p(x)=∫cx‍e-∫cv‍(2b(y)/σ2(y))dydv,c∈J.
If p(l+)>-∞, p(r-)=∞, then
(14)P(limt→∞X(t)=l)=P(supt≥0X(t)<r)=1.Theorem 3.
IfR0s:=2(βN-μ-γ)/σ2N2<1, then for any initial data I(0)=I0∈(0,N), the solution of SDE (7) has the following property:
(15)P(limt→∞I(t)=0)=1;
that is, the disease dies out with probability one.Proof.
Applying Lemma2 with b(x)=β(N-x)x/(1+αx2)-(μ+γ)x, σ(x)=σ(N-x)x/(1+αx2) and c∈J=(0,N), we can compute
(16)∫cx‍2b(τ)σ2(τ)dτ=2σ2∫cx‍β(N-τ)(1+ατ2)-(μ+γ)(1+ατ2)2τ(N-τ)2dτ=2σ2{βN-(μ+γ)N2lnx-[βN-(μ+γ)N2+α(N2βN+2(μ+γ)+3α(μ+γ)N2)βN-(μ+γ)N2]ln(N-x)-((μ+γ)/N)(N2α+1)2N-x-(μ+γ)α22x2-(2(μ+γ)α2N+αβ)x[βN-(μ+γ)N2}+C0.
Clearly, conditions (1) and (2) are satisfied. By calculation, we have
(17)p(x)=∫cx‍exp{-∫cs‍2b(τ)σ2(τ)dτ}ds=e-c0∫cx‍s-2(βN-(μ+γ))/σ2N2×(N-s)2[(βN-(μ+γ))/N2+α(βN+2(μ+γ)+3α(μ+γ)N2)]/σ2×e2(μ+γ)(N2α+1)2/Nσ2(N-s)×e((μ+γ)α2s2+(4(μ+γ)α2N+2αβ)s)/σ2ds.
We see that(18)p(N-)=e-c0∫cN‍s-2(βN-(μ+γ))/σ2N2×(N-s)2[(βN-(μ+γ))/N2+α(βN+2(μ+γ)+3α(μ+γ)N2)]/σ2×e2(μ+γ)(N2α+1)2/Nσ2(N-s)×e((μ+γ)α2s2+(4(μ+γ)α2N+2αβ)s)/σ2ds.
Let 1/(N-s)=t; we obtain that
(19)p(N-)=e-c0∫1/(N-c)∞‍(Nt-1)-2(βN-(μ+γ))/σ2N2t2(βN-(μ+γ))/σ2N2×t-2[(βN-(μ+γ))/N2+α(βN+2(μ+γ)+3α(μ+γ)N2)]/σ2×e2(μ+γ)(N2α+1)2t/Nσ2×e((μ+γ)α2(N-1/t)2+(4(μ+γ)α2N+2αβ)(N-1/t))/σ2t-2dt=e-c0∫1/(N-c)∞‍(Nt-1)-2(βN-(μ+γ))/σ2N2×t-(2α(βN+2(μ+γ)+3α(μ+γ)N2)/σ2)-2×e2(μ+γ)(N2α+1)2t/Nσ2×e((μ+γ)α2(N-1/t)2+(4(μ+γ)α2N+2αβ)(N-1/t))/σ2dt.
Clearly,
(20)p(N-)=∞.
When R0s<1, we have
(21)-p(0+)=e-c0∫0c‍s-2(βN-(μ+γ))/σ2N2×(N-s)2[(βN-(μ+γ))/N2+α(βN+2(μ+γ)+3α(μ+γ)N2)]/σ2×e2(μ+γ)(N2α+1)2/Nσ2(N-s)×e((μ+γ)α2s2+(4(μ+γ)α2N+2αβ)s)/σ2ds<∞;
that is,
(22)p(0+)>-∞.
It can be straightly shown by Lemma 2 that
(23)P(limt→∞I(t)=0)=1.
The proof is therefore completed.


## 4.  Ergodic Property

In this section, we show that the small perturbation forces the infective individuals to be ergodic.Theorem 4.
IfR0s:=2(βN-μ-γ)/σ2N2>1, then, for any initial data I(0)=I0∈(0,N), the solution of SDE (7) is ergodic.Proof.
IfR0s>1, we get
(24)∫0c‍exp{-∫cs‍2b(τ)σ2(τ)dτ}ds=e-c0∫0c‍s-2(βN-(μ+γ))/σ2N2×(N-s)2[(βN-(μ+γ))/N2+α(βN+2(μ+γ)+3α(μ+γ)N2)]/σ2×e2(μ+γ)(N2α+1)2/Nσ2(N-s)×e((μ+γ)α2s2+(4(μ+γ)α2N+2αβ)s)/σ2ds=∞.
It follows from (20) that
(25)∫cN‍exp{-∫cs‍2b(τ)σ2(τ)dτ}ds=p(N-)=∞.
Besides,
(26)m=∫0N‍1σ2(s)exp{∫cs‍2b(τ)σ2(τ)dτ}ds=ec0σ2∫0N‍(1+αs2)2×s(2(βN-(μ+γ))/σ2N2)-2×(N-s)-(2[(βN-(μ+γ))/N2+α(βN+2(μ+γ)+3αμN2)]/σ2)-2×e-2(μ+γ)(N2α+1)2/Nσ2(N-s)×e-((μ+γ)α2s2+(4(μ+γ)α2N+2αβ)s)/σ2ds.
Let t=N/(N-x)-1; then
(27)m=ec0σ2∫0∞‍[1+α(Ntt+1)2]2×(Ntt+1)(2(βN-(μ+γ))/σ2N2)-2×(t+1)(2[(βN-(μ+γ))/N2+α(βN+2(μ+γ)+3α(μ+γ)N2)]/σ2)+2×e-2(μ+γ)(N2α+1)2(t+1)/N2σ2×e-((μ+γ)α2(Nt/(t+1))2+(4(μ+γ)α2+2αβ)(Nt/(t+1)))/σ2×(t+1)-2dt=ec0σ2N(2(βN-(μ+γ))/σ2N2)-2×e-2(μ+γ)(N2α+1)2/N2σ2×∫0∞‍[αN2t2+(t+1)2]2t(2(βN-(μ+γ))/σ2N2)-2×(t+1)(2α(βN+2(μ+γ)+3α(μ+γ)N2)/σ2)-2×e-2(μ+γ)(N2α+1)2t/N2σ2×e-(((μ+γ)α2(Nt/(t+1))2+(4(μ+γ)α2+2αβ)(Nt/(t+1)))/σ2)dt.
When R0s>1, that is, 2-(2(βN-(μ+γ))/σ2N2)<1, we have
(28)m=∫0N‍1σ2(s)exp{∫cs‍2b(τ)σ2(τ)dτ}ds<∞.
The conditions of Theorem 1.16 in Kutoyants [19] follow from (24), (25), and (28). Therefore I(t) is ergodic, and the invariant density is given by
(29)π(x)=C(1+αx2)2×x(2(βN-(μ+γ))/σ2N2)-2×(N-x)-(2[(βN-(μ+γ))/N2+α(βN+2(μ+γ)+3α(μ+γ)N2)]/σ2)-2×e-2(μ+γ)(N2α+1)2/Nσ2(N-x)×e-((μ+γ)α2x2+(4(μ+γ)α2N+2αβ)x)/σ2,x∈(0,N),
where C is a constant.Remark 5.
Ifα=0, then
(30)π(x)=Cx(2[βN-(μ+γ)]/σ2N2)-2×(N-x)(-2[βN-(μ+γ)]/σ2N2)-2×e-2(μ+γ)/Nσ2(N-x),x∈(0,N).
By the property of density function, we have
(31)1=∫0N‍Cx(2[βN-(μ+γ)]/σ2N2)-2×(N-x)(-2[βN-(μ+γ)]/σ2N2)-2×e-2(μ+γ)/Nσ2(N-x)dx.
Let x=Nt/(t+1); then
(32)∫0NCx(2[βN-(μ+γ)]/(σ2N2))-2×(N-x)(-2[βN-(μ+γ)]/σ2N2)-2e-2(μ+γ)/Nσ2(N-x)dx=CN-3e-2(μ+γ)/σ2N2×∫0∞t2[βN-(μ+γ)]/σ2N2(t+1)2e-2(μ+γ)t/σ2N2dt=CN-3e-2(μ+γ)/σ2N2(N2σ22(μ+γ))2[βN-(μ+γ)]/σ2N2×[(βN-μ-γ)(2(βN-μ-γ)-σ2N2)(μ+γ)σ2N2+4[βN-(μ+γ)]-2σ2N2σ2N2+2(μ+γ)σ2N2(βN-μ-γ)(2(βN-μ-γ)-σ2N2)(μ+γ)σ2N2]×Γ(2(βN-μ-γ)σ2N2-1),
which implies
(33)1C=N-3e-2(μ+γ)/σ2N2×(N2σ22(μ+γ))2[βN-(μ+γ)]/σ2N2×[(βN-μ-γ)(2(βN-μ-γ)-σ2N2)(μ+γ)σ2N2+4(βN-μ-γ)-2σ2N2σ2N2+2(μ+γ)σ2N2(βN-μ-γ)(2(βN-μ-γ)-σ2N2)(μ+γ)σ2N2]×Γ(2(βN-μ-γ)σ2N2-1).
Next, we compute E(X):
(34)E(X)=∫0N‍Cx(2[βN-(μ+γ)]/σ2N2)-1×(N-x)(-2[βN-(μ+γ)]/σ2N2)-2e-2(μ+γ)/Nσ2(N-x)dx=CN-2e-2(μ+γ)/σ2N2(N2σ22(μ+γ))2[βN-(μ+γ)]/σ2N2×2(βN-μ-γ)-σ2N2σ2N2βN(μ+γ)×Γ(2(βN-μ-γ)σ2N2-1).
This together with (33) shows that
(35)E(X)=β[2(βN-μ-γ)-σ2N2]2β2-(μ+γ+βN)σ2.
Similarly, we have
(36)E(X2)=(βN-μ-γ)[2(βN-μ-γ)-σ2N2]2β2-(μ+γ+βN)σ2=(βN-μ-γ)E(X)β.
Consequently,
(37)Var(X)=(βN-μ-γ)E(X)β-[E(X)]2.
When α=0, model (7) becomes (4), and the mean and variance of the stationary distribution of model (4) are the same as the results in [10].


## 5. Simulations

We utilize the method developed in [20] to illustrate our findings. We have the corresponding discretization equation:
(38)Ik+1=Ik+Ik[β(N-Ik)1+αIk2-(μ+γ)Ik]Δt+σIkN-Ik1+αIk2εkΔt,
where εk are the Gaussian random variables N(0,1).Using the discretized equation and with the help of Matlab software, choosing the appropriate parametersN=2, β=0.5, μ=0.2, γ=0.3, and α=0.5, we get simulations of system (7) and the corresponding deterministic system.In Figure1, we choose initial values I(0)=1.8 and 0.4, respectively, σ=0.6, noting that 2(βN-μ-γ)/σ2N2=0.694<1, as the case in Theorem 3 expected, for any initial value I0∈(0,2), the large white noise leads to the extinction of I(t) and the solution of system (7) tends to zero; that is, the disease dies out. While the solution of corresponding deterministic system does not tend to zero.The solutionI(t) of SIS model (7) compared to the deterministic system with different initial value I(0)=1.8 and 0.4, σ=0.6. The blue line represents the solution of stochastic system, and the black line represents the solution of the corresponding deterministic system.
(a)
(b)In Figure2, we choose σ=0.2 and 0.05, respectively, with initial values I(0)=1.5, which satisfy the cases in Theorem 4; that is, 2(βN-μ-γ)/σ2N2>1. From the left picture in Figure 2, we can see that the solution of system (7) is fluctuating in a small neighborhood, and there is a stationary distribution (see the histogram in Figures 2(b) and 2(d)); the disease becomes endemic.The solutionI(t) of SIS model (7) and its histogram with differing value of σ=0.2,0.05.
(a)
(b)
(c)
(d)


## 6. Conclusion

In this paper, we analyze the dynamic behaviors of a stochastic SIS model with nonlinear incidence rate, under the assumption that the population lives in an environment subjected to random fluctuations which mainly affect the disease transmission term. First of all, we show that there exists a unique positive solution in system (7). Moreover, we obtain the threshold between prevalence and extinction of the disease; that is, if R0s=2(βN-μ-γ)/σ2N2<1, the disease will die out with probability one; if R0s=2(βN-μ-γ)/σ2N2>1, I(t) is ergodic, which means the disease will become endemic. Finally, we illustrate our results with computer simulations.We see that ifα=0, (7) becomes (4), which is studied in [10]; they show that if βN/(μ+γ)-σ2N2/2(μ+γ)<1 and σ2≤β/N, the disease I(t) will die out with probability one. Obviously, βN/(μ+γ)-σ2N2/2(μ+γ)<1 and R0s=2(βN-μ-γ)/σ2N2<1 are equivalent. The condition σ2≤β/N, however, is not necessary in our investigation. In addition, if R0s=2(βN-μ-γ)/σ2N2>1, the invariant density is obtained, which is not mentioned in [10].An extension of our work is to consider a stochastic SIS model with the general incidence rateβIpS/(1+αIq), and it is currently a work in progress.



---

*Source: 127321-2013-12-04.xml*