Bulletin of Mathematical Biology 67 (2005) 855 873
www.elsevier.com/locate/ybulm
Novel moment closure approximations in stochastic
epidemics
Isthrinayagy Krishnarajah
a,b,∗
,AlexCook
a,b
,GlennMarion
b
,
Gavin Gibson
a
a
Department of Actuarial Mathematics and Statistics, Heriot-Watt University, Edinburgh EH14 4AS,
United Kingdom
b
Biomathematics and Statistics Scotland, JCMB, The King s Buildings, Edinburgh EH9 3JZ, United Kingdom
Received 11 August 2003; accepted 11 November 2004
Abstract
Moment closure approximations are used to provide analytic approximations to non-linear
stochastic population models. They often provide insights into model behaviour and help validate
simulation results. However, existing closure schemes typically fail in situations where the population
distribution is highly skewed or extinctions occur. In this study we address these problems by
introducing novel second- and third-order moment closure approximations which we apply to the
stochastic SI and SIS epidemic models. In the case of the SI model, which has a highly skewed
distribution of infection, we develop a second-order approximation based on the beta-binomial
distribution. In addition, a closure approximation based on mixture distribution is developed in
order to capture the behaviour of the stochastic SIS model around the threshold between persistence
andextinction. This mixture approximation comprises a probability distribution designed to capture
the quasi-equilibrium probabilities of the system and a probability mass at 0 which represents the
probability of extinction. Two third-order versions of this mixture approximation are considered in
which the log-normal and the beta-binomial areused to model the quasi-equilibrium distribution.
Comparison with simulation results shows: (1) the beta-binomial approximation is exible in shape
and matches the skewness predicted by simulation as shown by the stochastic SI model and (2)
mixture approximations are able to predict transient and extinction behaviour as shown by the
∗
Corresponding author at: Department of Actuarial Mathematics and Statistics, Heriot-Watt University,
Edinburgh EH14 4AS, United Kingdom. Tel.: +44 131 650 7536; fax: +44 131 650 4901.
E-mail addresses: isthri@bioss.ac.uk (I. Krishnarajah), alex@bioss.ac.uk (A. Cook), glenn@bioss.ac.uk
(G. Marion), gavin@ma.hw.ac.uk (G. Gibson).
0092-8240/$30 ' 2004 Society for Mathem atical Biology. Published by Elsevier Ltd. All rights reserved.
doi:10.1016/j.bulm.2004.11.002

856 I. Krishnarajah et al. / Bulletin of Mathematical Biology 67 (2005) 855 873
stochastic SIS model, in marked contrast with existing approaches. We also apply our mixture
approximation to approximate a likelihood function and carry out point and interval parameter
estimation.
' 2004 Society for Mathematical Biology. Published by Elsevier Ltd. All rights reserved.
1. Introduction
Stochastic models are useful in epidemiology and ecology, and are used widely (Isham,
1991; Allenand Cormier, 1996; Bolker and Pacala, 1997; Filipe and Gibson, 1998; Marion
et al., 1998; Matis and Kiffe, 1999; Bauch and Rand, 2000; Keeling, 2000). Usually, the
transition probabilitiesexhibit non-linear dependence on population size or number of
infectives which makes the resultant stochastic processes analytically intractable. Hence,
techniques of approximation are needed to capture the underlying behaviour of the
stochastic processes. Linearisation is one such approximation, where the behaviour of
small stochastic uctuations can be examined around a xed point of the deterministic
dynamics (Bailey, 1963). An alternative approach is to analyse the quasi-equilibrium
probabilities which give a picture of the distribution independent of time and conditional on
extinction not having occurred (Renshaw, 1991). Both linearisation and quasi-equilibrium
probabilities are limited in their application to regions close to the xed points or at
equilibrium.
In contrast, closure methods are based on equations describing the temporal evolution of
moments or cumulants and in principle applytoboth transient and equilibrium dynamics.
One such widely used closure method is the cumulant truncation procedure (Matis and
Kiffe, 1996)wherethecumulant functions of say order k are approximated by setting
all cumulants of order higher than k to 0. Renshaw (1998) has shown that there is a
natural distribution associated with cumulant truncation. This saddlepoint approximation is
obtained by applying the method of steepest descents to the truncated cumulant generating
function and can be applied in multivariate situations (Renshaw, 2000). In our study, we
follow an alternative route using moment closure approximation based on distributional
assumptions, a technique introduced by Whittle (1957) which has been widely used
in recent years (Isham, 1991; Marion et al., 1998; Keeling, 2000; N sell, 2003). Most
commonly in these approximations, the population distribution is only described by the
rst- and second-order moments and descriptio nofextinction or bimodality is problematic.
Thus, we explore the use of moment closure using mixture approximations to population
distributions. Many existing methods of second-order approximation also have dif culties
in describing highly skewed distribution; hence we consider the application of a novel
second-order approximation based on the beta-binomial distribution.
Two generic epidemic models are studied: the stochastic SIS (susceptible infected
susceptible), as an example which exhibits extinction, and the stochastic SI
(susceptible infected), a special case of the SIS,usedasacase where the infected
population exhibits a highly skewed distribution but totality of infection is guaranteed.
Depending on the disease transmission rate, the SIS model exhibits meta-stable persistence
of disease, rapid extinction or a critical region corresponding to the border between the two