The threshold dynamics in a stochastic SIS epidemic model with vaccination and nonlinear incidence under regime switching

Physica A 529 (2019) 121555

Contents lists available at ScienceDirect

Physica A journal homepage: www.elsevier.com/locate/physa

The threshold dynamics in a stochastic SIS epidemic model with vaccination and nonlinear incidence under regime switching ∗

Junna Hu, Zhidong Teng , Zhiming Li, Buyu Wen College of Mathematics and Systems Science, Xinjiang University, Urumqi 830046, People’s Republic of China

highlights • Stochastic SIS epidemic model with vaccination and nonlinear incidence under regime switching is investigated. • New technique to deal with nonlinear incidence and vaccination for stochastic epidemic models is introduced. • Sufficient conditions on the extinction and permanence in the mean are established.

article

info

Article history: Received 18 February 2017 Received in revised form 28 February 2018 Available online 24 May 2019 Keywords: Stochastic SIS epidemic model Markov switching Extinction Permanence Lyapunov function

a b s t r a c t In this paper, a stochastic SIS epidemic model with vaccination and nonlinear incidence under regime switching is investigated. This model concerned about white noise and color noise. A new technique to deal with the nonlinear incidence and vaccination for the stochastic epidemic models is introduced. The sufficient conditions on the extinction and the permanence in the mean with probability one are established by constructing new Lyapunov functions with regime switching. © 2019 Elsevier B.V. All rights reserved.

1. Introduction It is well known, vaccination to the susceptible is an important strategy to prevent and control the spread of infectious disease. In [1–12], we see that the deterministic epidemic models with vaccination are investigated, and the basic reproduction number, persistence and extinction of the disease, stability of equilibria, and the bifurcation phenomena, etc. for these models are extensively discussed. On the other hand, the epidemic systems are inevitably influenced by the environmental noise, which made the parameters of epidemic model are not fixed constants and they fluctuate with some average as the surrounding environment fluctuation, such as white noise and color noise (See [13–20] and the references cited therein). The white noise means the sum of many small, independent random disturbances, which from a series of almost continuous small perturbation characteristic parameters such as infection rate (See [21–24]). And the color noise can cause the population system to switch from one environmental regime to another, such as floods, earthquakes and so on, this type of environmental disturbances, different from the white noise fluctuation caused by rainfall or predation, can caused a large impact on the growth of organisms (See [25,26]). Therefore, the study of stochastic epidemic models has great ∗ Corresponding author. E-mail address: [email protected] (Z. Teng). https://doi.org/10.1016/j.physa.2019.121555 0378-4371/© 2019 Elsevier B.V. All rights reserved.

2

J. Hu, Z. Teng, Z. Li et al. / Physica A 529 (2019) 121555

significance to research the spread of infectious disease in reality life. In [25], the authors studied the following stochastic SIS epidemic model with vaccination under regime switching

[ ⎧ dS(t) = (1 − qr(t) )Ar(t) − βr(t) S(t)I(t) − (µr(t) + pr(t) )S(t) ⎪ ⎪ ] ⎪ ⎨ + γr(t) I(t) + εr(t) V (t) dt + σ1r(t) S(t)dB1 (t), [ ] ⎪ dI(t) = βr(t) S(t)I(t) − (µr(t) + γr(t) + αr(t) )I(t) dt + σ2r(t) I(t)dB2 (t), ⎪ ⎪ [ ] ⎩ dV (t) = qr(t) Ar(t) + pr(t) S(t) − (µr(t) + εr(t) )V (t) dt + σ3r(t) V (t)dB3 (t).

(1.1)

The sufficient conditions on the existence of stationary distribution are established. However, we see that the stochastic extinction and persistence of disease for this model have not been considered. As is well known, the nonlinear incidence rate is an important substance in modeling the dynamics of epidemic systems. In recent years, many authors have investigated the deterministic and stochastic epidemic models with nonlinear incidence rate (See [27–34] and the references cited therein). However, we see also that stochastic epidemic model with nonlinear incidence rate and vaccination under regime switching still has not been investigated. Therefore, in this paper we investigate the following stochastic SIS epidemic model with vaccination and nonlinear incidence under regime switching

[ ⎧ dS(t) = (1 − qr(t) )Ar(t) − βr(t) f (S(t))g(I(t)) − (µr(t) + pr(t) )S(t) ⎪ ⎪ ] ⎪ ⎨ + γr(t) I(t) + εr(t) V (t) dt + σ1r(t) S(t)dB1 (t), [ ] ⎪ dI(t) = βr(t) f (S(t))g(I(t)) − (µr(t) + γr(t) + αr(t) )I(t) dt + σ2r(t) I(t)dB2 (t), ⎪ ⎪ [ ] ⎩ dV (t) = qr(t) Ar(t) + pr(t) S(t) − (µr(t) + εr(t) )V (t) dt + σ3r(t) V (t)dB3 (t).

(1.2)

Our main purpose in this paper is to study the threshold dynamics of global positive solutions. Particularly, we will establish almost sufficient and necessary conditions on the stochastic extinction and permanence in the mean of the disease for model (1.2). Here, we will propose a new threshold value which is different from the corresponding threshold value given in [25]. Furthermore, we also will introduce a new techniques to deal with the nonlinear incidence functions and vaccination for the stochastic epidemic models. It will be presented in Section 2 (See Lemmas 2, 3 and 7 and the proofs of main theorems in Section 3. The rest of this paper is organized as follows. In Section 2, we will introduce some basic assumptions for model (1.2) and introduce some basic lemmas which will be used in the proofs of main results. In Section 3, we will state and prove the main results on the stochastic extinction and permanence for model (1.2). In Section 4, the numerical examples are given to illustrate the main results obtained in this paper. Lastly, in Section 5 a brief conclusion is given. 2. Preliminaries Denote R+ = (0, ∞), R+0 = [0, ∞) and Rn+ = {(x1 , x2 , . . . , xn ) : xi > 0, i = 1, 2, . . . , n}. For a bounded function f (t) defined for t ∈ R+0 , we denote f u = supt ∈R+0 f (t) and f l = inft ∈R+0 f (t). Let N be a positive integer, and set M = {1, 2, . . . , N }. For a sequence g = {g1 , g2 · · · , gN }, we denote gˆ = mink∈M gk and g˘ = maxk∈M gk . In model (1.2), S(t), I(t) and V (t) denote the numbers of susceptible, infectious and immune, respectively; r(t) for t ≥ 0 be a right-continuous Markov chain with values in a finite space M ; Bi (t) (i = 1, 2, 3) are independent standard Brownian 2 represents the intensity of Bi (t). The parameters Ar(t) is the input of new members into the susceptible, motion and σir(t) qr(t) is a fraction of vaccinated for new members, βr(t) is the disease transmission coefficient, µr(t) is the natural death rate of total population, γr(t) is the recovery rate of infectious, pr(t) denotes the proportional coefficient of vaccinated for the susceptible, εr(t) is the rate of losing their immunity for vaccinated individuals, and αr(t) represents the disease-caused death rate of infectious. All these parameters are assumed to be nonnegative, and µr(t) and Ar(t) also are positive. Throughout this paper, we assume that model (1.2) is defined in a complete probability space (Ω , F , P) with a filtration {Ft }t ≥0 satisfying the usual conditions. The generator Γ = (γij )N ×N of Markov chain r(t) is given by

{ γij △ + o(△), if i ̸= j, P {r(t + △) = j|r(t) = i} = 1 + γii △ + o(△), if i = j, ∑N

where △ > 0 and γij ≥ 0 is the transition rate from i to j if i ̸ = j while j=1 γij = 0. We assume that Bi (t) (i = 1, 2, 3) and r(t) are independent, and Markov chain r(t) is irreducible, which means that the system can switch from one environmental regime to another. That is to say, the Markov chain ∑ r(t) has a unique stationary distribution π = {π1 , π2 , . . . , πN }. It can be determined by equation π Γ = 0 subject to Ni=1 πi = 1 and πi > 0 for all i ∈ M . For functions f (S) and g(I) we introduce the following assumptions. (H1 ) Functions f (S) and g(I) are continuously differentiable for S ≥ 0 and I ≥ 0, respectively. f (0) = g(0) = 0, and g ′ (0) > 0. g(I) g ′ (0) (H2 ) I is nonincreasing for I > 0, and maxI >0 { g(I) − 1I } < ∞. ′ (H3 ) f (S) is nondecreasing for S ≥ 0. (H4 ) f ′ (S) is nonincreasing for S ≥ 0.

J. Hu, Z. Teng, Z. Li et al. / Physica A 529 (2019) 121555

Remark 1. When f (S) = When f (S) = (H2 ) holds.

S 1+ω1 S

Sm 1+ω1 S

3

or f (S) = S with constants m ≥ 2 and ω1 ≥ 0, then (H1 ) and (H3 ) clearly are satisfied.

with constant ω1 ≥ 0, then (H1 ) and (H4 ) are satisfied. When g(I) =

I 1+ω2 I

with constant ω2 ≥ 0, then

Lemma 1. Assume that (H1 ) holds. For any initial value (S(0), I(0), V (0)) ∈ R3+ , model (1.2) has a unique solution (S(t), I(t), V (t)) defined on t ∈ R+0 satisfying (S(t), I(t), V (t)) ∈ R3+ for all t ≥ 0 with probability one. Lemma 1 can be easily proved by using the same method which is given in Theorem 2.1 by Zhao et al. in [14]. We hence omit it here. Lemma 2.

The following equation

⎧ ⎪ ⎪ ⎨ dS0 (t) = (1 − qr(t) )Ar(t) − (µr(t) + pr(t) )S0 (t) + εr(t) V0 (t), dt

(2.1)

⎪ ⎪ ⎩ dV0 (t) = q A + p S (t) − (µ + ε )V (t), r(t) r(t) r(t) 0 r(t) r(t) 0

dt has a unique positive stationary distribution solution (S0 (t), V0 (t)). Proof. It is clear that the existence of unique positive stationary distribution solution (S0 (t), V0 (t)) of Eq. (2.1) is equivalent to the following equation

⎧ ∑ ⎪ γkl S0 (l) = 0, ⎪ ⎨(µk + pk )S0 (k) − εk V0 (k) − (1 − qk )Ak + ∑ l∈M ⎪ γkl V0 (l) = 0, ⎪ ⎩(µk + εk )V0 (k) − pk S0 (k) − qk Ak +

(2.2)

l∈M

for k = 1, 2, . . . , N has a unique positive solution (S0 (k), V0 (k), k = 1, 2, . . . , N). Eq. (2.2) can be rewritten as AC = β,

(2.3)

where C = (S0 (1), . . . , S0 (N), V0 (1), . . . , V0 (N)), β = ((1 − q1 )A1 , . . . , (1 − qN )AN , q1 A1 , . . . , qN AN ) , and T

µ1 + p1 + γ11 ··· ⎢ ⎢ γN1 ⎢ A=⎢ −p1 ⎣ ···

⎡

0

··· ··· ··· ··· ··· ···

γ1N ···

−ε1 ···

µN + pN + γNN

0

µ1 + ε1 + γ11 ··· γN1

0

··· − pN

⎤

0

··· ··· ··· ··· ··· ···

··· ⎥ ⎥ −εN ⎥ ⎥ γ1N ⎦ ··· µN + εN + γNN

For k = 1, 2, . . . , N, the leading principal sub-matrixes of A are

⎡ µ1 + p1 + γ11 γ21 ⎢ Ak = ⎣ ··· γk1 ⎡ AN +k

µ1 + p1 + γ11 ··· ⎢ ⎢ γN1 =⎢ ⎢ −p1 ⎣ ··· 0

γ12

µ2 + p2 + γ22 ··· γk2 ··· ··· ··· ··· ··· ···

··· ··· ··· ···

γ1k γ2k ··· µk + pk + γkk

γ1N ···

−ε1 ···

µN + pN + γNN

0

0

··· −pk

⎤ ⎥ ⎦ ··· ··· ··· ··· ··· ···

µ1 + ε1 + γ11 ··· γk1

We see that each column of sub-matrix Ak has the sum µi + pi +

∑k

∑N

j=1

0

⎤

··· ⎥ ⎥ −εk ⎥ ⎥ γ1k ⎦ ··· µk + εk + γkk

γji ≥ µi > 0 for i = 1, 2, . . . , k, and for sub-matrix ∑k

AN +k , the sum of ith column is µi + pi + j=1 γji − pi = µi > 0 (1 ≤ i ≤ N) and µi + εi + j=1 γji − εi ≥ µi > 0 (N < i ≤ N + k). Hence, Lemma 5.3 in [35] (or see [36]) implies detAk > 0 for k = 1, 2, . . . , 2N. By Theorem 2.10 in [35] (or see [36]), A is a nonsingular M-matrix, Therefore, for vector β ∈ R2N , which β is nonnegative and β ̸ = 0, Eq. (2.3) has a unique positive solution (S0 (k), V0 (k), k = 1, 2, . . . , N). This completes the proof. □ Lemma 3.

The following equation

⎧ ⎪ ⎪ ⎨ dc1 (t) = (µr(t) + pr(t) )c1 (t) − pr(t) c2 (t) − βr(t) f ′ (S0 (t))g ′ (0), dt

⎪ ⎪ ⎩ dc2 (t) = (µ + ε )c (t) − ε c (t). r(t) r(t) 2 r(t) 1 dt

(2.4)

4

J. Hu, Z. Teng, Z. Li et al. / Physica A 529 (2019) 121555

has a unique positive stationary distribution solution (c1∗ (t), c2∗ (t)). The proof of the lemma is similar to Lemma 2, so we here omit it. 2 2 2 ) > 0 for all k ∈ M . Let (S(t), I(t), V (t)) be the ∨ σ3k ∨ σ2k Lemma 4. Assume that (H1 ) and (H2 ) hold, and µk − 12 (σ1k 3 solution of model (1.2) with initial value (S(0), I(0), V (0)) ∈ R+ , then

lim

1(

S(t) + I(t) + V (t) = 0 a.s.

)

t →∞ t This lemma can been easily proved by using the similar method given in [37]. Furthermore, we also can obtain the following conclusion: 1 2 2 2 )} > 0. Then, ∨ σ3k ∨ σ2k ∨ 0)(σ1k Let ω(t) = S(t) + I(t) + V (t). Choosing a constant θ > 0 such that mink∈M {µk − ( θ − 2 there is a constant M > 0 such that the expectation

E (1 + ω(t))θ ≤ M

]

[

t ≥ 0.

for all

(2.5)

2 2 2 Lemma 5. Assume that (H1 ) and (H2 ) hold, and µk − 12 (σ1k ∨ σ2k ∨ σ3k ) > 0 for all k ∈ M . Let (S(t), I(t), V (t)) be the 3 solution of model (1.2) with initial value (S(0), I(0), V (0)) ∈ R+ , then for any sequence {g1 , g2 , . . . , gN } we have

lim

t →∞

1

t

∫

t

gr(s) S(s)dB1 (s) = 0, lim

t →∞

0

1

lim

t

t →∞

t

∫

t

∫

1 t

gr(s) I(s)dB2 (s) = 0, 0

gr(s) V (s)dB3 (s) = 0 a.s. 0

This lemma can been easily proved by using the similar method given in [37]. We hence omit it. 2 2 2 Lemma 6. Assume that (H1 ) and (H2 ) hold, and µk − 21 (σ1k ∨ σ2k ∨ σ3k ) > 0 for all k ∈ M . Let (S(t), I(t), V (t)) be the 3 solution of model (1.2) with initial value (S(0), I(0), V (0)) ∈ R+ , then

lim

t →∞

1

t

∫

t

0

I(s) g(I(s))

σ2γ (s) dB2 (s) = 0 a.s.

∫t

I(s) 0 g(I(s))

Proof. Let Z (t) =

σ2γ (s) dB2 (s). Choose a constant θ > 0 satisfying 2 < θ < 1 + mink∈M { σ 2 ∨σ2µ2k∨σ 2 }. Using the 1k

2k

3k

Burkholder–Davis–Gundy inequality (See [38]), we obtain that there is a constant Cθ > 0 such that for all t ≥ 0 E sup |Z (s)|θ ≤ Cθ E

[

]

t

∫ [

0≤s≤t

( 0

Let M0 = maxI ≥0 {

g ′ (0) g(I)

g(I(s))

]

0≤s≤t

( 1 )θ g ′ (0)

)2 σ22γ (s) ds

I(t) g(I(t))

− 1I }. We have

E sup |Z (s)|θ ≤ Cθ

[

I(s)

σ˘ 2θ E

≤

] θ2

.

(2.6)

1 (M0 I(t) g ′ (0)

+ 1) for all t ≥ 0. From (2.5) and (2.6) we can obtain

t

∫ [

(M0 I(s) + 1)2 ds

] θ2

θ

≤ C˜ θ t 2 ,

(2.7)

0

where C˜ θ = Cθ ( g ′1(0) )θ σ˘ 2θ M0θ M. For any constant δ > 0 and integers m = 1, 2, . . ., from (2.7) we have θ

E |Z ((m + 1)δ )|θ ≤ C˜ θ ((m + 1)δ ) 2 .

[

]

(2.8)

By the Doob’s martingale inequality (See [38]), for any constant ε0 > 0 from (2.8) we further have P{

θ

sup

mδ≤t ≤(m+1)δ

|Z (t)|θ > (mδ )1+ε0 + 2 } ≤

≤

θ

C˜ θ ((m + 1)δ ) 2 1+ε0 + θ2

(mδ )

E [|Z ((m + 1)δ )|θ ] θ

θ

≤

(mδ )1+ε0 + 2

C˜ θ 2 2 (mδ )1+ε0

.

Therefore, by the Borel–Cantelli lemma (See [38]), we have that for almost all ω ∈ Ω there is an enough large integer K0 θ such that supmδ≤t ≤(m+1)δ |Z (t)|θ ≤ (mδ )1+ε0 + 2 for all k ≥ K0 . Thus, for almost all ω ∈ Ω , if m ≥ K0 and mδ ≤ t ≤ (m + 1)δ ln |Z (t)|θ ln t

≤

(1 + ε0 + θ2 ) ln(mδ ) ln(mδ )

= 1 + ε0 +

Therefore, lim sup t →∞

ln |Z (t)|θ ln t

≤

1+

θ

θ 2

=

1 2

+

1

θ

.

θ 2

.

J. Hu, Z. Teng, Z. Li et al. / Physica A 529 (2019) 121555

5

For any enough small constant 0 < ζ < 12 − θ1 , there exists a set Ωζ ⊂ Ω such that P(Ωζ ) ≥ 1 − ζ , and for any ω ∈ Ωζ there is a T˜ = T˜ (ω) > 0 such that for any t ≥ T˜ we have ln |Z (t)| ≤ ( 21 + θ1 + ζ ) ln t. This shows for any ω ∈ Ωζ lim sup t →∞

1

|Z (t)|

≤ lim

t

t →∞

1

t 2 + θ +ζ t

=0 |Z (t)|

which together with the lim inft →∞ t ≥ 0 implies limt →∞ Z (t) limt →∞ t = 0 a.s. This completes the proof. □

|Z (t)| t

= 0 for any ω ∈ Ωζ . Therefore, we finally have

Lemma 7. Assume that (H1 ) holds. Let (S(t), I(t), V (t)) be the solution of model (1.2) with initial value (S(0), I(0), V (0)) ∈ R3+ , S0 (t) is given in Lemma 2, and c1∗ (t), c2∗ (t) are given in Lemma 3. Then t

∫

βr(s) f ′ (S0 (s))g ′ (0)(S(s) − S0 (s))ds ∫ t ] [ = Φ (t) − βr(s) f ′ (S0 (s))g ′ (0) + pr(s) c2∗ (s) − pr(s) c1∗ (s) + αr(s) c1∗ (s) I(s)ds, 0

(2.9)

0

where

Φ (t) =

∗ −c1∗ (t)(S(t) + I(t) − S0 (t)) − c2∗ (t)(V ∫ (t) − V0 (t)) + c1 (0)(S(0) + I(0) t

−S0 (0)) + c2∗ (0)(V (0) − V0 (0)) + c1∗ (s)σ1r(s) S(s)dB1 (s) ∫ t 0 ∫ t c2∗ (s)σ3r(s) V (s)dB3 (s). c1∗ (s)σ2r(s) I(s)dB2 (s) + +

(2.10)

0

0

Proof. Using the Itoˆ ’s formula, by Lemma 3 we obtain d c1∗ (t)(S(t) + I(t)) + c2∗ (t)V (t)

[

]

=

(c1∗ (t))′ (S(t) + I(t))dt + c1∗ (t)d(S(t) + I(t)) + (c2∗ (t))′ V (t)dt + c2∗ (t)dV (t)

=

[ ] −pr(t) c2∗ (t) − βr(t) f ′ (S0 (t))g ′ (0) + (µr(t) + pr(t) )c1∗ (t) (S(t) + I(t))dt {[ +c1∗ (t) (1 − qr(t) )Ar(t) − (µr(t) + pr(t) )S(t) + ϵr(t) V (t) ] } −(µr(t) + αr(t) )I(t) dt + σ1r(t) S(t)dB1 (t) + σ2r(t) S(t)dB2 (t) [ ] {[ + −εr(t) c1∗ (t) + (µr(t) + εr(t) )c2∗ (t) V (t)dt + c2∗ (t) qr(t) Ar(t) ] } +pr(t) S(t) − (µr(t) + εr(t) ) dt + σ3r(t) V (t)dB3 (t) { ∗ c1 (t)(1 − qr(t) )Ar(t) + c2∗ (t)qr(t) Ar(t) − βr(t) f ′ (S0 (t))g ′ (0)S(t) [ ] } − βr(t) f ′ (S0 (t))g ′ (0) + pr(t) c2∗ (t) − pr(t) c1∗ (t) + αr(t) c1∗ (t) I(t) dt

=

+c1∗ (t)σ1r(t) S(t)dB1 (t) + c1∗ (t)σ2r(t) I(t)dB2 (t) + c2∗ (t)σ3r(t) V (t)dB3 (t) Similarly, by Lemmas 2 and 3 we also obtain d c1∗ (t)S0 (t) + c2∗ (t)V0 (t)

[

]

=

(c1∗ (t))′ S0 (t)dt + c1∗ (t)dS0 (t) + (c2∗ (t))′ V0 (t) + c2∗ (t)dV0 (t)

=

] βr(t) f ′ (S0 (t))g ′ (0) − pr(t) c2∗ (t) + (µr(t) + pr(t) )c1∗ (t) S0 (t)dt [ ] +c1∗ (t) (1 − qr(t) )Ar(t) − (µr(t) + pr(t) )S0 (t) + εr(t) V0 (t) dt

=

[

+[−εr(t) c1∗ (t) + (µr(t) + εr(t) )c2∗ (t)]V0 (t)dt [ ] +c2∗ (t) qr(t) Ar(t) + pr(t) S0 (t) − (µr(t) + εr(t) )V0 (t) dt [ ∗ ] c1 (t)(1 − qr(t) )Ar(t) + c2∗ (t)qr(t) Ar(t) − βr(t) f ′ (S0 (t))g ′ (0)S0 (t) dt

Therefore, we have d c1∗ (t)(S(t) + I(t) − S0 (t)) + c2∗ (t)(V (t) − V0 (t))

[

=

]

−βr(t) f ′ (S0 (t))g ′ (0)(S(t) − S0 (t))− βr(t) f ′ (S0 (t))g ′ (0) ] } +pr(t) c2∗ (t) − pr(t) c1∗ (t) + αr(t) c1 (t) I(t) dt − c1∗ (t)σ1r(t) S(t)dB1 (t)

{

[

−c1∗ (t)σ2r(t) I(t)dB2 (t) − c2∗ (t)σ3r(t) V (t)dB3 (t).

(2.11)

6

J. Hu, Z. Teng, Z. Li et al. / Physica A 529 (2019) 121555

For any t > 0, integrating (2.11) from 0 to t we obtain c1∗ (t)(S(t) + I(t) − S0 (t)) + c2∗ (t)(V (t) − V0 (t))

[

]

[ ] − c1∗ (0)(S(0) + I(0) − S0 (0)) + c2∗ (0)(V (0) − V0 (0)) ∫ t ∫ t [ − βr(s) f ′ (S0 (s))g ′ (0)(S(s) − S0 (s))ds − βr(t) f ′ (S0 (t))g ′ (0) 0 0 ∫ t ] c1∗ (t)σ1r(t) S(t)dB1 (t) +pr(t) c2∗ (t) − pr(t) c1∗ (t) + αr(t) c1∗ (t) I(t)dt − 0 ∫ t ∫ t c1∗ (t)σ2r(t) I(t)dB2 (t) − − c2∗ (t)σ3r(t) V (t)dB3 (t).

=

(2.12)

0

0

From (2.12) we can obtain conclusion (2.9) of the lemma. □

3. Main results We define the threshold for model (1.2)

( σ2 ) πk βk f (S0 (k))g ′ (0) − µk − γk − αk − 2k .

∑

Rs0 =

2

k∈M

Theorem 1. Assume that (H1 ), (H2 ) and (H4 ) hold and αk ≥ pk for all k ∈ M . Let (S(t), I(t), V (t)) be the solution of model (1.2) with initial value (S(0), I(0), V (0)) ∈ R3+ . If Rs0 < 0, then ln I(t)

lim sup

t

t →∞

≤ Rs0 < 0 a.s.

This shows that the disease I in model (1.2) is exponentially extinct with probability one when Rs0 < 0. Furthermore, we also have

∫ ] 1 t µr(s) (S(s) − S0 (s))ds + µr(s) (V (s) − V0 (s))ds = 0, t →∞ t t 0 ∫ 0 ∫ ] [1 t 1 t pr(s) (S(s) − S0 (s))ds − (µr(s) + εr(s) )(V (s) − V0 (s))ds = 0. lim lim

t →∞

[1

t

∫

t

t

0

0

Particularly, when µr(t) ≡ µ > 0, pr(t) ≡ p ≥ 0 and εr(t) ≡ ε ≥ 0 are constants, then it follows that lim

t →∞

1 t

t

∫

S(s)ds = 0

∑

πk S0 (k),

k∈M

lim

t →∞

1 t

t

∫

∑

V (s)ds = 0

πk V0 (k).

(3.1)

k∈M

Proof. Using the Itoˆ ’s formula to the seconded equation of model (1.2), we have d ln I =

≤

[ g(I) 1 2 ] βr(t) f (S) − (µr(t) + γr(t) + αr(t) ) − σ2r(t) dt + σ2r(t) dB2 (t) I

2

βr(t) g (0)(f (S0 (t)) + f (S) − f (S0 (t))) − (µr(t) + γr(t) + αr(t) ) 1 2 ] − σ2r(t) dt + σ2r(t) dB2 (t). ′

[

(3.2)

2

Using the mean value theorem and then by (H4 ), we have that there exist a ζ1 which is situated between S and S0 (t) such that f (S) − f (S0 (t)) = f ′ (ζ1 )(S − S0 (t)) ≤ f ′ (S0 (t))(S − S0 (t)),

(3.3)

Hence, from (3.2) we further obtain d ln I(t) ≤

βr(t) f ′ (S0 (t))g ′ (0)(S(t) − S0 (t)) + βr(t) g ′ (0)f (S0 (t)) 1 2 ] −µr(t) − γr(t) − αr(t) − σ2r(t) dt + σ2r(t) dB2 (t)

[

2

(3.4)

J. Hu, Z. Teng, Z. Li et al. / Physica A 529 (2019) 121555

7

For any t > 0, integrating (3.4) from 0 to t we have t

∫ ln I(t) ≤

βr(s) f ′ (S0 (s))g ′ (0)(S(s) − S0 (s))ds

ln I(0) + t

∫ +

∫0 t +

0

[ 1 2 ] ds βr(s) g ′ (0)f (S0 (s)) − µr(s) − γr(s) − αr(s) − σ2r(s) 2

σ2r(s) dB2 (s).

0

From (2.9) in Lemma 7, we further obtain

∫

t

[ 1 2 ] βr(s) g ′ (0)f (S0 (s)) − µr(s) − γr(s) − αr(s) − σ2r(s) ds 2 0∫ t [ − βr(s) f ′ (S0 (s))g ′ (0)S0 (s) + pr(s) c2∗ (s) − pr(s) c1∗ (s) 0 ∫ t ] +αr(s) c1∗ (s) I(s)ds + ln I(0) + Φ (t) + σ2r(s) dB2 (s) 0 ∫ t [ 1 2 ] βr(s) g ′ (0)f (S0 (s)) − µr(s) − γr(s) − αr(s) − σ2r(s) ds 2 0 ∫ t + ln I(0) + Φ (t) + σ2r(s) dB2 (s)

ln I(t) ≤

≤

(3.5)

0

∫t

According to the strong law of large numbers, we have limt →∞ 1t 0 σ2r(s) dB2 (s) = 0. From (2.10), by Lemma 5 we can obtain limt →∞ 1t Φ (t) = 0. From the ergodic property of the Markov chain (See, for example [39]), we also have 1

lim

t

t →∞

t

∫

∑ [ 1 2 ] βr(s) g ′ (0)f (S0 (s)) − µr(s) − γr(s) − αr(s) − σ2r(s) ds = πk R0k . 2

0

k∈M

Therefore, from (3.5) we finally have lim sup

ln I(t) t

t →∞

≤ Rs0 < 0 a.s.

Let N(t) = S(t) + I(t) + V (t) and N0 (t) = S0 (t) + V0 (t), then we have d(N(t) − N0 (t)) =

[ ] −µr(t) (N(t) − N0 (t)) − αr(t) I(t) dt + σ1r(t) S(t)dB1 (t) + σ2r(t) I(t)dB2 (t) + σ3r(t) V (t)dB3 (t).

Hence,

=

) 1( N(t) − N0 (t) − N(0) + N0 (0) t ∫ t ∫ 1 1 t − µr(s) (N(s) − N0 (s))ds − αr(s) I(s)ds t 0 t 0 ] 1[ + σ1r(s) S(s)dB1 (s) + σ2r(s) I(s)dB2 (s) + σ3r(s) V (s)dB3 (s) . t

Since limt →∞ I(t) = 0 a.s., by Lemmas 4 and 5, we obtain that lim

[1

t →∞

t

∫

t

µr(s) (S(s) − S0 (s))ds + 0

1

t

∫

t

] µr(s) (V (s) − V0 (s))ds = 0.

(3.6)

0

From the third equation of model (1.2), we have d(V (t) − V0 (t)) =

pr(t) (S(t) − S0 (t)) − (µr(t) + εr(t) )(V (t) − V0 (t)) dt + σ3r(t) V (t)dB3 (t).

[

]

Similarly to above, we easily obtain that lim

[1

t →∞

t

t

∫

pr(s) (S(s) − S0 (s))ds − 0

1 t

t

∫

(µr(s) + εr(s) )(V (s) − V0 (s))ds = 0.

]

0

(3.7)

8

J. Hu, Z. Teng, Z. Li et al. / Physica A 529 (2019) 121555

Fig. 1. Simulations of the solution (S(t), I(t), V (t)) to system (1.2) in Example 1. The graphs are the sample pathes of S(t), I(t), V (t) and their density functions.

It is clear that, when µr(t) ≡ µ > 0, pr(t) ≡ p ≥ 0 and εr(t) ≡ ε ≥ 0 are constants, from (3.6) and (3.7) we further have lim

t →∞

1

t

∫

S(s)ds = lim

t

t →∞

0

1 t

t

∫

S0 (s)ds, lim

t →∞

0

1

t

∫

t

V (s)ds = lim 0

t →∞

1 t

t

∫

V0 (s)ds. 0

In view of the ergodic property of the Markov chain, we can get that lim

t →∞

1 t

t

∫

S0 (s)ds = 0

∑ k∈M

πk S0 (k), lim

t →∞

1 t

t

∫

V0 (s)ds = 0

∑

πk V0 (k).

k∈M

Therefore, (3.1) holds. This completes the proof. □ Remark 2. We see in Theorem 1 the conclusion (3.1) is established only in special case of µr(t) ≡ µ > 0, pr(t) ≡ p ≥ 0 and εr(t) ≡ ε ≥ 0 are constants. Therefore, an interesting open problem is to prove (3.1) for the general case.

J. Hu, Z. Teng, Z. Li et al. / Physica A 529 (2019) 121555

9

Fig. 2. The first graph is the sample path of r(t). The second is the corresponding mean value functions of S(t), I(t), V (t) in Example 1.

Theorem 2. Assume that (H1 ) − (H3 ) hold. If Rs0 > 0, then there is a constant m > 0 such that for any (S(t), I(t), V (t)) of model (1.2) with initial value (S(0), I(0), V (0)) ∈ R3+ one has lim inf t →∞

t

∫

1 t

I(s)ds ≥ m a.s. 0

This shows that the disease I in model (1.2) is permanent in the mean with probability one. Proof. Define a Lyapunov function I

∫ V (I) =

g ′ (0) g(I)

1

dI .

By the Itoˆ ’s formula, we have

[

dV (I) =

βr(t) f (t)g ′ (0) − (µr(t) + γr(t) + αr(t) ) 1

2 g ′ (0)g ′ (I) − σ2r(t)

( I )2 ]

2

g(I)

dt +

g ′ (0) g(I)

g ′ (0) g(I)

I

I σ2 (k)dB2 (t)

[ 1 2 ] βr(t) f (S0 (t))g ′ (0) − (µr(t) + γr(t) + αr(t) ) − σ2r(t) dt

=

(3.8)

2

( I +βr(t) g ′ (0)[f (S) − f (S0 (t))]dt − (µr(t) + γr(t) + αr(t) ) g ′ (0)

g(I)

1

( I )2

2 g ′ (0)g ′ (I) − σ2r(t)

(

2

g(I)

− 1 dt + σ2r(t) g ′ (0)

I

)

g(I)

) − 1 dt

dB2 (t).

According to the mean value theorem and (H3 ), there is a ξ (t) which is situated between S and S0 (t) such that f (S) − f (S0 (t)) = f ′ (ξ (t))(S − S0 (t)) ≥ f ′ (S0 (t))(S − S0 (t)). I From (H3 ), we have g ′ (0) g(I) − 1 ≤ M0 I for all I > 0, where M0 = supI >0 {

for all I > 0. Hence, g

′

( I )2

g ′ (0)g ′ (I)

g(I)

I (I) g(I)

(3.9) g ′ (0) g(I)

≤ 1 for all I > 0. From this, we further have

− 1 ≤ g ′ (0)

I g(I)

− 1I }. By (H1 ), we obtain ( g(I) )′ = I

g ′ (I)I −g(I) I2

≤0

− 1 ≤ M0 I

for all I > 0. Thus, from (3.8) we can obtain dV (I) ≥

[

1 2 ] βr(t) f (S0 (t))g ′ (0) − (µr(t) + γr(t) + αr(t) ) − σ2r(t) dt 2

+βr(t) f ′ (S0 (t))g ′ (0)(S − S0 (t))dt − (µr(t) + γr(t) + αr(t) )M0 Idt 1

I

2

g(I)

2 − σ2r(t) M0 Idt + g ′ (0)

σ2r(t) dB2 (t).

(3.10)

10

J. Hu, Z. Teng, Z. Li et al. / Physica A 529 (2019) 121555

Fig. 3. Simulations of the solution (S(t), I(t), V (t)) to system (1.2) in Example 2. The graphs are the sample pathes of S(t), I(t), V (t) and their density functions.

For ant t > 0, integrating (3.10) from 0 to t we have t

∫

V (I(t)) ≥

βr(s) f (S0 (s))g ′ (0) − (µr(s) + γr(s) + αr(s) ) ∫ t ] ds + βr(s) f ′ (S0 (s))g ′ (0)(S(s) − S0 (s))ds

[

V (I(0)) +

0

1

2 − σ2r(s) 2 0 ∫ t ( 1 2 ) − µr(s) + γr(s) + αr(s) + σ2r(s) M0 I(s)ds 2 0 ∫ t I(s) +g ′ (0) σ2r(s) dB2 (s). 0

g(I(s))

(3.11)

J. Hu, Z. Teng, Z. Li et al. / Physica A 529 (2019) 121555

11

Fig. 4. The first graph is the sample path of r(t). The second is the corresponding mean value functions of S(t), I(t), V (t) in Example 2.

Substituting (2.9) into (3.11), we further have 1 t

1

V (I(t)) ≥

t

V (I(0)) + 1

1

t

∫

t

βr(s) f (S0 (s))g ′ (0) − (µr(s) + γr(s) + αr(s) )

[ 0

1

∫

t

βr(s) f ′ (S0 (s))g ′ (0) + pr(s) c2∗ (s) ∫ ] 1 t [ µr(s) + γr(s) + αr(s) +αr(s) c1∗ (s) I(s)ds − t 0 ∫ 1 2 ] 1 1 t I(s) ′ + σ2r(s) M0 I(s)ds + Φ (t) + g (0) σ2r(s) dB2 (s). 2 − σ2r(s)

]

2

ds −

t

0

2

Since

g ′ (0) g(i)

1 I

≤

[

t

t

0

g(I(s))

+ M0 , we obtain that V (I) ≤ ln I + M0 I for all I > 0. Therefore,

]1 [ 1 ˘ 0 )g ′ (0) + p˘ c˘2∗ + α˘ c˘1∗ + (µ β˘ f ′ (S ˘ + γ˘ + α˘ + σ˘ 22 )M0 2

≥

1

t

∫

t

(3.12)

[

t

t

∫

I(s)ds 0

1 2 ] βr(s) f (S0 (s))g ′ (0) − (µr(s) + γr(s) + αr(s) ) − σ2r(s) ds 2

0

1

1

(3.13)

1

+ V (I(0)) − (ln I(t) + M0 I(t)) + Φ (t) t t t ∫ t 1 I(s) +g ′ (0) σ2r(s) dB2 (s) t

0

g(I(s))

1 t 1 lim supt →∞ t

By Lemma 6 we have limt →∞ limt →∞ 1t I(t) lim inf t →∞

= 0 and 1

σ2r(s) dB2 (s) = 0. From Lemmas 4 and 5 we also have limt →∞ 1t Φ (t) = 0, ln I(t) ≤ 0. Therefore, from (3.13) we finally obtain I(s) 0 g(I(s))

t

∫

t

∫t

I(s)ds ≥ 0

Rs0

˘ 0 )g ′ (0) + p˘ c˘ ∗ + α˘ c˘ ∗ + (µ β˘ f ′ (S ˘ + γ˘ + α˘ + 12 σ˘ 22 )M0 2 1

This completes the proof. □

Particularly, for model (1.1), we have f (S) = S, g(I) = I and Rs0 =

∑ k∈M

( 1 2) . πk βk S0 (k) − µk − γk − αk − σ2k 2

a.s.,

12

J. Hu, Z. Teng, Z. Li et al. / Physica A 529 (2019) 121555

Fig. 5. Simulations of the solution (S(t), I(t), V (t)) to system (1.2) in Example 3. The graphs are the sample pathes of S(t), I(t), V (t) and their density functions.

As a consequence of Theorems 1 and 2 we have the following result. Corollary 1. For model (1.1), we have that (1) if Rs0 < 0 and αk ≥ pk for all k ∈ M , then the disease I is exponentially extinct with probability one; (2) if Rs0 > 0, then the disease I is permanent in the mean with probability one. Remark 3. In Theorem 1, we see that the function f (S) is required to satisfy (H4 ). But, in Theorem 2 f (S) is required to satisfy (H3 ). However, (H3 ) and (H4 ) are incompatibly each other. In the numerical examples given in next section, we will see that the conclusions obtained in Theorems 1 and 2 many be true even if (H3 ) and (H4 ) do not satisfy. Therefore, an interesting open problem is to prove Theorems 1 and 2 without (H3 ) and (H4 ).

J. Hu, Z. Teng, Z. Li et al. / Physica A 529 (2019) 121555

13

Fig. 6. The first graph is the sample path of r(t). The second is the corresponding mean value functions of S(t), I(t), V (t) in Example 3.

4. Numerical examples In this paper, we only consider the r(t) is a right-continuous Markov chain taking values in M = {1, 2}. By Milstein’s higher order method (See [40]), we derive the corresponding discretization equations of model (1.2)

[ ] ⎧ Sj+1 =Sj + (1 − qk )Ak − βk f (Sj )g(Ij ) − (µk + pk )Sj + γk Ij + εk Vj ∆t ⎪ ⎪ ⎪ ⎪ 2 1 ⎪ σ1k Sj ∆t (ζ1j2 − 1) ⎪ ⎪ , + σ1k Sj ∆t2 ζ1j + ⎪ ⎪ ⎪ 2 ⎪ ⎪ 1 ⎨ ] [ Ij+1 =Ij + βk f (Sj )g(Ij ) − (µk + γk + αk )Ij ∆t + σ2k Ij ∆t2 ζ2j ⎪ 2 ⎪ Ij ∆t (ζ2j2 − 1) σ2k ⎪ ⎪ ⎪ , + ⎪ ⎪ 2 ⎪ ⎪ ⎪ 2 2 ⎪ 1 ⎪ ⎩ V =V + [q A + p S − (µ + ε )V ]∆ + σ V ∆ 2 ζ + σ3k Vj ∆t (ζ3j − 1) , j+1 j k k k j k k j t 3k j t 3j

(4.1)

2

where ξij (i = 1, 2, 3, j = 1, 2, . . .) are N(0, 1)-distributed independent Gaussian random variables and time increment ∆t > 0. Example 1. Take f (S) =

S , 1+ω1 S

g(I) =

I , 1+ω2 I

ω1 = 1, ω2 = 1.5, A = [0.65 0.7], q = [0.85 0.9], β = [0.6 0.8],

µ = [0.1 0.2], α = [0.65 0.6], p = [0.35 0.3], γ = [0.3 0.4], ε = [0.2 0.4], σ1 = [0.1 0.3], σ2 = [0.5 0.8] and 5 7 σ3 = [0.2 0.4]. The unique stationary of r(t) is given by π = (π1 , π2 ) = ( 12 , 12 ). It is easily verified that (H1 ), (H2 ) and (H4 ) hold. By computing, we have Rs0 = −0.9315 < 0. Therefore, by Theorem 1 we obtain that the disease I is extinct with probability one. The numerical simulations given in Fig. 1 illustrate the validity of theoretical results. Fig. 2 shows the sample path of Markov chain r(t) and the corresponding mean value function. Example 2. Take f (S) =

S2 , 1+ω1 S

g(I) =

I , 1+ω2 I

A = [2.6 2.7], q = [0.2 0.15], β = [4.8 4.5], µ = [0.65 0.5], α = [0.45 0.5],

p = [0.5 0.7], γ = [0.2 0.35], ε = [0.25 0.35], σ1 = [0.25 0.35], σ2 = [0.2 0.35], σ3 = [0.5 0.45], ω1 = 1.5 and ω2 = 1.5. Given the unique stationary distribution of r(t) is π = (π1 , π2 ) = ( 57 , 27 ). It is easily verified that (H1 ), (H2 ) and (H3 ) hold. By computing, we have Rs0 = 10.6106 > 0. Therefore, by Theorem 2 we obtain that the disease I is permanent in the mean with probability one. The numerical simulations given in Fig. 3 illustrate the validity of theoretical results. Fig. 4 shows the sample path of Markov chain r(t) and the corresponding mean value function. Example 3. Take f (S) =

S3 , 1+ω1 S 2

g(I) =

I , 1+ω2 I 2

A = [2.6 2.7], q = [0.8 1], β = [0.15 0.2], µ = [0.2 0.18], α = [0.55 0.7],

p = [0.4 0.6], γ = [0.2 0.25], ε = [0.15 0.2], σ1 = [0.1 0.3], σ2 = [0.6 0.6], σ3 = [0.1 0.35], ω1 = 1 and ω2 = 1.5. The unique stationary distribution of r(t) is π = (π1 , π2 ) = ( 45 , 15 ). It is found that (H2 ) − (H4 ) are not satisfied. Therefore, Theorems 1 and 2 are not applicable. By computing, we have Rs0 = −0.5636 < 0. However, the numerical simulations given in Fig. 5 illustrate that the disease I is extinct with probability one. Fig. 6 shows the sample path of Markov chain r(t) and the corresponding mean value function.

14

J. Hu, Z. Teng, Z. Li et al. / Physica A 529 (2019) 121555

Fig. 7. Simulations of the solution (S(t), I(t), V (t)) to system (1.2) in Example 4. The graphs are the sample pathes of S(t), I(t), V (t) and their density functions.

Example 4. Take f (S) and g(I) as in Example 3, ω1 = 2, ω2 = 1, A = [2.6 2.7], q = [0.2 0.15], β = [4.8 4.6], µ = [0.25 0.2], α = [0.35 0.45], p = [0.55 0.6], γ = [0.15 0.17], ε = [0.35 0.48], σ1 = [0.2 0.3], σ2 = [0.2 0.45] and 10 1 σ3 = [0.6 0.4]. The unique stationary distribution of r(t) is π = (π1 , π2 ) = ( 11 , 11 ). By computing, we have Rs0 = 6.5734 > 0. The numerical simulations given in Fig. 7 illustrate that the disease I is permanent in the mean with probability one. Fig. 8 shows the sample path of Markov chain r(t) and the corresponding mean value function. From above Examples 3 and 4, we see that even if nonlinear incidence functions f (S) and g(I) do not satisfy assumptions (H2 ) − (H4 ), then when Rs0 < 0, model (2) still is extinct and when Rs0 > 0, model (2) still is permanent in the mean with probability one. Therefore, we can propose the following conjecture. Conjecture. Assume that only (H1 ) holds, then for model (2) we have (1) when Rs0 < 0 and αk ≥ pk for all k ∈ M , then the disease I is exponentially extinct with probability one; (2) when Rs0 > 0, then the disease I is permanent in the mean with probability one.

J. Hu, Z. Teng, Z. Li et al. / Physica A 529 (2019) 121555

15

Fig. 8. The first graph is the sample path of r(t). The second is the corresponding mean value functions of S(t), I(t), V (t) in Example 4.

5. Conclusion In this paper, we propose a stochastic SIS epidemic model with vaccination and nonlinear incidence under regime switching. It is assumed that the epidemic system is stochastically influenced by the independent standard Brownian motions Bi (t) (i = 1, 2, 3) and Markov chain r(t) with values in a finite space. The new criteria on the extinction and permanence in the mean of the disease with probability one are established by constructing the new stochastic Lyapunov functions and introducing the new techniques to deal with the nonlinear incidence and vaccination. We see that, in order to establish our main results, some assumptions for nonlinear incidence functions f (S) and g(I) are given, see (H1 ) − (H4 ). We also observe that for many incidence functions, such as bilinear incidence and saturated incidence, see Remark 2.1, these assumptions are satisfied. However, there are some nonlinear incidence functions which I S3 do not satisfy assumptions (H2 ) − (H4 ), for example f (S) = 1+ω 2 and g(I) = 1+ω I 2 . But, the numerical simulations 1S 2 given in Examples 3 and 4 show that the conclusions established in Theorems 1 and 2 are still right. Therefore, an interesting open problem is to establish the stochastic extinction and permanence of the disease for model (1.2) in more weak assumptions than (H1 ) − (H4 ). It will be considered in the future. Acknowledgments We are very grateful to the editor and the reviewers for their very helpful comments and careful reading of our manuscript. This research is supported by the Natural Science Foundation of Xinjiang Province of China (Grant Nos. 2016D03022). References [1] C.M. Kribs-Zaleta, J.X. Velasco-Hernandez, A simple vaccination model with multiple epidemic ststes, Math. Biosci. 164 (2000) 183–201. [2] J. Li, Z. Ma, Qualitative analyses of SIS epidemic model with vaccination and varying total population size, Math. Comput. Modell. 35 (2002) 1235–1243. [3] J. Arino, C.C. McCluskey, P. van den Driessche, Global results for an epidemic model with vaccination that exhibits backward bifurcation, SIAM J. Appl. Math. 64 (2003) 260–276. [4] I. Moneim, D. Greenhalgh, Threshold and stability results for an SIRS epidemic model with a general periodic vaccination strategy, J. Biol. Syst. 13 (2005) 131–150. [5] F. Chen, A susceptible-infected epidemic model with voluntary vaccinations, J. Math. Biol. 53 (2006) 253–272. [6] T. Zhang, Z. Teng, An SIRVS epidemic model with pulse vaccination strategy, J. Theoret. Biol. 250 (2008) 375–381. [7] A. Lahrouz, L. Omari, D. Kiouach, Complete global stability for an SIRS epidemic model with generalized non-linear incidence and vaccination, Appl. Math. Comput. 218 (2012) 6519–6525. [8] X.-B. Zhang, H.-F. Huo, H. Xiang, An SIRS epidemic model with pulse vaccination and non-monotonic incidence rate, Nonlinear Anal. Hybrid Syst. 8 (2013) 13–21. [9] X. Feng, Z. Teng, K. Wang, Backward bifurcation and global stability in an epidemic model with treatment and vaccination, Disc. Cont. Dyn. Syst. B 19 (2014) 999–1025. [10] J. Yang, M. Martcheva, L. Wang, Global threshold dynamics of an SIVS model with waning vaccine-induced immunity and nonlinear incidence, Math. Biosci. 268 (2015) 1–8. [11] M. Shen, Y. Xiao, Global stability of a multi-group SVEIR epidemiological model with the vaccination age and infection age, Acta Appl. Math. 144 (2016) 137–157. [12] M.A. Pires, N. Crokidakis, Dynamics of epidemic spreading with vaccination: Impact of social pressure and engagement, Physica A 467 (2017) 167–179. [13] A. Gray, D. Greenhalgh, L. Hu, X. Mao, J. Pan, A stochastic differential equation SIS epidemic model, SIAM J. Appl. Math. 71 (3) (2011) 876–902. [14] Y. Zhao, D. Jiang, D. O’Regan, The extinction and peresistence of the stochastic SIS epidemic model with vaccination, Physica A 392 (2013) 4916–4927. [15] X. Mao, G. Marion, E. Renshaw, Environmental noise suppresses explosion in population dynamics, Stochastic Process. Appl. 97 (1999) 45–67.

16

J. Hu, Z. Teng, Z. Li et al. / Physica A 529 (2019) 121555

[16] R. Farnoosh, M. Parsamanesh, Stochastic differential equation systems for an SIS epidemic model with vaccination and immigration, Comm. Statist. Theory Methods 46 (2017) 8723–8736. [17] M. Slatkin, The dynamics of a population in Markovian environment, Ecology 59 (1978) 249–256. [18] B. Cao, M. Shan, Q. Zhang, W. Wang, A stochastic SIS epidemic model with vaccination, Physica A 486 (2017) 127–143. [19] X. Zhang, D. Jiang, T. Hayat, B. Ahmad, Dynamical behavior of a stochastic SVIR epidemic model with vaccination, Physica A 483 (2017) 94–108. [20] A. Gray, D. Greenhalgh, X. Mao, J. Pan, The SIS epidemic model with Markovian switching, J. Math. Anal. Appl. 394 (2012) 496–516. [21] Y. Cai, Y. Kang, M. Banerjee, W. Wang, A stochastic SIRS epidemic model with infectious fore under intervention strategies, J. Differential Equations 259 (2015) 7463–7502. [22] D. Zhao, T. Zhang, S. Yuan, The threshold of a stochastic SIVS epidemic model with nonlinear saturated incidence, Physica A 443 (2016) 372–379. [23] A. Lahrouz, A. Settati, A. Akharif, Effects of stochastic perturbation on the SIS epidemic system, J. Math. Biol. 74 (2017) 469–498. [24] Q. Liu, D. Jiang, N. Shi, T. Hayat, A. Alsaedi, Nontrivial periodic solution of a stochastic non-autonomous SISV epidemic model, Physica A 462 (2016) 837–845. [25] X. Zhang, D. Jiang, A. Alsaedi, T. Hayat, Stationary distribution of stochastic SIS epidemic model with vaccination under regime switching, Appl. Math. Lett. 59 (2016) 87–93. [26] D. Greenhalgh, Y. Liang, X. Mao, Modelling the effect of telegraph noise in the SIRS epidemic model using Markovian switching, Physica A 462 (2016) 684–704. [27] Q. Tang, Z. Teng, H. Jiang, Global behaviors for a class of multi-group SIRS epidemic models with nonlinear incidence rate, Taiwan. J. Math. 19 (2015) 1509–1532. [28] X. Feng, Z. Teng, F. Zhang, Global dynamics of a general class of multi-group epidemic models with latency and relapse, Math. Biosci. Eng. 12 (2015) 99–115. [29] A. Korobeinikov, Global asymptotic properties of virus dynamics models with dose-dependent parasite reproduction and virulence and nonlinear incidence rate, Math. Medic. Biol. IMA 26 (2009) 225–239. [30] Z. Shuai, P. van den Driessche, Global stability of infectious disease models using Lyapunov functious, SIAM J. Appl. Math. 73 (2013) 1513–1532. [31] Z. Teng, L. Wang, Persistence and extinction for a class of stochastic SIS epidemic models with nonlinear incidence rate, Physica A 451 (2016) 507–518. [32] D. Li, J. Cui, M. Liu, The evolutionary dynamics of stochastic epidemic model with nonlinear incidence rate, Bull. Math. Biol. 77 (2015) 1705–1743. [33] D. Zhao, T. Zhang, S. Yuan, The threshold of a stochastic SIVS epidemic model with nonlinear saturated incidence, Physica A 443 (2016) 372–379. [34] R. Rifhat, L. Wang, Z. Teng, Dynamics for a class of stochastic SIS epidemic models with nonlinear incidence and periodic coefficients, Physica A 481 (2017) 176–190. [35] X. Mao, C. Yuan, Stochastic Differential Equations with Markovian Switching, Imperial College Press, Lodon, 2006. [36] T. Huang, C. Yang, Special Matrix Analysis and Applications, Kexue Publishing Company, Beijing, 2007, in Chinese. [37] Y. Zhao, D. Jiang, The threshold of a stochastic SIS epidemic model with vaccination, Appl. Math. Comput. 243 (2014) 718–727. [38] X. Mao, Stochastic Differential Equations and Applications, second ed., Horwood Chichester, UK, 2008. [39] W.J. Anderson, Continuous-Time Markov Chain, Springer-Verlag, Berlin-Heidelberg, 1991. [40] D.J. Higham, An algorithmic introduction to numerical simulation of stochastic differential equations, SIAM Rev. 43 (2001) 525–546.