Expected Time to Extinction of SIS Epidemic Model Using Quasy–
Stationary Distribution

arXiv:1803.06096v1 [math.ST] 16 Mar 2018

Kurnia Susvitasari1
1

School of Mathematical Science, University of Nottingham, UK
ARTICLE HISTORY
Compiled March 19, 2018
ABSTRACT
We study that the breakdown of epidemic depends on some parameters, that is
expressed in epidemic reproduction ratio number. It is noted that when R0 exceeds
1, the stochastic model have two different results. But, eventually the extinction
will be reached even though the major epidemic occurs. The question is how long
this process will reach extinction. In this paper, we will focus on the Markovian
process of SIS model when major epidemic occurs. Using the approximation of quasi–
stationary distribution, the expected mean time of extinction only occurs when the
process is one step away from being extinct. Combining the theorm from Ethier
and Kurtz, we use CLT to find the approximation of this quasi distribution and
successfully determine the asymptotic mean time to extinction of SIS model without
demography.
KEYWORDS
quasi–stationary distribution, central limit theorem, SIS

1. Introduction
According to the threshold theorem, epidemic can only occur if the initial number
of susceptibles is larger than some critical values, which depends on the parameters
under cosideration [4] . Usually, it is expressed in terms of epidemic reproductive ratio
number, R0 .
Susvitasari [22] showed that both deterministic and stochastic models performed
similar results when R0 ≤ 1, i.e. the disease–free stage in the epidemic. But then,
when R0 > 1, the deterministic and stochastic models had different interpretations.
In the deterministic model, both SIS and SIR showed an outbreak of the disease, and
after some time t, the disease persisted and reached equilibrium stage, i.e. endemic.
The stochastic model, on the other hand, had different interpretation. There are two
possible outcomes of this approach. First, the infection may die out in the first cycle.
If it did, it would happen very quickly since the time of the disease removed must less
than the time of infectee–susceptible contact. Second, if it survives the first cycle, the
outbreak was likely to occur, but after some time t, it would reach equilibrium as the
deterministic version. In fact, the stochastic model would mimic the deterministic path
and scattered randomly around its equilibrium. Furthermore, by letting population size
K. Susvitasari. Email: susvitasari@sci.ui.ac.id

be large and ignoring the initial value of infectious individuals, as t → ∞, the empirical
distribution of n−1 Yn (t), where Yn (t) is the number of infectious individuals at time t
with n population size followed normal distribution in endemic phase as seen in Figure
2.
In this paper, we will focus on the endemic stage, the stage where R0 > 1 and the
outbreak occurs in SIS epidemic model without demography to determine the expected
mean time to extinction of this model.

2. The Birth and Death Process in SIS
Suppose that {Yn (t) : t ≥ 0} denotes the Markovian process that represents the
number of infected individuals in the population at time t. Then we define Sn =
{0, 1, . . . , n} as the state space of the process, where state zero is an absorbing state.
Figure 1 is the illustration of this process.

Figure 1. Birth–death of SIS without demography transition process

Consider the case where R0 > 1 and the outbreak occurs. Eventually, given a
necessary large n and some time t, the process will not visit state zero. So, let us
redefine the birth–death process before by only considering the transient states. Let
{Ŷn (t) : t ≥ 0} be the modified process of {Yn (t) : t ≥ 0} with state space Sn0 =
Sn \ {0}. The intuitive explaination is that whenever the SIS process goes extinct,
it immediately will restart again in state 1. It is not hard to show that the process
{Ŷn (t) : t ≥ 0} is time–reversible Markov chain. Therefore, there exist colections of
equilibrium distributions πj for all j ∈ Sn0 that sums to unity and satisfy the detailed–
balanced conditions:
λ
(1) (n − 1) π̂i = 2γ π̂2 ,
n


λ
λ
i = 2, 3, . . . , n−1,
(2) i(n − i) + γi π̂i = (i−1)(n−i+1) π̂i−1 +(i+1)γ π̂i+1
n
n
λ
(3) nγ π̂n = (n − 1) π̂n−1 ,
n
X
with constraint
π̂i = 1. Thus the solution is
i∈Sn0

1
π̂i =
i



λ
nγ

i−1

(n − 1)!
π̂1
(n − i)!

i = 1, 2, . . . , n.

(1)

Therefore,

π̂1 =

 i−1
n
X
(n − 1)!
λ
i(n − i)! nγ
i=1

2

!−1
.

(2)

(n)

Now, let Text denote the time to extinction in SIS epidemic model. Recall that the
process {Ŷn (t) : t ≥ 0} is time-reversible on state Sn0 . Then, the time to extinction
(n)
Text is almost surely finite since the state Sn0 is visited finitely often.
Suppose that {Ŷn (t) : t ≥ 0} is in equilibrium, then the mean time to extinction
has intensity π̂1 times the rate of the process getting absorbed to state {0}.



(n)
E Text = E(interarival time of jump 1 → 0)
 i−1
n
1 X (n − 1)!
λ
1
=
.
=
γ π̂1
γ
i(n − i)! nγ

(3)

i=1



log (1 − λ)
(n)
As n → ∞, E Text → −
when we set R0 < 1. Ball et al. [6] showed the
λ
approximation of mean time extinction when R0 < 1, R0 = 1, and R0 > 1.
In this paper, we will focus on the time to extinction when the process reaches
endemic–stage, i.e. when R0 > 1 and the process survives the first cycle of extinction
using quasi– statioary distribution and property of central limit theorm.

3. Expected Time of Extinction of SIS in Endemic–Stage
Before we start this section, we will introduce first the quasi–stationary distribution.
3.1. Quasi–Stationary Distribution
Recall the process {Yn (t) : t ≥ 0} with finite state space Sn , where 0 is an absorbing
state. For j ∈ T = Sn \ {0}, let
qi,j (t) = P (Yn (t) = j | Yn (t) > 0, Yn (0) = i)
Pi,j (t)
i = 1, 2, . . . , n.
=
1 − Pi,0 (t)

(4)

Since we concerns on the distribution of the process in endemic stage, we wish to find
q̃j = lim qi,j (t).
t→∞

We define A and T as sets of absorbing and transient state. Suppose


Q=

QAA
0
QT A QT T



where Q is transition rate matrix, defined as Q = (qij )i,j∈Sn . Let W = diagn×n (π̂).
Since W 1/2 QT T W −1/2 is symmetry, according to spectral decomposition theorem,
W 1/2 QT T W −1/2 =

n
X

λi ui uTi ,

i=1

where λ1 , λ2 , . . . , λn are eigenvalues of W 1/2 QT T W −1/2 and u1 , u2 , . . . , un are corre-

3

sponding orthonormal vectors of right-eigenvectors. Therefore,
QT T =

n
X

λ i E i ⇒ QT T · t =

i=1

n
X

λi t · E i ⇒ exp(QT T t) = exp

i=1

n
X

!
λi t · E i

,

(5)

i=1



W−1/2 ui uTi W1/2 and we can see in Appendix ......... that
n
X
exp(QT T t) =
exp(λi t)E i . Recall that forward Kolmogorov equation, P(t) must
where Ei =

i=1

satisfy P0 (t) = P(t)Q, where P(0) = In×n with solution P(t) = exp(Qt). So, for
i, j ∈ T ,
Pi,j (t) = P {Yn (t) = j | Yn (0) = i}
" n
#
X
 Q t
= e T T ij =
eλk t E k .
k=1

(6)

ij

According to Darroch and Seneta [12] there exists a simple, real eigenvalue λ1 <
0 such that all other eigenvalues of W 1/2 QT T W −1/2 have real part less than λ1 .
Therefore, eq. (4) becomes
qi,j (t) =

eλ1 t (E 1 )ij + eλ2 t (E 2 )ij + . . . + eλn t (E n )ij
n
X
{eλ1 t (E 1 )il + eλ2 t (E 2 )il + . . . + eλn t (E n )il }
l=1

=

(E 1 )ij + e(λ2 −λ1 )t (E 2 )ij + . . . + e(λn −λ1 )t (E n )ij
.
n
X
{(E 1 )il + e(λ2 −λ1 )t (E 2 )il + . . . + e(λn −λ1 )t (E n )il }

(7)

l=1

(E1 )ij
Consequently, as t → ∞, q̃j = Pn
.

 l=1 (E1 )il 


Note that E 1 = W −1/2 u1 uT1 W 1/2 and it turns out that uT1 W 1/2 is left
eigenvector of QT T (see Appendix B). In the stationary scenario, when t → ∞, each
row of probability transition matrix P (t) has nearly similar entries. So, it follows
that the quasi-distribution q̃ is given by the left-eigenvector of QT T corresponding to
n
X
eigenvalue λ1 with constraint q̃1 =
q̃k = 1.
k=1

3.2. Expected Time to Extinction by Quasi–Stationary Distribution
Suppose
that
(Yn (t) | Yn (t) > 0, Yn (0) = i ∼ q̃)
and
let
TQ
=
inf {t ≥ 0; Yn (t) = 0, Yn (t) ∼ q̃}. Since two events (infection and recovery) in
the epidemic model follows the Poisson process, there is zero probability that two or
more events occur at the same time. So, the extinction only occurs when the process

4

is one step away from being extinct. Therefore,
P (not extinct at time t) = P (TQ > t)
n
X
=
lim P (Yn (t) = j | Yn (0) = i)
=

=

j=1
n
X
j=1
n
X

t→∞

q̃j · P (Yn (t) > 0)
q̃j · eQT T t = q̃ · eQT T t · 1.

(8)

j=1

Recall that q̃ is the left-eigenvector of QT T , corresponding to the eigenvalue λ1 . So,
q̃ · QT T = λ1 q̃

⇒

q̃ · QkT T = λk1 q̃

⇒

q̃ · eQT T t = eλ1 t q̃.

Hence, the equation (8) becomes
P (not extinct at time t) = eλ1 t .

(9)

Since the process is Poisson process, then the inter arrival time must follow exponential distribution. Therefore, in eq. (9), TQ follows exponential distribustion with
intensity −λ1 . But, recall the memoryless property of the Poisson process and that the
extinction only occurs if the process is one step away from being extinct. Intuitively, it
also means that the extinction happens with intensity q̃1 times the rate of the process
to jump into absorption state 0, which is γ. On the other hand, λ1 = −γ q̃1 (note that
λ1 < 0).
Proof. Using the fact that q̃ is the left-eigenvector of QT T ,
q̃ · QT T = λ1 q̃ ⇒ q̃ · QT T · 1 = λ1 (q̃1) ⇒ −γ · q̃1 = λ1 .
| {z }
|{z}
=1

(−γ,0,...,0)T

Therefore, the mean time to extinction of the SIS process is
E (TQ ) =

1
.
γ q̃1

(10)

where q̃1 = (q̃)1 .
The quasi–stationary is a very powerful approximation when we let population size
large enough. Unfortunately, consider that by using this method, we need to capture
all the possible transition rate for matrix Q. Considering the large n, the size of matrix
Q will also big and leads to inefficient computation. Therefore, we should approximate
the distribution of q̃1 .
Now, suppose that we only concern on the size of infectious individuals in population
when the outbreak occurs, i.e. {Yn (t)|Yn (t) > 0, Yn (t) ∼ q̃}. By ignoring the initial
value of infectives, the empirical distribution of n−1 Yn (t) follows normal distribution
with mean converges to equilibrium point in deterministic model as seen in Figure
5

2. It turns out that by applying Central Limit Theorem, we could approximate the
quantity of q̃1 .

4. Central Limit Theorem
We define {X n (t); t ≥ 0} to be the CTMC on Zd with finite state space and
P finite number of possible transition l ∈ ∆ ⊆ Zd in such a way that βl : Zd → [0, ∞), l∈∆ βl (x) <
P
∞ for all x ∈ Zd . Then let F (x) = l∈∆ lβl (x) as the drift function of the process
{X n (t) : t ≥ 0}.
Suppose that F (x) is a continuous function and differentiable in Zd lattice such
that ∂F (x) exists and is continuous and let G(x) be d × d matrix function defined as
G(x) =

X

lT lβl (x).

l∈∆

Further, let Φ(t, s), where 0 ≤ s ≤ t, be the solution of the matrix differential equation
∂
Φ(t, s) = ∂F (x(t)) · Φ(t, s)
∂t

(11)

where Φ(s, s) = I.
Theorem 4.1. (Ethier and Kurtz) Suppose that


Z t
√
X n (0)
(1) lim n
− x0 = v 0 is a constant, where x(t) = x0 +
F (x(u))du,
n→∞
n
0
(2) βl and ∂F (x) are continuous functions.
Then, for t ≥ 0,
√


n


X n (t)
D
− x(t) −−−−−→ V (t)
weakly
n

as n → ∞, where {V (t) : t ≥ 0} is a zero–mean Gaussian process with covariance
function given as

Z
cov(V (t), V (s)) =

min (t,s)

Φ(t, u) · G(x(u)) · Φ(s, u)T du.

(12)

0

The functions Φ(t, s) and G(x) are as defined previously.
The proof of this theorem can be seen in Britton and Andersson [9] or in Ethier
and Kurtz [14].
Note that the covarian function in equation (12) is time-dependent. It implies that
the variance function is also time-dependent. Using equation (12), the variance function
of {V (t) : t ≥ 0} is
Σ(t) = cov(V (t), V (t))
Z t
=
Φ(t, u) · G(x(u)) · Φ(t, u)T du
0

6

(13)

where Φ(t, s) (0 ≤ s ≤ t) must satisfy equation (11). By differentiating equation (13)
with respect to t,

Z t
dΣ(t)
∂
T
= Φ(t, t) · G(x(t)) · Φ(t, t) +
Φ(t, u) · G(x(u)) · Φ(t, u)T du
dt
∂t 0
Z t
= G(x(t)) + ∂F (x(t))
Φ(t, s) · G(x(u)) · Φ(t, u)T du+
0
Z t
Φ(t, u) · G(x(u)) · Φ(t, s)T du · [∂F (x(t))]T
0

= G(x(t)) + ∂F (x(t)) · Σ(t) + Σ(t) · [∂F (x(t))]T ,

(14)

where Σ(0) = 0.
d Sigma(t)
Note that in the epidemic case, when we let t → ∞,
tends to zero
dt
because if the epidemic survives and takes off, it will reach the endemic stage and the
infection process converges to a certain distribution. Otherwise, if the epidemic dies
out, we will see in the later section that the stable stage is reached in the disease-free
stage. Both scenarios are in equilibrium level. Therefore, once one of these two cases
dΣ(t)
is attained,
definitely goes to zero as t → ∞.
dt
Now let Σ̂ and x̂ be the variance-covariance matrix and the stationary infection
point when the epidemic reaches endemic stage. Hence, as t → ∞,
dΣ(t)
= G(x̂) + ∂F (x̂) · Σ̂ + Σ̂ · [∂F (x̂)]T
dt
−G(x̂) = ∂F (x̂) · Σ̂ + Σ̂ · [∂F (x̂)]T .

0=

(15)

Figure 2 illustrates the distribution of {n−1 Yn (t)|Yn (t) > 0, Yn (t) ∼ q̃}.
4.1. Expected Mean Time to Extinction by CLT
Suppose that R0 > 1 and {Yn (t) : t ≥ 0} is in quasi– equilibrium (meaning that Yn (t)
enters endemic level). Then, the process will definitely mimic the deterministic model
and with equilibrium point ŷ = 1 − R0−1 .
Now, recall a drift function
F (y), where we define y ∈ Z in such a way that βl :
P
Z → [0, ∞) , l ∈ ∆ ⊆ Z, l βl (y) < ∞. In SIS, ∆ = {1, −1} represents the possible
transitions of the process, i.e., 1 for infection to occur and −1 for removal to occur.
Therefore,
F (y) = β1 (y) + (−1) · β−1 (y)
= λy(t) · (1 − y(t)) − γy(t)

(16)

Furthermore,




∂F (y)
∂F (y) =
= λ(1 − 2y(t)) − γ
∂y

7

(17)

(a) (n, m, R0 ) = (1000, 2, 5)

(b) (n, m, R0 ) = (3000, 2, 5)

(c) (n, m, R0 ) = (1000, 2, 8)

(d) (n, m, R0 ) = (3000, 2, 8)

Figure 2. Empirical distribution of SIS endemic level with CLT approximation

and
G(y) = β1 (y) + β−1 (y)
= λy(t) · (1 − y(t)) + γy(t).

(18)

Thus since we assume that the process is in quasi–equilibrium, y(t) = ŷ for all
t ∈ {t ≥ 0 : Yn (t) | Yn (0) 6= 0 ∼ q̃}. Therefore,
∂F (ŷ) = −(λ − γ),
γ
G(ŷ) = 2 · (λ − γ).
λ

(19)
(20)

Since (19) is continuous, then according to Theorem 4.1, for any t ≥ 0,
√


n

Yn (t)
− (1 − R0−1 )
n



D

−−−−→ N (0, σ 2 (t))

(21)

as n → ∞ and σ 2 (t) must satisfy
dσ 2 (t)
γ
= 2 · (λ − γ) − 2 · (λ − γ)σ 2 (t)
dt
λ

(22)

with initial value σ 2 (0) = 0. Therefore, the solution of equation (22), for t → ∞, is
σ 2 (t) =


γ
γ
1 − e−2(λ−γ)t → = R0−1 .
λ
λ

8

(23)

The eq.(21) is equivalent to
√


n


Yn (t)
−1
− (1 − R0 ) ≈ N (0, R0−1 )
n


Yn (t)
1
− (1 − R0−1 ) ≈ N 0,
n
nR0


n
−1
Yn (t) ≈ N (1 − R0 )n,
R0

(24)

for large n and large t.
Suppose that µn = (1 − R0−1 )n and σn2 = Rn for simpler notations. Let fn and
0
Fn denote the pdf and cdf of N (µn , σn2 ). Then using a continuity correction, we can
approximate quasi-stationary distribution in (10) as follows
fn (1)

 
1
1
Fn n +
− Fn
2
2
2
1
− (x−µn )
p
e 2σn2
2πσn2
= Z n+ 1
Z 1
2
2
2
2
1
1
− (x−µn )
− (x−µn )
p
p
e 2σn2 dx −
e 2σn2 dx
2πσn2
2πσn2
−∞
−∞


1−µn
1
σn φ
σn
1

=  n+ 1 −µ 
−µn
n
2
2
Φ
−
Φ
σn
σn

q̃1 ≈



(25)

where φ(x) and Φ(x) are pdf and cdf of standardised normal distribution. Therefore,
the approximate mean time to extinction of the SIS epidemic model by applying
Theorem 4.1 is
E(TQ ) =

1
γ q̃1




≈



Φ

γ
φ
σn
n+

1
2



− µn

σn

1 − µn
σ
! n
−Φ

−1


1
2



!

− µn 
σn

(26)

for large n and large t.

5. Conclusion
We have successfully determined the approximation of expected mean time to extinction of SIS model using quasi–stationary distribution . An interesting result is that
the empirical distribution of the SIS in endemic–stage followed normal distribution
as the population size went to infinity. By using the CLT of Ethier and Kurtz, it
turned out that the mean parameters of the SIS is the equilibrium points derived from
9

solving deterministic’s ODE. This result was not surprising since we have showed in
previous work that the stochastic models moved randomly around the deterministic
model’s equilibrium point. Furthermore, the variances of the models converged to some
positive constants since intuitively, as the population size went sufficiently large, the
stochastic paths moved randomly around its endemic–equilibrium point.

6. Appendix A

(
Ei
First, we need to show that for i, j ∈ T , E i E j =
0

if i 6= j
.
if i =
6 j

Proof. Note that since ui is orthonormal, then



E i E j = W −1/2 ui uTi W 1/2 W −1/2 uj uTj W 1/2




= W −1/2 ui uTi uj uTj W 1/2
(
Ei
if i = j
=
0
if i 6= j
(
1
since uTi uj =
0

(27)

if i = j
.
if i =
6 j

And secondly, we need to show that

n
X

E i = In×n .

i=1

Proof. Note that,
n
X







E i = W −1/2 u1 uT1 W 1/2 + . . . + W −1/2 un uTn W 1/2

i=1


= W −1/2 u1 uT1 + . . . + un uTn W 1/2
= W −1/2 In×n W 1/2 = In×n .

Note that exp (Qt) =

∞
X
Qk tk
k=0

exp (λi tE i ) =

k!

. Then, for any i ∈ T ,

∞
X
(λi t)k E k
i

k=0

(28)

k!

= Ei

∞
X
(λi t)k
k=0

10

k!

= exp (λi t)E i .

7. Appendix B




Notice E i = W −1/2 ui uTi W 1/2 . Then,
QT T =

n
X




λi W −1/2 ui uTi W 1/2

i=1





uTi W 1/2 QT T =

n
X





λi uTi W 1/2 W −1/2 ui uTi W 1/2
{z
}
|
i=1
=1

=

n
X





λi uTi W 1/2 .

i=1



Hence, uTi W 1/2 is the left-eigenvector of QT T .

8. References
References
[1] H. Andersson and T. Djehiche, “Stochastic Epidemics in Dynamic Populations: QuasiStationary and Extinction”, J. Math. Bio. 41, 559–580 (2000).
[2] H. Andersson and B. Djehiche, “A Threshold Limit Theorem for the Stochastic Logistic
Epidemic”, J. Appl. Prob. 35, 662–670 (1998).
[3] N. T. Bailey, The Mathematical Theory of Infectious Diseases and its Applications, 2nd
edition (Griffin, London, 1975).
[4] F. G. Ball, “The Threshold Behaviour of Epidemic Models”, J. Appl. Prob. 20, 227–241
(1983).
[5] F. G. Ball, Epidemic Thresholds. Encyclopedia of Biostatistics. 3. (John Wiley & Sons,
Chichester, 1998).
[6] F. G. Ball, T. Britton, and P. Neal, “On Expected Durations of Birth-Death Processes,
with Application to Branching Processes and SIS Epidemics”, Submitted J. Appl. Prob.
53, 203–215 (2016).
[7] F. G. Ball and O. D. Lyne, “Optimal Vaccination Policies for Stochastic Epidemics among
a Population of Households”, Submitted Math. Bioscie. 177–178, 333–354 (2002).
[8] P. G. Ballard, N. G. Bean, J. V. Ross, “The Probability of Epidemic Fade-Out is NonMonotonic in Transmission Rate for Markovian SIR Model with Demography”, Submitted
J. Theoretical Bio. 393, 170–178 (2016).
[9] H. Andersson and T. Britton, Stochastic Epidemic Models and Their Statistical Analysis
(Springer, New York, 2000).
[10] T. Britton and P. Neal “The Time to Extinction for an SIS-houshold-epidemic Model”,
Submitted J. Math. Bio. 61, 763–779 (2010).
[11] D. Clancy and P. K. Pollett, “A Note on Quasi-Stationary Distributions of Birth-Death
Processes and the SIS Logistic Epidemic”, Submitted J. Appl. Prob. 40, 821–825 (2003).
[12] J.N. Darroch and E. Seneta, “On Quasy-Stationary Distributions in Absorbing DiscreteTime Finite Markov Chains”, Submitted J. Appl. Prob. 2, 88–100 (1965).
[13] J.N. Darroch and E. Seneta, “On Quasy-Stationary Distributions in Absorbing
Continuous-Time Finite Markov Chains”, Submitted J. Appl. Prob. 4, 192–196 (1967).
[14] S. N. Ethier and T. G. Kurtz, Markov Processes Characterization and Convergence (Wiley,
New York, 1982).

11

[15] s.‘Karlin, A First Course in Stochastic Processes. 2nd ed. (Academic Press Inc., California, 1975).
[16] F. P. Kelly, Reversibility and Stochastic Networks (Cambridge University Press, Cambridge, 2011).
[17] D. G. Kendall, “On the Generalized Birth and Death Process”, The Annals of Mathematical Statistics 19, 1–15 (1948).
[18] R. J. Kryscio and C. Lefévre, “On the Extinction of SIS Stochastic Logistic Epidemic”,
J. Appl. Prob. 26, 685–694 (1989).
[19] I. Nåsell, “Stochastic Models of Some Endemic Infections”, Math. Bioscie. 179, 1–19
(2002).
[20] R.H. Norden, “On the Distribution of The Time to Extinction in the Stochastic Logistic
Population Model”, Adv. in Appl. Prob. 14, 687–708 (1982).
[21] S.M. Ross, Introduction to Probability Models. 9th ed. (Elsevier Inc., London, 2007).
[22] K. Susvitasari, “Comparing the Behaviour of Deterministic and Stochastic Model of SIS
Epidemic”, arXiv:1803.01497 [q-bio.PE]

12

