Physica A 462 (2016) 684–704

Contents lists available at ScienceDirect

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

Modelling the effect of telegraph noise in the SIRS epidemic

model using Markovian switching
D. Greenhalgh ∗ , Y. Liang, X. Mao
Department of Mathematics and Statistics, University of Strathclyde, Glasgow G1 1XH, UK

highlights
• We introduce telegraph noise with Markov Switching into the SIRS epidemic model.
• Establish extinction and persistence conditions.
• SIR model a special case.
• Analytical results confirmed by simulation.

article info abstract

Article history: We discuss the effect of introducing telegraph noise, which is an example of an
Received 30 June 2015 environmental noise, into the susceptible–infectious–recovered–susceptible (SIRS) model
Received in revised form 26 May 2016 by examining the model using a finite-state Markov Chain (MC). First we start with a
Available online 25 June 2016
two-state MC and show that there exists a unique nonnegative solution and establish
the conditions for extinction and persistence. We then explain how the results can be
Keywords:
generalised to a finite-state MC. The results for the SIR (Susceptible–Infectious–Removed)
Markov Chain
SIRS epidemic model
model with Markovian Switching (MS) are a special case. Numerical simulations are
Extinction produced to confirm our theoretical results.
Persistence © 2016 The Author(s). Published by Elsevier B.V. This is an open access article under the
Environmental noise CC BY license (http://creativecommons.org/licenses/by/4.0/).
Lyapunov stability

1. Introduction

The dynamics of population systems are often influenced by different types of environmental noise for example white
noise or Brownian motion. The effects of white noise have already been considered by various authors (e.g. Refs. [1–4]).
Environmental noise has the potential to have a huge impact on the population dynamics of a system. For example it has
been shown that sufficiently large white noise can cause a population that would otherwise explode or tend to a unique
endemic equilibrium to die out [2,5]. In this paper, we will focus on another type of environmental noise, namely telegraph
noise (or burst noise). This consists of sudden instantaneous transitions between two or more sets of parameter values in
the underlying model corresponding to two or more different environments or regimes (e.g. Refs. [6–8]). Switching between
environments follows a finite state continuous time Markov Chain (MC) with state space S = {1, 2, . . . , M }, where M is the
number of different environments. Hence the switching times are memoryless and follow an exponential distribution.

∗ Corresponding author.
E-mail address: david.greenhalgh@strath.ac.uk (D. Greenhalgh).

http://dx.doi.org/10.1016/j.physa.2016.06.125
0378-4371/© 2016 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/
by/4.0/).
 D. Greenhalgh et al. / Physica A 462 (2016) 684–704 685

There are already various papers which have looked at the effect of telegraph noise in a population system model. As an
example, let us consider a Lotka–Volterra predator–prey model,

ẋ(t ) = x(t )(a − by(t )),

(1.1)
ẏ(t ) = y(t )(−c + dx(t )),

where x is the number of prey and y is the number of predators at time t. a, b, c and d are positive constants.
So if r (t ) is a continuous time Markov Chain with state space {1, 2} this Lotka–Volterra model with Markov Switching
becomes

ẋ(t ) = x(t )(ar (t ) − br (t ) y(t )),

(1.2)
ẏ(t ) = y(t )(−cr (t ) + dr (t ) x(t )),

where ai , bi , ci and di are positive constants.

This Lotka–Volterra predator–prey model has been looked at in various papers. For example Takeuchi et al. [9] studied
the behaviour of this system. Note that under each environment the numbers of predators and prey follow a deterministic
predator–prey equation, switching between environments according to telegraph noise. Takeuchi et al. have shown that if
the two equilibrium states of the two subsystems differ, all positive trajectories of the system always exit from any compact
set in R2+ with probability one. On the other hand, if the two equilibrium states coincide, then the trajectory either exits from
any compact set in R2+ or converges to the common equilibrium point. These properties imply that the population system
(1.1) under telegraph noise is neither permanent nor dissipative [9]. Du et al. [6] investigated the impact that telegraph
noise has on the behaviour of Lotka–Volterra competition systems. The oscillatory behaviour of the solution to the systems
with telegraph noise was observed. Li et al. [10] looked at a more generalised Lotka–Volterra model with n interacting
species described by an n-dimensional system of ordinary differential equations. In their paper they looked at the effect that
two different types of environmental noise have on the system. First of all they introduced white noise into the model in
the form of Brownian motion. They then took a further step by considering telegraph noise using a finite-state MC. They
obtained the existence conditions for the system to have global positive solutions as well as the conditions for the solutions
to be ultimately bounded and permanent. Furthermore, they also established the extinction conditions.
In this paper, we want to look at the effect that telegraph noise has on the dynamics of an SIRS epidemic model. In the
1920s, Kermack and McKendrick [11] constructed the SIR and SIRS epidemic models to illustrate respectively diseases where
there is a permanent acquired immunity such as measles [12] and where there is a temporary acquired immunity such as
rubella. The SIR model is a special case of the SIRS model. We will explain later on in this paper that the results for our SIRS
model also apply to the SIR model. There has been much research done on different aspects of both SIR and SIRS epidemic
models. For example Hethcote [12] shows that the behaviour of each model is determined by a threshold parameter (the
basic reproduction number R0 ). If R0 is less than or equal to one or no disease is initially present then the system tends to
the unique disease-free equilibrium (DFE), but if R0 exceeds one then the system tends to a unique endemic equilibrium.
Tornatore et al. [13] propose a stochastic SIR model with environmental white noise added into the disease transmission
term with or without distributed time delay and study the stability of the DFE. The numerical simulation of the stochastic SIR
model shows that the introduction of noise modifies the threshold of the system for an epidemic to occur and the threshold
value is found. Lu [14] later extended their results into an SIRS model.
Yang et al. [15] introduce stochastic environmental noise into the death rates for SIR and SEIR (susceptible–exposed–
infectious–removed) epidemic models with different death rates for different population classes. They investigated the
dynamics of the models depending on the basic reproduction number R0 . The long-term behaviour of the two stochastic
systems is studied. The authors mainly use stochastic Lyapunov functions to show that under certain conditions, the
solutions are ergodic if R0 > 1, and that they are exponentially stable when R0 ≤ 1. Finally they show numerically that
the analytical results are true.
Zhao and Jiang [16] studied the dynamics of a stochastic SIRS epidemic model with saturated incidence. The disease
transmission term is β SI /(1 + α I ), where α is a constant, and there are deaths due to the disease. They introduce
environmental white noise into the disease transmission parameter β in a similar way to Gray et al. [17]. They obtain a
threshold value of the stochastic system which determines the extinction and persistence of the epidemic. They also show
that large noise will suppress the epidemic.
O’ Regan et al. [18] constructed a new Lyapunov function for a variety of deterministic SIR and SIRS models in
epidemiology. They considered the SIR and SIRS models with proportional disease incidence and deaths due to the disease.
They used this to establish the global stability of the endemic equilibrium states in these models. On the other hand
Korobeinikov [19] constructed different Lyapunov functions for two-dimensional SIR and SIRS compartmental epidemic
models with nonlinear transmission rate of a very general form f (S , I ) subject to a few biologically relevant conditions. The
models included some with vertical and horizontal transmission. Korobeinikov shows that provided that the population size
is constant and f (S , I ) is concave in I, the number of infectious individuals, then the positive endemic equilibrium state is
globally stable.
Vargas de Leon [20] establishes the global stability conditions for classic deterministic SIS, SIR and SIRS epidemic
models with constant recruitment, disease-induced death and standard incidence rate. He uses novel methods to construct
686 D. Greenhalgh et al. / Physica A 462 (2016) 684–704

Lyapunov functions and shows that for the SIRS model the unique endemic equilibrium is globally stable under certain
parameter conditions.
Liu and Stechlinski [21] consider pulse and constant control schemes for deterministic SIR epidemic models with
seasonality in the contact rate. A constant treatment scheme is applied to the model. Easily verifiable conditions on the
basic reproduction number of the infectious disease are established which ensure disease eradication under these constant
control strategies. Later both pulse vaccination and pulse treatment models are applied to an SIR model with time-varying
contact rate. Further, a vaccine failure model as well as a model with a reduced infective class are considered with pulse
control schemes. Again conditions on the basic reproduction number are developed which ensure disease eradication.
Nasell [22] considers stochastic models of some endemic infections with demography. Approximations of quasi-
stationary distributions and times to extinction are derived for stochastic versions of the SI (susceptible–infectious), SIS,
SIR and SIRS epidemic models. The approximations are valid for sufficiently large population sizes. Conditions for validity of
the approximations are given for each of the models. There are also conditions for validity of the corresponding deterministic
model. It is noted that some deterministic models are unacceptable approximations of the stochastic models for a large range
of realistic parameter values.
Chen and Li [23] introduced the effect of white noise into the SIR epidemic model and the time delayed SIR epidemic
model. This was done by adding a separate independent Brownian motion term to each of the per capita susceptible and
infectious death rates. They showed that the system has a positive global solution. They then linearised the stochastic delayed
SIR model and studied the exponential mean square stability of the linearised system with and without delay.
A more recent paper written by Shrestha et al. [24] looked at a different aspect of the SIRS model. They developed a
new dynamic message-passing (DMP) algorithm, namely rDMP for recurrent epidemic models such as the SIRS model on
networks. They have shown that the rDMP algorithm provides a good approximation to the results obtained from Monte
Carlo simulation, its accuracy is often better than the pair approximation and that rDMP is more user friendly.
The well-known SIS (Susceptible–Infectious–Susceptible) epidemic model [12,25] is used to model diseases which
do not develop immunity once infected individuals recover, for example gonorrhea, Ref. [25], meningitis [12] and
pneumococcus [26,27]. Inspired by the work done by Takeuchi et al. [9], Gray et al. [2] introduced the effect of telegraph
noise into this model using a finite state MC. They established the conditions required for almost sure (a.s.) extinction and
persistence for their solution to the stochastic SIS model with finite state Markovian Switching (MS).
Consequently motivated by Refs. [2,9], in this paper, we will extend the results given in Ref. [2] to a more complicated
three-dimensional SIRS epidemic model as well as the SIR epidemic model. Note that Wei et al. [28] also looked at the
stochastic SIR model under regime switching, but the two models and the results obtained differ. In Ref. [28] they considered
an SIR epidemic model with a nonlinear incidence term different to ours and different per capita death rates for susceptible,
infectious and removed individuals. Similarly to Ref. [10], Wei et al. simultaneously introduced the effect of white noise into
the deterministic model in the form of Brownian motion as well as telegraph noise using a continuous time finite state space
MC. The results that Wei et al. obtain are different to ours and although their model is more general some of their results are
obtained under quite restrictive conditions, whereas our results are not. So their model and results are different. Our model
follows the same basic idea as in Refs. [2,9]. We have a group of deterministic SIS models with different parameter values
corresponding to different environmental regimes. The switching between regimes occurs according to a continuous time
finite state space MC. As far as we know, although there have been various types of work done on the SIRS and SIR models,
ours is the first paper that gives a detailed analysis on the effect that telegraph noise has on the SIRS model, and thus we
have filled a gap in the existing literature.
The paper is organised as follows. In Section 2, we will introduce the MS SIRS epidemic model. A recap of some of the
fundamental concepts of finite state MCs will also be given. In Section 3, the existence of a unique nonnegative solution
will be proven. In Section 4, we will look at the conditions needed for extinction for the MS SIRS model. In Section 5, we
will obtain the conditions needed for persistence. In Section 6, by using the Lyapunov Theorem, we examine the persistence
conditions on the stochastic SIRS model. In Section 7, we will summarise our results and explain how the results for the
SIR model are a special case of the SIRS model. Numerical simulations are produced throughout the paper to support our
theoretical results.

2. MS SIRS epidemic model

Unless stated otherwise, we let (Ω , F , {Ft }t ≥0 , P) be a complete probability space with filtration {Ft }t ≥0 satisfying the
usual conditions (i.e. it is increasing and right continuous while F0 contains all P-null sets). Let us consider the following
deterministic SIRS epidemic model:

dS (t )
= −β I (t )S (t ) + µN − µS (t ) + υ R(t ),
dt
dI (t )
= β I (t )S (t ) − (µ + γ )I (t ), (2.1)
dt
dR(t )
= γ I (t ) − µR(t ) − υ R(t ),
dt
 D. Greenhalgh et al. / Physica A 462 (2016) 684–704 687

where S , I and R denote respectively the number of susceptible, infectious and recovered individuals. N is the total size of
the population, β is the disease transmission coefficient and β = λ/N where λ is the disease contact rate, that is the rate at
which susceptibles come into potentially infectious contact with infecteds. µ is the per capita birth and death rate and γ is
the rate at which an infected becomes cured and thus moves to the recovery group. υ is the rate of loss of immunity. Those
individuals who lose immunity immediately re-enter into the susceptible class.
Note that in this paper S , I and R denote respectively the absolute numbers of individuals in the population as opposed
to the proportions. The total population size remains constant so if s, i and r denote the fractions of individuals in each of
these categories they satisfy the differential equations:

ds(t )
= −β Ni(t )s(t ) + µ − µs(t ) + υ r (t ),
dt
di(t )
= β Ni(t )s(t ) − (µ + γ )i(t ), (2.2)
dt
dr (t )
= γ i(t ) − µr (t ) − υ r (t ).
dt
So the equations have the same functional form, the only differences are that the per capita disease transmission
coefficient in the model using absolute numbers is replaced by the per capita disease contact rate in the model using
proportions and the population input term changes from µN in the absolute numbers model to µ in the proportional model.
Whilst there are many mathematical epidemic models which use the proportions of individuals in each class there are
also many papers that use numbers, for example Refs. [29–31,23,32,33,2,17,34–37] and many more. It is more natural to
use absolute numbers rather than proportions as the differential equations are usually derived by considering the changes
in absolute numbers and then converted to proportions.
As our results are concerned with persistence of solutions and lower and upper bounds for the lim supremum and
lim infimum of the variables there is no qualitative difference between the results obtained for our model using absolute
numbers and the results which would have been obtained from a model using proportions. It is straightforward to convert
the results for our model with absolute numbers into those for the proportional model and vice-versa.
Also note that our models include births and deaths. Again there is a dichotomy of models in the literature. Whilst there
are some models that do not include births and deaths these are generally appropriate only for modelling relatively short
epidemic outbreaks. If epidemics are modelled over a long period of time births and deaths must be included. Whilst there
are some standard models that do not include births and deaths many standard epidemic models do [29,33,38,12,39]. As
many of our results are about the long-term behaviour of the system we feel that it is appropriate to include population
demographics, that is births and deaths, into the model.
We shall now recall some of the fundamental theories of MCs. Let r (t ), t ≥ 0, be a right-continuous MC on the probability
space taking values in finite state space S = {1, 2, . . . , M } with generator Γ = (νij )M ×M defined as

νij δ + o(δ), if i ̸= j,

P{r (t + δ) = j|r (t ) = i} = (2.3)
1 + νii δ + o(δ), if i = j,

where δ > 0, νij ≥ 0 is the transition rate from state i to j for i ̸= j and νii = − 1≤j≤M ,j̸=i νij . Almost every sample path of

r (·) is a right-continuous step function with a finite number of sample jumps in any finite subinterval of R+ = [0, ∞) [40].
There is a sequence {τk }k≥0 of finite-valued Ft -stopping times such that 0 = τ0 < τ1 < · · · < τk → ∞ a.s. and
∞

r (t ) = r (τk )1[τk ,τk+1 ) (t ) (2.4)
k=0

where 1A denotes the indicator function of set A. The switching is memoryless and if r (τk ) = i, the random variable τk+1 − τk
will have an exponential distribution with parameter −νii . In addition, let us define 5 = (π1 , π2 , . . . , πM ) to be the unique
stationary distribution of this MC. If M = 2
ν21 ν12
π1 = and π2 = . (2.5)
ν12 + ν21 ν12 + ν21
Now we will introduce two-state MS into (2.1) which becomes

dS (t )
= −βr (t ) I (t )S (t ) + µr (t ) N − µr (t ) S (t ) + υr (t ) R(t ),
dt
dI (t )
= βr (t ) I (t )S (t ) − (µr (t ) + γr (t ) )I (t ), (2.6)
dt
dR(t )
= γr (t ) I (t ) − µr (t ) R(t ) − υr (t ) R(t ),
dt
where r (t ) is a right-continuous MC with state space S = {1, 2}. We will focus on analysing this model.
688 D. Greenhalgh et al. / Physica A 462 (2016) 684–704

3. Existence of unique nonnegative solution

Theorem 3.1. For any given initial value S (0) = S0 ∈ (0, N ), I (0) = I0 ∈ (0, N ) and R(0) = R0 ∈ (0, N ), there exists a unique
and nonnegative solution for the MS SIRS model (2.6) for all t.
Proof. The proof is straightforward and so we will not discuss this in detail here.

4. Extinction

In this section we will focus on discussing the conditions for extinction for our MS SIRS model (2.6). For the deterministic
SIRS model, the criterion used to determine whether a disease will go extinct or persist is called the basic reproduction
βN
number RD0 = µ+γ . This represents the expected number of secondary infections caused by an infected individual entering
the DFE (e.g. Refs. [29,33,38,12]). If RD0 > 1 then we expect that the disease will persist while RD0 ≤ 1 indicates that it will
die out. We will use another type of threshold to determine whether the disease will die out or persist a.s., namely
π1 β1 N + π2 β2 N
T0S = , (4.1)
π1 (µ1 + γ1 ) + π2 (µ2 + γ2 )
where
• (π1 , π2 ) is the unique stationary distribution,
• β1 is the disease transmission coefficient in state 1,
• β2 is the disease transmission coefficient in state 2,
• µ1 is the per capita birth and death rates in state 1,
• µ2 is the per capita birth and death rates in state 2,
• γ1 is the recovery rate in state 1,
• γ2 is the recovery rate in state 2.
This notation is used by Gray et al. [2] to analyse extinction and persistence for their MS SIS model. By working with the
same threshold we will extend their results to our more complex three-dimensional SIRS model (2.6). Let us recall that r (t )
is a MC with state space S = {1, 2}. If r (t ) = 1, then we are in state 1 and if r (t ) = 2 then we are in state 2.

Proposition 4.1. Let us define αr (t ) = βr (t ) N − µr (t ) − γr (t ) , then we have the following alternative ways of interpreting T0S :
• T0S < 1 ⇔ π1 α1 + π2 α2 < 0,
• T0S = 1 ⇔ π1 α1 + π2 α2 = 0,
• T0S > 1 ⇔ π1 α1 + π2 α2 > 0.
Proof. The proof is straightforward.

Theorem 4.2. If T0S < 1, then for any given initial value (S0 , I0 , R0 ) ∈ (0, N )3 , the solution of the stochastic SIRS epidemic
model (2.6) obeys
(i) lim supt →∞ 1t log(I (t )) ≤ α1 π1 + α2 π2 < 0 a.s.,
(ii) limt →∞ R(t ) = 0 a.s.,
(iii) limt →∞ S (t ) = N a.s.
By Proposition 4.1, we hence conclude that I (t ) tends to zero exponentially and R(t ) tends to zero as t → ∞, thus making
S (t ) tend to N as t → ∞ a.s. In other words, the disease will die out with probability one and the solution will tend to its DFE
(N , 0, 0).
Proof. (i) It is a straightforward modification of the proof of Theorem 4.2 in Ref. [2].
(ii) Suppose that lim supt →∞ R(t ) > 0 on a set Ω1 where P(Ω1 ) = δ for some δ > 0. Then I (t ) → 0 as t → ∞ on a set
Ω2 where P(Ω2 ) ≥ 1 − 2δ . For ω ∈ Ω2 then given ϵ > 0 let us choose ϵ1 small enough so that

ε1 max(γ1 , γ2 ) ε
< , (4.2)
min(µ1 + υ1 , µ2 + υ2 ) 2
where υ1 and υ2 represent the rate of loss of immunity in state 1 and state 2 respectively. The other parameter values are
defined as before.
∃t0 such that for t ≥ t0 , 0 ≤ I (t ) ≤ ϵ1 . By integrating the R(t ) equation in (2.6), we have that for t ≥ t0 ,
 t
R(t ) = R(t0 )e−Q (t ) + e−Q (t ) γr (s)I (s)eQ (s) ds,
t0
 t t
≤ Ne−Q (t ) + γr (s) ε1 e− s (µr (u) +υr (u) )du
ds, (4.3)
t0
 D. Greenhalgh et al. / Physica A 462 (2016) 684–704 689

(µr (s) + υr (s))ds. By choosing t1 ≥ t0 such that for t ≥ t1 , Ne−Q (t ) ≤ 21 ε we have that for t ≥ t1 ,
t
where Q (t ) = t0

ε t

R(t ) ≤ + ε1 max(γ1 , γ2 ) e− min(µ1 +ν1 ,µ2 +ν2 )(t −s) ds. (4.4)
2 t0

Using (4.2) we have that for t ≥ t1 , R(t ) ≤ ε . Hence for ω ∈ Ω2 , lim supt →∞ R(t ) = 0. This is a contradiction. Hence as
t → ∞, R(t ) → 0 a.s.
Theorem 4.2(iii) is obvious by using the fact that S + I + R = N. 
Note that if both α1 < 0 and α2 < 0, then clearly the corresponding RD0,i values for both subsystems (state 1 and state 2)
are less than one, thus both subsystems will die out. However, the readers may wonder what would happen if one subsystem,
say state 1, has α1 < 0 while in state 2 α2 > 0? In other words, one subsystem will go extinct whilst the other will persist.
It turns out that if the time it takes for the MC to switch from state 2 to state 1 is relatively faster than from state 1 to 2, so
that π1 α1 + π2 α2 < 0, then the effect from state 1 will predominate, thus making the overall system die out.
Throughout the paper we shall use numerical simulations to illustrate our results. We shall assume that the unit of time
is one day, and the population sizes are measured in units of one million. We now illustrate Theorem 4.2 using a numerical
example.

Example 4.3. Let us define the parameters to be

µ1 = 0.65, µ2 = 0.10, γ1 = 0.45, γ2 = 0.25, υ1 = 0.15, υ2 = 0.75
β1 = 0.002, β2 = 0.005, ν12 = 0.5, ν21 = 0.8 and N = 100.
By using the definition of αr (t ) in Proposition (2.5) and (4.1), we deduce that α1 = −0.90, α2 = 0.15, π1 = 8/13 and
π2 = 5/13. Thus π1 α1 + π2 α2 = −0.4962 < 0 to four d.p. Similarly, by using Theorem 4.2, we expect that for any initial
value (S (0), I (0), R(0)) ∈ (0, N )3 , the solution to our stochastic SIRS model (2.6) satisfies
1. lim supt →∞ 1t log(I (t )) ≤ −0.4962 < 0 a.s.,
2. limt →∞ R(t ) = 0 a.s.,
3. limt →∞ S (t ) = N a.s.
In other words, the disease will die out a.s.
Again, the numerical simulations produced by using the Euler method support our results in Theorem 4.2, namely the
disease dies out a.s. Note that in this case, α1 < 0 while α2 > 0. The biological meaning of this is that one subsystem will
die out while the other subsystem will persist. Similarly, the numerical simulations were repeated numerous times with
different parameter values and initial values and all support our results (see Fig. 1).

5. Persistence

Apart from extinction, the aspect of persistence of a disease is very important when analysing an epidemic model for
a particular disease. As a result, in this section we will be looking at different types of conditions on persistence for the
MS SIRS model (2.6) when T0S > 1. Note that there are two possible cases that could arise from the condition T0S > 1,
i.e. π1 α1 + π2 α2 > 0, namely:
α α
(a) Both α1 and α2 are positive. Without loss of generality, we will assume that 0 < β1 ≤ β2 .
1 2
α α
(b) One of α1 and α2 is positive. Without loss of generality, we will assume that β1 ≤ 0 < β2 .
1 2

First, we will examine in detail the persistence condition > 1 by looking at the above two cases separately in order to
T0S
give us a better understanding of the persistence results. Before we begin with the main theorems, we will look at another
aspect of persistence which is given by using the uniform persistence theorem (e.g. Refs. [41,42]). We will prove that our
solution I (t ) for our stochastic SIRS model (2.6) under either case for T0S > 1 is uniformly strong persistent. Recall that
αr (t ) = βr (t ) N − µr (t ) − γr (t ) .

Theorem 5.1 (Uniform Strong Persistence). Suppose that I (0) > 0.

α α
Case (a): If 0 < β1 ≤ β2 , ∃ε ′ > 0 independent of the initial conditions such that
1 2

lim inf I (t ) ≥ ε > 0 ′

a.s. (5.1)
t →∞

In other words the MS SIRS model is a.s. uniformly persistent.

α α
Case (b): If β1 < 0 < β2 , given δ1 > 0, ∃ε ′ > 0 such that ∀t1 > 0, I (t ) ≥ ε ′ for some t ≥ t1 on a set Ω1 where
1 2
P(Ω1 ) ≥ 1 − δ1 . To put this another way,
lim inf I (t ) > 0 a.s.
t →∞
690 D. Greenhalgh et al. / Physica A 462 (2016) 684–704

Fig. 1. Numerical simulation for our solution to (2.6) with T0S < 1 and the corresponding MC r (t ) using the parameter values given in Example 4.3 with
initial values S (0) = 20, I (0) = 60, R(0) = 20 and r (0) = 1.

Proof. Case (a): Let us choose ε > 0 small enough such that
α1 min(µ1 + υ1 , µ2 + υ2 )
ε< .
β1 min(µ1 + υ1 , µ2 + υ2 ) + 2 max(γ1 , γ2 )
Suppose that I (t ) < ε for all t ≥ t0 and I (0) > 0. Then from the third equation in (2.6) for t ≥ t0

dR(t )
≤ max(γ1 , γ2 )ε − min(µ1 + υ1 , µ2 + υ2 )R(t ). (5.2)
dt
By integrating (5.2), it is easy to obtain the following expression:

max(γ1 , γ2 )ε
R(t ) ≤ Ne− min(µ1 +υ1 ,µ2 +υ2 )(t −t0 ) + . (5.3)
min(µ1 + υ1 , µ2 + υ2 )
Let us choose t1 > t0 such that for t ≥ t1 , we have
max(γ1 , γ2 )ε
Ne− min(µ1 +υ1 ,µ2 +υ2 )(t −t0 ) ≤ . (5.4)
min(µ1 + υ1 , µ2 + υ2 )
By using (5.4), (5.3) becomes

2 max(γ1 , γ2 )ε
R(t ) ≤ , (5.5)
min(µ1 + υ1 , µ2 + υ2 )
for t ≥ t1 .
We have that
1 dI (t )
= αr (t ) − βr (t ) (I (t ) + R(t )),
I (t ) dt
≥ αr (t ) − βr (t ) (ε + R(t )). (5.6)
By substituting (5.5) into the above equation, we have that

1 dI (t ) 2 max(γ1 , γ2 )
  
≥ min αr − βr ε 1 + = K1 > 0. (5.7)
I (t ) dt r ={1,2} min(µ1 + υ1 , µ2 + υ2 )
 D. Greenhalgh et al. / Physica A 462 (2016) 684–704 691

This implies that I (t ) is an increasing function and it must eventually increase above ε . Moreover, from our argument we
2 max(γ ,γ )ε
know that by time t1 , R(t ) must drop to a level at most min(µ +υ 1,µ2 +υ ) , where
1 1 2 2

max(γ1 , γ2 )ε max(γ1 , γ2 )ε
  
−1
 log , if < N,
min(µ1 + υ1 , µ2 + υ2 ) N min(µ1 + υ1 , µ2 + υ2 ) min(µ1 + υ1 , µ2 + υ2 )


t1 − t0 = (5.8)
max(γ1 , γ2 )ε
0, ≥ N.

 if
min(µ1 + υ1 , µ2 + υ2 )
− max(µ1 +γ1 ,µ2 +γ2 )t1
Furthermore, from the second equation of (2.6) I (t1 ) ≥ I (0)e > 0. For t ≥ t1 , I (t ) ≥ I (t1 )eK1 (t −t1 ) ≥
I (0)e− max(µ1 +γ1 ,µ2 +γ2 )t1 eK1 (t −t1 ) , hence I (t ) must reach level ε by a time at most t2 where

ε
   
1
t2 = t1 + log + max(µ1 + γ1 , µ2 + γ2 )t1 . (5.9)
K1 I (0)
As a result, we have shown that if I (0) < ε , then I (t ) will reach the level ε by at most time t2 . In other words, I (t )
will always be greater than ε at some time in the future provided that we start below it. However it is possible for I (t )
subsequently to go below ε again later. Consequently, we will now assume that I (0) = ε and from the above if I (t ) does go
below ε , it will eventually rise back up again by time at most
 
1
′
t2 = t1 1 +′
max(µ1 + γ1 , µ2 + γ2 ) , (5.10)
K1
where t1′ is defined by (5.8) with t0 = 0.
In general, let us define t ∗ with I (t ∗ ) = ε to be the first time that I (t ) drops beneath ϵ . Then a similar argument as before
will show that for t ≥ t ∗ ,

I (t ) ≥ ε e− max(µ1 +γ1 ,µ2 +γ2 )t2 = ε ′ > 0.

′
(5.11)
So we have shown that our solution I (t ) to (2.6) is uniformly strong persistent in case (a).
Case (b): In this case, we have that α1 < 0, which indicates that RD0,1 < 1 in state 1 while in state 2 we have RD0,2 > 1.
In other words, if we stay in state 1 long enough, I (t ) will tend to 0 thus making our solution (S (t ), I (t ), R(t )) for (2.6) tend
towards the DFE (N , 0, 0). As a result, unlike in case (a), the uniform strong persistence result will not hold for all the domain
as there will be a region where it is possible for I (t ) to approach 0 arbitrarily closely with a small but non-zero probability.
However, we can make the probability of that happening as small as we want.
Choose ε small enough so that
2 max(γ1 , γ2 )
  
π1 α1 − β1 ε 1 +
min(µ1 + υ1 , µ2 + υ2 )
2 max(γ1 , γ2 )
  
+ π2 α2 − β2 ε 1 + > K2 > 0. (5.12)
min(µ1 + υ1 , µ2 + υ2 )
Now suppose that lim inft →∞ I (t ) = 0 on a set Ω1 where P(Ω1 ) = δ1 > 0. By the ergodic theory of the MC, ∃ T independent
of the initial state such that on a set Ω2 where P(Ω2 ) ≥ 1 − 2δ for t ≥ T ,

2 max(γ1 , γ2 )
 t  
1
lim αr (s) − βr (s) ε 1 + ds
t →∞ 0t min(µ1 + ν1 , µ2 + ν2 )
2 max(γ1 , γ2 )
  
= π1 α1 − β1 ε 1 +
min(µ1 + υ1 , µ2 + υ2 )
2 max(γ1 , γ2 )
  
+ π2 α2 − β2 ε 1 + > K2 . (5.13)
min(µ1 + υ1 , µ2 + υ2 )
Consider any ω ∈ Ω1 ∩ Ω2 . Suppose that ∃ t0 (ω) such that I (t ) ≤ ε for all t ≥ t0 (ω). Similarly to case (a) we have that R(t )
2 max(γ ,γ )ε
falls beneath a level at most min(µ +υ 1,µ2 +υ ) from time t1 (ω) = t0 (ω) + t1′ onwards. Again similarly to case (a), I (t1 ) > 0.
1 1 2 2
By integrating (5.6) and substituting R(t ) by its upper bound given by (5.5), we have that for t ≥ t1 (ω),

I (t ) 2 max(γ1 , γ2 )
   t  
log ≥ αr (s) − βr (s) ε 1 + ds. (5.14)
I (t1 ) t1 min(µ1 + υ1 , µ2 + υ2 )
Hence using (5.13) for t ≥ t1 (ω) + T ,

I (t )
 
1
lim log ≥ K2 > 0, (5.15)
t →∞ t − t1 (ω) I (t1 (ω))
692 D. Greenhalgh et al. / Physica A 462 (2016) 684–704

so for t ≥ t1 (ω) + T , I (t ) ≥ I (t1 )eK2 (t −t1 ) . In other words, from time t1 + T onwards, I (t ) is bounded below by an
increasing unbounded function and thus we have a contradiction and I (t ) must rise above the level ε by a time at most
max(t2 (ω), t1 (ω) + T ) where ε = I (t1 (ω))eK2 (t2 (ω)−t1 (ω)) .
Starting at max(t2 (ω), t1 (ω) + T ), ∃ t3 (ω) > max(t2 (ω), t1 (ω) + T ) with I (t3 (ω)) = ϵ . Moreover arguing as previously
every time that I (t ) drops beneath ε it must rise up again to this level by time at most t4′ where
 
1
t4′ = (t1′ + T ) 1 + max(µ1 + γ1 , µ2 + γ2 ) > t1′ .
K2
Then similarly to (a) we have that lim inft →∞ I (t ) ≥ ε ′ for some ε ′ > 0, contradicting ω ∈ Ω1 . This completes the proof of
Theorem 5.1. 
Let us now look at more conditions on persistence for our MS SIRS model (2.6).

Theorem 5.2. If T0S > 1, then for any given initial value (S (0), I (0), R(0)) ∈ (0, N )3 , then the solution S (t ) of the stochastic
SIRS model has the properties that:
π α +π α
(a) lim inft →∞ S (t ) ≤ N − π1 β1 +π2 β2 a.s.,
1 1 2 2
π α +π α
(b) lim supt →∞ S (t ) ≥ N − π1 β1 +π2 β2 a.s.
1 1 2 2
π α +π α
In other words, the number of susceptibles will reach the neighbourhood of the level N − π1 β1 +π2 β2 infinitely many times a.s.
1 1 2 2

Proof. Case (a): Assume the statement given in Theorem 5.2(a) is not true, then ∃ε > 0 sufficiently small such that
P(Ω1 ) > 0 where
π1 α1 + π2 α2
 
Ω1 = ω ∈ Ω : lim inf S (t ) > N − + 2ε .
t →∞ π1 β1 + π2 β2
In addition, by the ergodic theory of the MC, we have that P(Ω2 ) = 1 where for any ω ∈ Ω2 ,
π1 α1 + π2 α2
 t  
1
lim αr (s) − βr (s) − ε ds
t →∞ t 0 π1 β1 + π2 β2
π1 α1 + π2 α2 π1 α1 + π2 α2
     
= π1 α1 − β1 −ε + π2 α2 − β2 −ε ,
π1 β1 + π2 β2 π1 β1 + π2 β2
= (π1 β1 + π2 β2 )ε. (5.16)
Now consider any ω ∈ Ω1 ∩ Ω2 . Then there is a positive number T = T (ω) such that
π1 α1 + π2 α2
S (t ) ≥ N − + ε, ∀t ≥ T (ω), (5.17)
π1 β1 + π2 β2
which we can easily rearrange to get
π1 α1 + π2 α2
I (t ) + R(t ) ≤ − ε, ∀t ≥ T (ω). (5.18)
π1 β1 + π2 β2
By integrating (5.6) and using (5.18), we have that for all t ≥ T (ω),
T
π1 α1 + π2 α2
  t  
log(I (t )) ≥ log(I (0)) + [αr (s) − βr (s) (I (s) + R(s))]ds + αr (s) − βr (s) −ε ds.
0 T π1 β1 + π2 β2
Dividing both sides by t and letting t → ∞, we could simplify the above expression to
1
lim inf log(I (t )) ≥ (π1 β1 + π2 β2 )ε > 0 (5.19)
t →∞ t
by using (5.16). So I (t ) → ∞ as t → ∞, which clearly contradicts (5.18). As a result, it is obvious that our assumption at
the beginning is false and Theorem 5.2(a) follows.
Case (b): As above we will assume that there exists ε > 0 sufficiently small such that P(Ω3 ) > 0 where
π1 α1 + π2 α2
 
Ω3 = ω ∈ Ω : lim sup S (t ) < N − − 2ε .
t →∞ π1 β1 + π2 β2
Consider any ω ∈ Ω2 ∩ Ω3 . Then there is a positive number T = T (ω) such that
π1 α1 + π2 α2
S (t ) ≤ N − − ε, ∀t ≥ T (ω). (5.20)
π1 β1 + π2 β2
The remainder of the proof follows the same lines as (a). 
 D. Greenhalgh et al. / Physica A 462 (2016) 684–704 693

As a result we have proved that the number of susceptibles will persist and it will reach the neighbourhood of the level
π α +π α
N − π1 β1 +π2 β2 infinitely many times a.s.
1 1 2 2
Before we look at the persistence theorem for I (t ) we need the following lemma:

Lemma 5.3. Given ε1 > 0,

(a) If I (t ) ≥ ξ for t ≥ t0 , ∃t1 ≥ t0 such that for t ≥ t1 ,
γ1 γ2
 
R(t ) ≥ ξ min , (1 − ε1 ).
µ1 + υ1 µ2 + υ2
(b) If I (t ) ≤ ξ for t ≥ t0 , ∃t1 ≥ t0 such that for t ≥ t1 ,
γ1 γ2
 
R(t ) ≤ ξ max , (1 + ε1 ).
µ1 + υ1 µ2 + υ2
Proof. Case (a): Let us define a sequence of stopping times t0 = τ0 < τ1 < · · · < τm < t where τm+1 is interpreted as t.
dR(t )
Then for the case I (t ) ≥ ξ , the equation dt for our stochastic SIRS model defined in (2.6) gives:

d 
R(t )eF (t ) ≥ γr (t ) ξ eF (t ) ,

(5.21)
dt
m

where F (t ) = (µr (τk ) + υr (τk ) )(τk+1 − τk ). (5.22)
k=0

By integrating Eq. (5.21), replacing the term F (t ) with (5.22) and some rearranging, we deduce that:
 t
F (t )
R(t )e − R(t0 ) ≥ γr (s) ξ exp (µr (τ0 ) + υr (τ0 ) )(τ1 − τ0 ) + · · · + (µr (τm′ ) + υr (τm′ ) )(s − τm′ ) ds,
 
t0

where t0 = τ0 < τ1 < · · · < τm′ ≤ s . . . ≤ τm ≤ t ,

m 
 τk+1
γr (s) ξ exp (µr (τ0 ) + υr (τ0 ) )(τ1 − τ0 ) + · · · + (µr (τk ) + υr (τk ) )(s − τk ) ds,
 
=
k=0 τk
γr (τ0 ) γr (τm )
ξ eF (τ1 ) − eF (τ0 ) + · · · + ξ eF (t ) − eF (τm ) ,
   
=
µr (τ0 ) + υr (τ0 ) µr (τm ) + υr (τm )
where eF (τ0 ) = 1,
γ1 γ2
 
≥ ξ min , (eF (t ) − 1). (5.23)
µ1 + υ 1 µ2 + υ 2
As a result,

γ1 γ2
 
R(t ) ≥ R(t0 )e−F (t ) + ξ min , (1 − e−F (t ) ). (5.24)
µ1 + υ 1 µ2 + υ 2
Given ε1 > 0 by choosing t large enough, we have that for t ≥ t1 ,

γ1 γ2
 
R(t ) ≥ ξ min , (1 − ε1 ). (5.25)
µ1 + υ 1 µ2 + υ 2
We have thus completed the proof for Lemma 5.3(a). The proof for case (b) follows similarly. 

Theorem 5.4. If T0S > 1, then for any given initial value (S (0), I (0), R(0)) ∈ (0, N )3 , the solution I (t ) of the stochastic SIRS
model has the properties that:
 
π1 α1 +π2 α2
(a) lim inft →∞ I (t ) ≤ π1 β1 +π2 β2
 1
γ1 γ2
 a.s.,
1+min µ +υ ,
  1 1 µ2 +υ2
π1 α1 +π2 α2
(b) lim supt →∞ I (t ) ≥ π1 β1 +π2 β2
 1
γ1 γ2
 a.s.
1+max µ +υ ,
1 1 µ2 +υ2

So given ϵ > 0 the number of infectives will enter between the levels

π α + π α  1 π α + π α  1
1 1 2 2 1 1 2 2
−ϵ  and +ϵ
π1 β1 + π2 β2 π β π β
  
γ1 γ2 + γ1 γ2
1 + max µ +υ , µ +υ 1 1 2 2 1 + min µ +υ , µ2 +υ2
1 1 2 2 1 1

infinitely often a.s.

694 D. Greenhalgh et al. / Physica A 462 (2016) 684–704

Proof. Case (a): Suppose that the assertion is false. Then there exists ε > 0 such that P(Ω5 ) > 0 where
 
π1 α1 + π2 α2
 
1
 
Ω5 = ω ∈ Ω : lim inf I (t ) >  + 2ε .
π1 β1 + π2 β2 1 + min

γ1
 t →∞
, γ2
µ1 +υ1 µ2 +υ2


Now by considering any ω ∈ Ω5 , there is a positive number T = T (ω) such that

π1 α1 + π2 α2
 
1
I (t ) ≥  + ε, (5.26)
π1 β1 + π2 β2

γ γ
1
1 + min µ +υ
1 1
, µ2 +υ
2
2

for all t ≥ T (ω). From Lemma 5.3(i), given ε1 > 0 and I (t ) ≥ ξ + ε , ∃T1 (ω) ≥ T (ω) such that for t ≥ T1 (ω) ≥ T (ω),

γ1 γ2
 
R(t ) ≥ (ξ + ε) min , (1 − ε1 ), (5.27)
µ1 + υ 1 µ2 + υ 2
 
π1 α1 +π2 α2
where ξ = π1 β1 +π2 β2
 1
γ1 γ2
. By using S (t ) + I (t ) + R(t ) = N, (5.27) becomes
1+min µ +υ ,
1 1 µ2 +υ2

γ1 γ2
   
S (t ) ≤ N − (ξ + ε) 1 + min , (1 − ε1 ) , (5.28)
µ1 + υ 1 µ2 + υ 2
whence
γ1 γ2
   
lim sup S (t ) ≤ N − (ξ + ε) 1 + min , (1 − ε1 ) . (5.29)
t →∞ µ1 + υ 1 µ2 + υ 2
Now let ε1 → 0. By using Theorem 5.2 we arrive at the following contradiction

γ1 γ2
  
0 ≤ −ε 1 + min , , (5.30)
µ1 + υ1 µ2 + υ2
and we must therefore have
π1 α1 + π2 α2
 
1
lim inf I (t ) ≤  a.s. (5.31)
π1 β1 + π2 β2

γ γ
t →∞ 1
1 + min µ +υ
1 1
, µ2 +υ
2
2

Case (b): Similarly, we will assume that there exists ε > 0 sufficiently small such that P(Ω6 ) > 0 where
 
π1 α1 + π2 α2
 
1
 
Ω6 = ω ∈ Ω : lim sup I (t ) <  − 2ε .
π1 β1 + π2 β2

γ γ
 t →∞
1
1
1 + max µ +υ
1
, µ2 +υ
2
2


By using a similar method as in Case (a), but this time using Lemma 5.3(b), the result follows easily. 

Theorem 5.5. If T0S > 1, then for any given initial value (S (0), I (0), R(0)) ∈ (0, N )3 , then the solution R(t ) of the stochastic
SIRS model (2.6) has the properties that:
(a) lim inft →∞ R(t ) > 0 a.s.,
N max(γ ,γ2 )
(b) lim supt →∞ R(t ) < max(γ ,γ )+min(µ1 +υ ,µ +υ )
< N a.s.
1 2 1 1 2 2
In other words, the limiting value of the number of recovered individuals will be strictly positive and will not ultimately exceed
N max(γ1 ,γ2 )
max(γ ,γ )+min(µ +υ ,µ +υ )
a.s.
1 2 1 1 2 2

Proof. Case (a): If lim inft →∞ R(t ) = 0 on a set Ω1 where P(Ω1 ) ≥ δ > 0 then by Theorem 5.1, ∃ε > 0 and t0 such
that I (t ) ≥ ε > 0 for t ≥ t0 on a set Ω2 where P(Ω2 ) ≥ 1 − 2δ > 0. By Lemma 5.3 ∃ε ′ > 0 and t1 > t0 such that for
t ≥ t1 , R(t ) ≥ ε ′ > 0 on Ω2 . This is a contradiction proving the result.
Case (b): Let us choose
N max(γ1 , γ2 )
ξ= < N.
max(γ1 , γ2 ) + min(µ1 + υ1 , µ2 + υ2 )
From (2.6),
dR(t )
≤ max(γ1 , γ2 )(N − R) − min(µ1 + υ1 , µ2 + υ2 )R(t ), (5.32)
dt
= N max(γ1 , γ2 ) − [max(γ1 , γ2 ) + min(µ1 + υ1 , µ2 + υ2 )]R(t ). (5.33)
The result of Theorem 5.5(b) follows. 
 D. Greenhalgh et al. / Physica A 462 (2016) 684–704 695

We will continue to investigate persistence by looking at the two cases that could possibly arise from T0S > 1.

α α
Theorem 5.6. Assume that T0S > 1 and let I (0) ∈ (0, N ) be arbitrary. If β1 ≤ 0 < β2 , then the following statements hold a.s.:
1 2
  
α γ γ
(i) lim inft →∞ S (t ) ≥ N − β2 1
1 + max µ +υ , µ2 +υ
2
.
2 1 1 2
α
(ii) lim supt →∞ I (t ) ≤ β2 .
2
 
α2 γ γ
(iii) lim supt →∞ R(t ) ≤ β2
1
max µ +υ , µ2 +υ
2
.
1 1 2

α
Proof. Note that I (t ) > 0 for all t. Suppose that lim supt →∞ I (t ) > β2 . Then using Theorem 5.4(a), ∃t1 and t2 with t1 < t2 ,
2
α
such that β2 < I (t1 ) < I (t2 ) and I (t ) is strictly monotonic increasing in [t1 , t2 ]. Let us now choose t3 ∈ (t1 , t2 ), not a jump
2
dI (t )
point of the MC such that dt
> 0. For r (t ) = 1 and r (t ) = 2, from (2.6):

dI (t ) 

1
 = αi − βi (I (t3 ) + R(t3 )) < 0, for i = 1, 2. (5.34)
I (t3 ) dt t
3

dI (t )
Clearly, for both states we have that dt
< 0 which is a contradiction. Theorem 5.6(ii) follows. Subsequently, by using
Lemma 5.3(b), we have that

α2 γ1 γ2
 
lim sup R(t ) ≤ max , , (5.35)
t →∞ β2 µ1 + υ1 µ2 + υ2
whence by using the fact that S (t ) + I (t ) + R(t ) = N, we obtain the desired result that

α2 γ1 γ2
  
lim inf S (t ) ≥ N − 1 + max , . 
t →∞ β2 µ1 + υ1 µ2 + υ2

Example 5.7. Let us now define the parameters to be

µ1 = 0.65, µ2 = 0.40, γ1 = 0.45, γ2 = 0.20, υ1 = 0.15, υ2 = 0.75

β1 = 0.009, β2 = 0.012, ν12 = 0.5, ν21 = 0.8 and N = 100.
We see that α1 = −0.2, α2 = 0.60, π1 = 8/13 and π2 = 5/13, where clearly π1 α1 + π2 α2 = 0.1077 > 0 to four d.p. From
Theorem 5.6, we expect that for any initial value (S (0), I (0), R(0)) ∈ (0, N )3 , the solution to (2.6) satisfies:
  
α γ γ
(i) lim inft →∞ S (t ) ≥ N − β2 1 + max µ +υ
1
, µ2 +υ
2
= 21.875,
2 1 1 2
α
(ii) lim supt →∞ I (t ) ≤ β2 = 50,
2
 
α γ γ
(iii) lim supt →∞ R(t ) ≤ β2 max µ +υ
1
, µ2 +υ
2
= 28.125,
2 1 1 2

to three d.p.
Again, the numerical simulations generated by the Euler method illustrated in Fig. 2 support our results in Theorem 5.6.

α α
Theorem 5.8. (a) Assume that T0S > 1 and let I (0) ∈ (0, N ) be arbitrary. If 0 < β1 ≤ β2 , then the following statements hold
1 2
a.s.:
  
α γ γ
(i) lim inft →∞ S (t ) ≥ N − β2 1
1 + max µ +υ , µ2 +υ
2
.
2 1 1 2
α
(ii) lim supt →∞ I (t ) ≤ β2 .
2
 
α γ γ
(iii) lim supt →∞ R(t ) ≤ β2 max µ +υ
1
, µ2 +υ
2
.
2 1 1 2

(b) If I (0) > 0 under the same conditions the following statements hold a.s.:
     
α α γ1 γ2 γ1 γ2
(i) lim supt →∞ S (t ) ≤ N − β1 − β2 max µ +υ , µ +υ
× 1 + min ,
µ1 +υ1 µ2 +υ2
.
1
2 1 1
 2 2
γ1
α α
(ii) lim inft →∞ I (t ) ≥ β1 − β2 max µ +υ , µ γ+υ
2
.
1 2
1 1 2 2   
α1 α2 γ1 γ2 γ1 γ2
(iii) lim inft →∞ R(t ) ≥ β − β max µ +υ , µ +υ × min µ +υ , µ +υ
.
1 2 1 1 2 2 1 1 2 2
696 D. Greenhalgh et al. / Physica A 462 (2016) 684–704

Fig. 2. Numerical simulation for our solution to (2.6) T0S > 1 and its corresponding MC r (t ) using the parameter values given in Example 5.7 with initial
values S (0) = 60, I (0) = 20, R(0) = 20 and r (0) = 1.

Proof. The proof for case (a) follows as in Theorem 5.6. In order to prove Theorem 5.8(b), without loss of generality we may
assume that

α1 α2 γ1 γ2
 
> max , .
β1 β2 µ1 + υ1 µ2 + υ2

Suppose that Theorem 5.8(bii) is false and choose ε > 0 such that

α1 α2 γ1 γ2
 
lim inf I (t ) < − max , −ε (5.36)
t →∞ β1 β2 µ1 + υ 1 µ2 + υ 2

on a set Ω1 where P(Ω1 ) = δ1 > 0. Moreover by Theorems 5.4(b) and 5.8(aii)

π1 α1 + π2 α2
 
1
lim sup I (t ) ≥
π β π β
 
+ γ1 γ2
t →∞ 1 1 2 2 1 + max µ +υ
1 1
, µ 2 +υ 2

α2 γ1 γ2
 
and lim sup R(t ) ≤ max , ,
t →∞ β2 µ1 + υ 1 µ2 + υ 2
on a set Ω2 where P(Ω2 ) = 1.
For ω ∈ Ω1 ∩ Ω2 , ∃t4 (ω) such that for t ≥ t4 ,

α2 γ1 γ2
 
R(t ) < max , + ε.
β2 µ1 + υ1 µ2 + υ2
α1 γ1 γ2
  
Also lim sup I (t ) > 1 − max , ,
t →∞ β1 µ1 + υ 1 µ2 + υ 2
α1 α2 γ1 γ2
 
> − max , .
β1 β2 µ1 + υ1 µ2 + υ2
 D. Greenhalgh et al. / Physica A 462 (2016) 684–704 697

Fig. 3. Numerical simulation for our solution (S (t ), I (t ), R(t )) to (2.6) T0S > 1 and its corresponding MC r (t ) using the parameter values given in
Example 5.9 with initial values S (0) = 20, I (0) = 45, R(0) = 35 and r (0) = 1.

Hence from (5.36) there must exist some t5 and t6 where t4 < t5 < t6 such that

α1 α2 γ1 γ2
 
I (t6 ) < I (t5 ) < − max , − ε,
β1 β2 µ1 + υ1 µ2 + υ2
α2 γ1 γ2
  
≤ 1 − max , − ε, (5.37)
β2 µ1 + υ 1 µ2 + υ 2
and I (t ) is strictly monotonic decreasing in [t5 , t6 ].
dI (t )
Let us now choose t7 ∈ (t5 , t6 ), not a jump point of the MC, such that  < 0. Similar to the proof for Theorem 5.6,

dt t7
dI (t )
 > 0 which again is a contradiction proving Theorem 5.8(bii). Again, by using Lemma 5.3 and

it is easy to show that dt t7
that S (t ) + I (t ) + R(t ) = N, we obtained the required results (bi)–(biii). 
α α
Therefore for the case 0 < β1 ≤ β2 , we have obtained both an upper and lower bound for our solution (S , I , R) for (2.6),
2 2
α α
which is a better result than in the case β1 ≤ 0 < β2 .
1 2

Example 5.9. Let us define the parameter values to be

µ1 = 0.85, µ2 = 0.50, γ1 = 0.55, γ2 = 0.20, υ1 = 0.15, υ2 = 0.75

β1 = 0.02, β2 = 0.012, ν12 = 0.5, ν21 = 0.8 and N = 100.
By using the definition of αr (t ) defined in Proposition (2.5) and (4.1), we deduce that α1 = 0.6, α2 = 0.5, π1 = 8/13 and
π2 = 5/13, where clearly π1 α1 + π2 α2 = 0.5615 > 0 to four d.p. By substituting the appropriate parameter values into
Theorem 5.8, we would expect that for any initial value (S (0), I (0), R(0)) ∈ (0, N )3 ,

(a) 35.4167 ≤ lim inft →∞ S (t ) ≤ lim supt →∞ S (t ) ≤ 91.7833,

(b) 7.0833 ≤ lim inft →∞ I (t ) ≤ lim supt →∞ I (t ) ≤ 41.6667,
(c) 1.1333 ≤ lim inft →∞ R(t ) ≤ lim supt →∞ R(t ) ≤ 22.9167,

to four d.p. a.s. This implies that regardless of whatever the initial values, the solution (S (t ), I (t ), R(t )) asymptotically
approaches the appropriate region above.
Once again, we could conclude from Fig. 3 that the numerical simulations support our results proved in Theorem 5.8. The
numerical simulations were repeated many times with various initial values and the same conclusion was obtained.
698 D. Greenhalgh et al. / Physica A 462 (2016) 684–704

In the next section, we will continue to investigate persistence, but we will be using Lyapunov stability (e.g. Refs.
[18,35,20]) as well as Theorem 5.1, to obtain results on the convergence of the solution (S , I , R) to its corresponding endemic
α α α α
and disease-free equilibria in each of state 1 and state 2 under the persistence conditions 0 < β1 ≤ β2 and β1 ≤ 0 < β2 .
1 2 1 2

6. Lyapunov stability

When analysing the behaviour of a dynamical system, one of the significant aspects would be the stability of the solution.
There are various types of stability, but the most important one is the stability of a solution near its equilibrium point, in
other words will the solution converge to its equilibrium point or will it diverge? This aspect of stability could be discussed
by using a Lyapunov Theorem, which is what we shall look at here. By using the results from Theorem 5.1, we have obtained
some results which further enhance our understanding. It is easy to see that the DFE is (N , 0, 0) while the endemic equilibria
for state 1 and 2 are:
N
Si∗ = , (6.1)
R0,i
µi + υ i αi µi + υ i
   
1
Ii∗ = 1− N = , (6.2)
µi + υi + γi R0,i βi µi + υi + γi
γi γi αi
   
1
R∗i = 1− N = , (6.3)
µi + υi + γi R0,i µi + υi + γi βi
βN
where RD0,i = µ +γ
i
is the basic reproduction number when the MC is in state i for i = 1, 2.
i i

α α
Theorem 6.1. Assume that T0S > 1 and 0 < β1 ≤ β2 and let (S (0), I (0), R(0)) ∈ (0, N )3 be arbitrary and let the switching
1 2
times of the MC be 0 = τ0 < τ1 < · · · < τk where τk → ∞ as k → ∞. Define the Lyapunov function to be:

βi
 
I
Vi (x) = I − Ii∗ − Ii∗ log + (R − R∗i )2 , (6.4)
Ii∗ 2γi

where x = (S (t ), I (t ), R(t )), for i = 1, 2.

Note that by considering the Taylor series expansion about I = Ii∗ for ϵ small enough, say ε ≤ ε1 then

(I − Ii∗ )2
 
1 I
∗ (I − Ii ) ≤ I − Ii − Ii log ,
∗ 2 ∗ ∗
≤ (6.5)
4Ii Ii∗ Ii∗

in (Ii∗ − ε, Ii∗ + ε), for i = 1, 2.

For any ε ≤ ε1 sufficiently small, the Lyapunov function (6.4) for our MS SIRS model has the properties that:

ε2 β1 2(µ1 + υ1 + γ1 )
  
P lim inf V1 (t ) < 1+ ≥ e−ν12 T1 (ε) , (6.6)
t →∞ 2(µ1 + υ1 ) α1
and

ε2 β2 2(µ2 + υ2 + γ2 )
  
P lim inf V2 (t ) < 1+ ≥ e−ν21 T2 (ε) , (6.7)
t →∞ 2(µ2 + υ2 ) α2
where T1 (ε) = ε2Wβ > 0 and T2 (ε) = ε2Wβ > 0 for some constant W .
1 2

Proof. By differentiating the Lyapunov function (6.4), we have that

βi
 
dVi
= (βi S − µi − γi )(I − Ii∗ ) + (R − R∗i )(γi (I − Ii∗ ) − (µi + υi )(R − R∗i )). (6.8)
dt γi
Now by substituting (µi + γi ) = βi (N − Ii∗ − R∗i ) and S = N − I − R, consequently, (6.8) becomes

dVi (µi + υi )βi

= −βi (I − Ii∗ )2 − (R − R∗i )2 < 0. (6.9)
dt γi
Thus, Vi (x) ≥ 0 and V̇i (x) ≤ 0 with equality if and only if I = Ii∗ and R = R∗i . If there is no switching then Vi (x) is a Lyapunov
function and the endemic equilibria (Si∗ , Ii∗ , R∗i ) given by (6.1)–(6.3) are globally asymptotically stable, i.e. S → Si∗ , I → Ii∗
and R → R∗i as t → ∞, whatever the initial condition.
 D. Greenhalgh et al. / Physica A 462 (2016) 684–704 699

We shall now prove the results with switching time involved. The proof will split into two parts, corresponding to the
Lyapunov functions for state 1 and state 2. First we will show that the result holds in state 1. By Theorem 5.1, ∃t1 , W < ∞
such that for t ≥ t1 ,

βi
 
I
Vi (x) = I − Ii∗ − Ii∗ log + (R − R∗i )2 ≤ W < ∞, (6.10)
Ii∗ 2γi
and max(V1 (t ), V2 (t )) ≤ W .
Define a stopping time
σ1 = inf {t ≥ t1 : r (t ) = 1} .
Clearly, P(σ1 < ∞) = 1, and by the right-continuity of the MC, r (σ1 ) = 1. Define

V1 (σ1 )
T1′ (ε) = < ∞, (6.11)
ε 2 β1
and note that
V1 (σ1 ) W
T1′ (ε) = ≤ T1 (ε) = a.s. (6.12)
ε 2 β1 ε 2 β1
The probability that the MC will not jump to state 2 before σ1 + T1′ (ε) is

P(Ω1 ) = e−ν12 T1 (ε) ,

′

where Ω1 = ω : r (σ1 + t ) = 1, for all t ∈ [0, T1′ (ε)] . Consider any ω ∈ Ω1 on [σ1 , σ1 + T1′ (ε)] and suppose that
 

(µ1 + υ1 )β1
− β1 (I − I1∗ )2 − (R − R∗1 )2 ≤ −ε 2 β1 , (6.13)
γ1
in this region which by rearranging implies that

(µ1 + υ1 )
(I − I1∗ )2 + (R − R∗1 )2 ≥ ε 2 > 0, (6.14)
γ1
for t ∈ [σ1 , σ1 + T1′ (ε)]. As a result, for t ∈ [σ1 , σ1 + T1′ (ε)], (6.9) becomes

dVi
≤ −ε 2 β1 . (6.15)
dt
Thus, after integrating:

0 ≤ V1 (σ1 + T1′ (ε)) ≤ V1 (σ1 ) − ε 2 β1 (T1′ (ε)), (6.16)

from which by substituting T1 (ε) by its definition in (6.11),
′

V1 (σ1 + T1′ (ε)) = 0. (6.17)

However, if we recall the Lyapunov function given by (6.4), it is equal to zero if and only if I (σ1 + T1 (ε)) = I1∗ and ′

R(σ1 + T1′ (ε)) = R∗1 . This clearly contradicts (6.14) for t ∈ [σ1 , σ1 + T1′ (ε)]. Thus, we must have instead

(µ1 + υ1 )β1
β1 (I − I1∗ )2 + (R − R∗1 )2 < ε2 β1 , (6.18)
γ1
for some s ∈ [σ1 , σ1 + T1′ (ε)]. Note that at time s,

(I − I1∗ )2 ε2 ε 2 γ1
< , and (R − R∗1 )2 < . (6.19)
∗
I1 ∗
I1 µ1 + υ 1
Therefore, if ε ≤ ε1 , then by using (6.5)

(I − I1∗ )2 ε2
 
I
0 ≤ I − I1∗ − I1∗ log ∗ ≤ ≤ . (6.20)
I1 I∗ I1∗
By using (6.19) and (6.20), the Lyapunov function (6.4) at time s is bounded above by

β1
 
1
V1 (s) < ε 2 + . (6.21)
I1∗ 2(µ1 + υ1 )
700 D. Greenhalgh et al. / Physica A 462 (2016) 684–704
 
µ +υ α1
Recall from (6.2) that I1∗ = µ +υ
1 1
β1
hence (6.21) becomes
1 1 +γ1

(µ1 + υ1 + γ1 )β1 β1
 
V1 ( s ) < ε 2 + . (6.22)
(µ1 + υ1 )α1 2(µ1 + υ1 )
Consequently, if T ≥ 0,
ε2 β1 2(µ1 + υ1 + γ1 )
  
≥ P(Ω1 ) = e−ν12 T1 (ε) ,
′
P inf V1 (t ) < +1
T ≤t <∞ 2(µ1 + υ1 ) α1
≥ e−ν12 T1 (ε) , (6.23)
where T1 (ε) = ε2 β defined as before.
W
1
Note that
ε 2 β1 2(µ1 + υ1 + γ1 )
  
lim inf V1 (t ) < +1
t →∞ 2(µ1 + υ1 ) α1
ε 2 β1 2(µ1 + υ1 + γ1 )
   
= inf V1 (t ) < +1 . (6.24)
0<T <∞
T ≤t <∞ 2(µ1 + υ1 ) α1
By letting T → ∞ in (6.23), we have obtained (6.6). (6.7) is proven similarly. 
Theorem 6.1 shows that our solution (S (t ), I (t ), R(t )) can approach either endemic equilibria (Si∗ , Ii∗ , R∗i ) arbitrarily
closely with strictly positive probability.

Corollary 6.2. If ε ≤ ε1 , then

   
2α1 γ1 2(µ1 + υ1 + γ1 )

P lim inf max{|S − S1 |, |I − I1 |, |R − R1 |} < ε
∗ ∗ ∗
4+ + +1
t →∞ µ1 + υ1 + γ1 µ1 + υ 1 α1

≥ e−ν12 T1 (ε) , (6.25)

and
   
2α2 γ2 2(µ2 + υ2 + γ2 )

P lim inf max{|S − S2 |, |I − I2 |, |R − R2 |} < ε
∗ ∗ ∗
4+ + +1
t →∞ µ2 + υ2 + γ2 µ2 + υ 2 α2

≥ e−ν21 T2 (ε) , (6.26)

where T1 (ε) = ε2 β , T2 (ε) = ε2 β , and ε1 is defined as in Theorem 6.1. Recall that αi = βi N − µi γi .
W W
1 2

Proof. (i) We shall begin by looking at state 1. Recall from (6.5) that for I ∈ (Ii∗ − ε, Ii∗ + ε) and ε < ε1

(I − Ii∗ )2
 
1 I
∗ (I − Ii ) ≤ Ii − Ii − Ii log ,
∗ 2 ∗ ∗
≤ (6.27)
4Ii Ii∗ Ii∗
which implies that if (6.18) holds for t ∈ [σ1 , σ1 + T1′ (ε)] then for some s ∈ [σ1 , σ1 + T1′ (ε)],

ε 2 β1 2(µ1 + υ1 + γ1 )
 
1
∗ (I − I1 ) ≤ V1 (s) ≤ +1 .
∗ 2
(6.28)
4I1 2(µ1 + υ1 ) α1
By rearranging the above expression, and taking the square root we deduce that

2α1
|I − I1 | ≤ ε 4 +
∗
. (6.29)
µ1 + υ1 + γ1
Recall again from (6.4) that

β1
 
I
V1 (s) = I − I1∗ − I1∗ log + (R − R∗1 )2 , (6.30)
I1∗ 2γ1
which by using (6.28) and some simple rearrangement we have that

γ1 2(µ1 + υ1 + γ1 )
 
|R(s) − R1 | ≤ ε∗
+1 . (6.31)
µ1 + υ 1 α1
 D. Greenhalgh et al. / Physica A 462 (2016) 684–704 701

By using S (s) = N − I (s) − R(s) and S1∗ = N − I1∗ − R∗1 , then

  
2α1 γ1 2(µ1 + υ1 + γ1 )

|S (s) − S1 | ≤ ε
∗
4+ + +1 . (6.32)
µ1 + υ1 + γ1 µ1 + υ 1 α1

Arguing as in the proof of Theorem 6.1 it is easy to see that (6.25) holds.
(ii) The proof for state 2 follows similarly. 

Corollary 6.2 shows similarly to Theorem 6.1, but using the Euclidean metric instead of the metric induced by the
Lyapunov function, that the solution (S (t ), I (t ), R(t )) can approach either endemic equilibrium (Si∗ , Ii∗ , R∗i ) arbitrarily closely
with strictly positive probability.
α α
In Theorem 6.1 and Corollary 6.2 we have been focusing on analysing the persistence condition where 0 < β1 ≤ β2 by
1 2
using Lyapunov stability. We will now complete the results on persistence by obtaining results on the convergence of the
α1 α2
solution (S , I , R) to its corresponding disease-free and endemic equilibria under the condition β ≤ 0 < β .
1 2

α α
Theorem 6.3. Assume that T0S > 1 (namely π1 α1 + π2 α2 > 0) and β1 ≤ 0 < β2 . Let (S0 , I0 , R0 ) ∈ (0, N ) be arbitrary. Then
1 2
the solution to (2.6) has the properties that
(i) If ε > 0, then
2γ1
  
P lim inf max(|N − S |, |I |, |R|) ≤ ε 1+ ≥ e−ν12 T1 (ε) , (6.33)
t →∞ µ1 + υ 1
where T1 (ε) = t̄1 (ε) + t̄2 (ε) and t̄1 (ε) and t̄2 (ε) are defined as:

γ1 ε (µ1 + υ1 )N
  
−1
, ≥ ε,
  
 1 log N , log if

if N ≥ ε, µ1 + υ 1 (µ1 + υ1 )N γ1


t̄1 (ε) = β1 ε ε and t̄2 (ε) = (6.34)
(µ1 + υ1 )N
0, if N < ε. 0 , < ε,
 
 if
γ1
respectively.
(ii) If ε > 0 is small enough such that
2 max(γ1 , γ2 ) 2 max(γ1 , γ2 )
     
π1 α1 − β1 2ε 1 + + π2 α2 − β2 2ε 1 + > 0, (6.35)
min(µ1 + υ1 , µ2 + υ2 ) min(µ1 + υ1 , µ2 + υ2 )
ε 2 β2 2(µ2 + υ2 + γ2 )
  
then P lim inf V2 (t ) ≤ 1+ ≥ e−ν21 T2 (ε) , (6.36)
t →∞ 2(µ2 + υ2 ) α2
W (ε)
  β
where T2 (ε) = β ε2 and W (ε) = max N − I2∗ − I2∗ log IN∗ , ε − I2∗ − I2∗ log Iε∗  + 2γ2 N 2 < ∞. Vi (x) denotes the Lyapunov
   
2 2 2 2
function which is defined as in (6.4) for i = 1, 2.

Proof. (i) Suppose that ε > 0. Define a stopping time such that

σ1 = inf {t ≥ 0 : r (t ) = 1} .
Clearly, P(σ1 < ∞) = 1 and by the right-continuity of the MC, r (σ1 ) = 1. The probability that the MC will not jump to state
2 before σ1 + T1 (ε) is

P(Ω1 ) = e−ν12 T1 (ε) ,

where Ω1 = {ω : r (σ1 + t ) = 1, for all t ∈ [0, T1 (ε)]}. Consider any ω ∈ Ω1 on [0, T1 (ε)], it is easy to see that
dI (t )
≤ −β1 I (t )2 ≤ −β1 ε I (t ),
dt
provided I ≥ ε > 0, which after integration becomes

I (σ1 + t ) ≤ I (σ1 )e−β1 εt ≤ Ne−β1 εt . (6.37)

If N ≥ ε then (6.37) shows that by time t 1 (ε), I (t ) must drop to a level at most ε where
 
1 N
t̄1 (ε) = log . (6.38)
β1 ε ε
702 D. Greenhalgh et al. / Physica A 462 (2016) 684–704

On the other hand if N < ε then I (0) < N < ε . Arguing as in Theorem 5.1, we know that for t ≥ t̄1 (ε) + t̄2 (ε),

2γ1 ε
R(σ1 + t ) ≤ , (6.39)
µ1 + υ 1
γ1 ε (µ1 + υ1 )N
  
−1
 log , if ε ≤ ,
µ1 + υ 1 N (µ1 + υ1 ) γ1


where t̄2 (ε) = (6.40)
(µ1 + υ1 )N
0, if ε > .


γ1
Hence for t ≥ t̄1 (ε) + t̄2 (ε), we have that

2γ1
 
|N − S (σ1 + t )| = I (σ1 + t ) + R(σ1 + t ) ≤ ε 1 + . (6.41)
µ1 + υ 1
Thus we could see from (6.41) that

2γ1
 
max N − S (σ1 + t̄1 (ε) + t̄2 (ε)), I (σ1 + t̄1 (ε) + t̄2 (ε)), R(σ1 + t̄1 (ε) + t̄2 (ε)) ≤ ε 1 + .
 
(6.42)
µ1 + υ 1
As this is true for each ω ∈ Ω1 , we have that

P max |N − S (σ1 + t̄1 (ε) + t̄2 (ε))|, I (σ1 + t̄1 (ε) + t̄2 (ε)), R(σ1 + t̄1 (ε) + t̄2 (ε))
 

2γ1
 
≤ε 1+ ≥ e−ν12 T1 (ε) , (6.43)
µ1 + υ 1
where T1 (ε) = t̄1 (ε) + t̄2 (ε). Consequently, if T ≥ 0, then

2γ1
  
P inf max (|N − S (t )|, I (t ), R(t )) ≤ ε 1+ ≥ e−ν12 T1 (ε) . (6.44)
T ≤t <∞ µ1 + υ 1
Theorem 6.3(i) follows by arguing as in the proof of Theorem 6.1.
(ii) Recall that (6.35) holds, which is inequality (5.12) with ε replaced by 2ε . Note also that when r (t ) = 1, RD0,1 ≤ 1 and
also dt1 < 0. Hence, if Ω denotes the whole sample space, given ω ∈ Ω and t3 (ω) > 0, for t ≥ t3 (ω), I (t ) must rise up and
dI

over the level ε at some time t4 (ω) > t3 (ω). So ∃t5 (ω) > t3 (ω) with I (t5 (ω)) = ε and r (t5 (ω)) = 2. Also V2 (t5 (ω)) ≤ W (ε)
where V2 (t ) denotes the Lyapunov function in state 2 given by (6.4) in Theorem 6.1 and W (ε) is a constant.
Now arguing as in the proof of Theorem 6.1 define a new stopping time

t5 (ω) = inf{t ≥ 0 : r (t5 (ω)) = 2, I (t5 (ω)) = ε}

we will have the required result namely,

ε 2 β2 2(µ2 + υ2 + γ2 )
  
P lim inf V2 (t ) < +1 ≥ e−ν21 T2 (ε) , (6.45)
t →∞ 2(µ2 + υ2 ) α2
W (ε) β
where T2 (ε) = ε2 β and W (ε) = max N − I2∗ − I2∗ log IN∗ , ε − I2∗ − I2∗ log Iε∗  + 2γ2 N 2 < ∞, which could be easily
     
2 2 2 2
derived from (6.10). 
In this theorem, we have obtained interesting probabilistic results on the convergence of the solution (S (t ), I (t ), R(t ))
of the stochastic SIRS model (2.6) to its corresponding disease-free and endemic equilibria.

Corollary 6.4. If ε < ε1 , then:

   
2α2 γ2 2(µ2 + υ2 + γ2 )

P lim inf max{|S − S2∗ |, |I − I2∗ |, |R − R∗2 |}ε 4+ + +1
t →∞ µ2 + υ2 + γ2 µ2 + γ2 α2

≥ e−ν21 T2 (ε) , (6.46)

W (ε)
where T2 (ε) = ε2 β and ε1 is defined as in Theorem 6.1.
2

Proof. Similar to the proof for Corollary 6.2. 

 D. Greenhalgh et al. / Physica A 462 (2016) 684–704 703

6.1. T0S = 1 case

So far, we have looked into great detail on the dynamic behaviour of (2.6) under the thresholds T0S < 1 and T0S > 1.
The reader may ask what about the case when T0S = 1? We are unable to prove analytically the behaviour of our solution
(S (t ), I (t ), R(t )) in this situation. However numerical simulations indicated that the disease would always ultimately die
out whatever the initial conditions.

7. Summary and discussion

There are many environmental factors that could affect the behaviour of a population system such as the availability of
food and temperature [43]. As discussed in the introduction some previous work on ecological models with telegraph noise
has been done by Takeuchi et al. [9] who obtained results on boundedness or convergence to the equilibrium of trajectories
of the Lotka–Volterra model with telegraph noise. Li et al. [10] have obtained results on the behaviour of an n-dimensional
Lotka–Volterra model under regime switching. Motivated by Refs. [2,9] we have examined the effect of environmental noise
on a more complicated model, the SIRS model, by using the concept of MS to include telegraph noise to give the MS SIRS
model (2.6). In our model we have linked several deterministic models in different environmental regimes using a MC. We
have obtained the conditions needed for almost surely extinction and persistence using the threshold T0S which was also
used in Ref. [2]. In Theorem 4.2, we showed that if T0S < 1 then the disease will go extinct almost surely. On the other hand if
T0S > 1, then the disease will persist almost surely (Theorems 5.2, 5.4 and 5.5). In Theorems 5.6 and 5.8 we obtained two sets
α α α α
of persistence conditions for the two possible cases in which T0S > 1, namely β1 ≤ 0 < β2 and 0 < β1 ≤ β2 . Furthermore
1 2 1 2
by using the uniform strong persistence result for I (t ), (Theorem 5.1) and the Lyapunov stability theorem, we managed to
obtain probabilistic results on convergence of our solution to the disease-free and endemic equilibria in Section 6. Numerical
simulations were produced to support and illustrate our theoretical results.
Note that it is easy to see that the two-state MS Susceptible–Infectious–Removed (SIR) model is a special case of the MS
SIRS model (2.6). In fact, we could easily derive the corresponding extinction and persistence results for the MS SIR model
by setting υr (t ) , in (2.6) to zero. Furthermore, the results obtained for the two-state MC can be easily extended into a finite
state MC with state space S = {1, 2, . . . , M }, similarly to the corresponding extension in Ref. [2].
Remember that we have chosen to model the absolute numbers of individuals in each category as opposed to the
proportions. However the results which we have obtained have all been concerned with persistence of solutions and upper
and lower limits for the lim supremum and lim infimum of the variables. The results do not depend on our choice to model
the absolute numbers of individuals in the population rather than the proportions. There is no essential difference between
the models in the two formats. It is straightforward to convert the results for the model in one format into the results for
the model with the other format.
Note also that we chose to include births and deaths into the model because we felt that this was appropriate as we
were considering modelling the disease over a long timescale. If we exclude births and deaths from the model it will be
appropriate only for short term disease outbreaks. For the SIRS model (υ1 > 0 or υ2 > 0) all of the results go through
(just set µ1 = µ2 = 0). Our proofs break down if µ1 = µ2 = 0, so that the population is closed with no births and
deaths, and either or both of υ1 and υ2 are zero, so that at least one of the models is the closed population SIR model. If
µ1 = µ2 = υ1 = υ2 = 0 then it is clear that Theorem 4.2 is not true as if T0S < 1 in this case I (t ) → 0 but S (t ) decreases
monotonically to a limiting value S ∗ < N whilst R(t ) increases monotonically to a limiting value R∗ > 0.
If we consider the SIRS epidemic model without MS then provided that there is some mechanism for generating new
susceptibles, whether that is through births and deaths in the population, or through immune individuals losing immunity
and returning to the susceptible class, then the qualitative behaviour of the population is the same. There is a threshold
βN
value R0 = µ+γ . If R0 ≤ 1 then the disease ultimately dies out. If R0 > 1 then there is a unique endemic equilibrium which
the system ultimately approaches [12].
On the other hand for the SIR model in a closed population, with no births and deaths where immunity is permanent
βN
then there is a qualitative difference in the behaviour. Here R0 = γ . The initial infective replacement number is defined as
the expected number of cases generated by each separate infective individual at the start of the epidemic. So this is
S (0) β S (0)
R0 = .
N γ
If this number exceeds unity then the number of infectious individuals rises up to a limit and then drops down to zero. If
this number is less than or equal to one then the number of infectious decreases to zero. The disease transmission ceases
because the infective replacement number (the expected number of cases generated by a single infective individual) is less
than one when the number of susceptibles is small, however the limiting number of susceptibles is strictly positive [12].
So even without MS there is a fundamental qualitative difference between an SIRS epidemic model that has a mechanism
for introducing new susceptibles and one that does not. The results in this paper are true provided that all of the individual
SIRS models corresponding to different states of the MC have some mechanism for generating new susceptibles, whether
that is births and deaths, or immunity being temporary. The behaviour of SIRS models with MS where at least one of the
704 D. Greenhalgh et al. / Physica A 462 (2016) 684–704

individual SIRS models is a closed population (without births and deaths) SIR model remains an interesting open question
for further study.

Acknowledgements

Yanfeng Liang is grateful to the EPSRC for a studentship (EPSRC DTP Grant Reference no. EP/K503174/1) to support
this work. We are grateful to Drs A. Gray and J. Pan for passing on a program that was developed to use in the numerical
simulations. David Greenhalgh is grateful to the Leverhulme Trust for support from a Leverhulme Research Fellowship (RF-
2015-88). Xuerong Mao is grateful to the Leverhulme Trust (RF-2015-385) and the Royal Society of London (IE131408) for
their financial support.

References

[1] N.H. Du, V.H. Sam, Dynamics of a stochastic Lotka–Volterra model perturbed by white noise, J. Math. Anal. Appl. 324 (1) (2006) 82–97.
[2] A. Gray, D. Greenhalgh, X. Mao, J. Pan, The SIS epidemic model with Markovian switching, J. Math. Anal. Appl. 394 (2) (2012) 496–516.
[3] X. Mao, Stochastic Differential Equations and Applications, second ed., Horwood, Chichester, UK, 2008.
[4] X. Mao, S. Sabanis, E. Renshaw, Asymptotic behaviour of the stochastic Lotka-Volterra model, J. Math. Anal. Appl. 287 (1) (2003) 141–156.
[5] X. Mao, G. Marion, E. Renshaw, Environmental Brownian noise suppresses explosions in population dynamics, Stochastic Process. Appl. 97 (1) (2004)
95–110.
[6] N.H. Du, R. Kon, K. Sato, Y. Takeuchi, Dynamical behavior of Lotka–Volterra competition systems: Non-autonomous bistable case and the effect of
telegraph noise, J. Comput. Appl. Math. 170 (2004) 399–422.
[7] M. Slatkin, The dynamics of a population in a Markovian environment, Ecology 9 (2) (1978) 249–256.
[8] Y. Takeuchi, Global Dynamical Properties of Lotka–Volterra Systems, World Scientific, 1996.
[9] Y. Takeuchi, N.H. Du, N.T. Hieu, K. Sato, Evolution of predator–prey systems described by a Lotka–Volterra equation under random environment,
J. Math. Anal. Appl. 323 (2006) 938–957.
[10] X. Li, D. Jiang, X. Mao, Population dynamical behavior of Lotka-Volterra system under regime switching, J. Comput. Appl. Math. 232 (2009) 427–448.
[11] W.O. Kermack, A.G. McKendrick, A contribution to the mathematical theory of epidemics, Proc. R. Soc. Lond. Ser. A 115 (1927) 700–721.
[12] H.W. Hethcote, Qualitative analyses of communicable disease models, Math. Biosci. 28 (1976) 335–356.
[13] E. Tornatore, S.M. Buccellato, P. Vetro, Stability of a stochastic SIR system, Physica A 354 (2005) 111–126.
[14] Q. Lu, Stability of SIRS system with random perturbations, Physica A 388 (18) (2009) 3,677–3,686.
[15] Q. Yang, D. Jiang, N. Shi, C. Ji, The ergodicity and extinction of stochastically perturbed SIR and SEIR epidemic models with saturated incidence, J. Math.
Anal. Appl. 388 (1) (2012) 248–271.
[16] Y. Zhao, D. Jiang, The threshold of a stochastic SIRS epidemic model with saturated incidence, Appl. Math. Lett. 34 (2014) 90–93.
[17] 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.
[18] S.N. O’Regan, T.C. Kelly, A. Korobeinikov, M.J.A. O’Callaghan, A.V. Pokrovskii, Lyapunov functions for SIR and SIRS epidemic models, Appl. Math. Lett.
23 (4) (2010) 446–448.
[19] A. Korobeinikov, Lyapunov functions and global stability for SIR and SIRS epidemiological models with non-linear transmission, Bull. Math. Biol. 68
(2006) 615–626.
[20] C. Vargas-de-León, On the global stability of SIS, SIR and SIRS epidemic models with standard incidence, Chaos Solitons Fractals 44 (12) (2011)
1,106–1,110.
[21] X. Liu, P. Stechlinski, Pulse and constant control schemes for epidemic models with seasonality, Nonlinear Anal. RWA 12 (2) (2011) 931–946.
[22] I. Nasell, Stochastic models of some endemic infections, Math. Biosci. 179 (2002) 1–19.
[23] G. Chen, T. Li, Stability of stochastic delayed SIR model, Stoch. Dyn. 9 (2) (2009) 231–252.
[24] M. Shrestha, S.V. Scarpino, C. Moore, A message-passing approach for recurrent-state epidemic models on networks, Phys. Rev. E 92 (2) (2015)
0220821.
[25] H.W. Hethcote, J.A. Yorke, Gonorrhea Transmission Dynamics and Control, in: Lecture Notes in Biomathematics, vol. 56, Springer-Verlag, Berlin, 1984.
[26] K.E. Lamb, D. Greenhalgh, C. Robertson, A simple mathematical model for genetic effects in pneumococcal carriage and transmission, J. Comput. Appl.
Math. 235 (7) (2011) 1,812–1,818.
[27] M. Lipsitch, Vaccination against colonising bacteria with multiple serotypes, Proc. Natl. Acad. Sci. USA 94 (1997) 6,571–6,576.
[28] Q. Wei, Z. Xiong, F. Wang, Dynamics of a stochastic SIR model under regime switching, J. Inf. Comput. Sci. 10 (2013) 2,727–2,734.
[29] R.M. Anderson, R.M. May, Infectious Diseases of Humans, Oxford University Press, Oxford, 1992.
[30] L.J.S. Allen, An introduction to stochastic epidemic models, in: F. Brauer, P. van den Driessche, J. Wu (Eds.), Mathematical Epidemiology, in: Lecture
Notes in Biomathematics, Mathematical Biosciences Subseries, vol. 1,945, Springer-Verlag, Berlin, 2008.
[31] N.T.J. Bailey, The Mathematical Theory of Infectious Diseases and its Applications, Griffin, London, 1975.
[32] N. Dalal, D. Greenhalgh, X. Mao, A stochastic model of AIDS and condom use, J. Math. Anal. Appl. 325 (2007) 36–53.
[33] O. Diekmann, J.A.P. Heesterbeek, Mathematical Epidemiology of Infectious Diseases, Vol. 146, Wiley, Chichester, 2000.
[34] H.W. Hethcote, P. van den Driessche, An SIS epidemic model with variable population size and a delay, J. Math. Biol. 34 (1995) 177–194.
[35] A. Korobeinikov, G.C. Wake, Lyapunov functions and global stability for SIR, SIRS, and SIS epidemiological models, Appl. Math. Lett. 15 (2002) 955–960.
[36] A. Lahrouz, A. Settati, Asymptotic properties of switching diffusion epidemic model with varying population size, Appl. Math. Comput. 219 (2013)
11,134–11,148.
[37] I. Nasell, Extinction and Quasi-Stationarity in the Stochastic Logistic SIS Model, in: Lecture Notes in Mathematics, vol. 2,022, Springer-Verlag, 2011.
[38] H.W. Hethcote, The mathematics of infectious diseases, SIAM Rev. 42 (2000) 599–653.
[39] M.J. Keeling, P. Rohani, Modelling Infectious Diseases in Humans and Animals, Princeton University Press, Princeton, New Jersey, 2008.
[40] W.J. Anderson, Continuous-Time Markov Chains, Springer-Verlag, Berlin-Heidelberg, 1991.
[41] G.J. Butler, P. Waltman, Persistence in dynamical systems, J. Differential Equations 63 (1986) 255–263.
[42] P. Waltman, A brief survey of persistence in dynamical systems, in: Delay Differential Equations and Dynamical Systems, in: Lecture Notes in
Mathematics, vol. 1,475, Springer Berlin, Heidelberg, 1991, pp. 31–40.
[43] J.W. Tang, The effect of environmental parameters on the survival of airborne infectious agents, J. R. Soc. Interface 6 (2009) S737–S746.