This article appeared in a journal published by Elsevier. The attached
copy is furnished to the author for internal non-commercial research
and education use, including for instruction at the authors institution
and sharing with colleagues.
Other uses, including reproduction and distribution, or selling or
licensing copies, or posting to personal, institutional or third party
websites are prohibited.
In most cases authors are permitted to post their version of the
article (e.g. in Word or Tex form) to their personal website or
institutional repository. Authors requiring further information
regarding Elsevier’s archiving and manuscript policies are
encouraged to visit:
http://www.elsevier.com/copyright

Author's personal copy
Moment closure approximations for stochastic kinetic models with rational
rate laws
Peter Milner, Colin S. Gillespie ⇑, Darren J. Wilkinson
School of Mathematics & Statistics, Newcastle University, Newcastle upon Tyne, NE1 7RU, United Kingdom
a r t i c l e i n f o
Article history:
Received 3 September 2010
Received in revised form 27 January 2011
Accepted 11 February 2011
Available online 19 February 2011
Keywords:
Moment closure
Rational
Stochastic
Hill
Michaelis–Menten
Holling
a b s t r a c t
Stochastic models are often used when modelling chemical species that have low numbers of molecules.
However, as these models become large, it can become computationally expensive to simulate even a sin-
gle realisation of the system since even efficient simulation techniques have a high computational cost.
One possible technique to approximate the stochastic system is moment closure. The moment closure
approximation is used to provide analytic approximations to non-linear stochastic models. Until now,
this approximation has only been applied to models with polynomial rate laws. In this paper we extend
the moment closure method to cover models with rational rate laws.
 2011 Elsevier Inc. All rights reserved.
1. Introduction
Stochastic population models are becoming increasingly useful
in a variety of fields see for example [1–4]. For most stochastic
models of interest the transition probabilities are non-linear,
which makes the resultant stochastic process intractable.
To gain insight into the stochastic process we simulate from the
model. Exact realisations from the stochastic model can be ob-
tained using the discrete event simulation method described by
Kendall and later by Gillespie [5,6]. However, for large models this
leads to computational problems even for simulating a single real-
isation, let alone the many repeated realisations that are needed in
order to understand the stochastic variation inherent in the
process.
Other simulation methods have been proposed. The next-reac-
tion method [7] is a modification of a variant of the first reaction
method [8]. This algorithm can increase the speed of the simula-
tion. This increase in speed is achieved by minimising the recalcu-
lation of hazards and reducing the random numbers needed per
step. However this increase depends on the model and the random
number generator used. Approximate methods are also possible.
The s-leap method was developed by Gillespie [9] as an approxi-
mate method to accelerate the simulation process with an accept-
able loss of accuracy, by moving ahead by a variable amount of
time s and approximating the number of reactions in the interval
s with a Poisson random variable. For each leap, s is selected to
give a trade-off between speed and accuracy, where an increase
in speed is given by a large s and the accuracy is measured by
the difference between the reaction hazards at each end of the
interval, as this affects the assumption of constant hazard over
the interval s. The original s-leap algorithm has been improved
over the years [10,11].
In contrast, moment closure methods derive moment or
cumulant ODEs which can be solved numerically to obtain esti-
mates of the moments at a particular time [12–18]. Typically for
non-linear models, the ODEs for the moments form an infinite
set, i.e. the ith moment depends on at least the (i + 1)th moment.
This results in a set of ODEs that can not be solved analytically
or numerically. Moment closure is a technique to ‘close’ this set,
by eliminating the dependence on higher order moments, allowing
us to get a solvable set of coupled ODEs. This is usually done by set-
ting moments or cumulants above a certain order to zero, and solv-
ing the remaining coupled ODEs. For instance [19,20] set all third
and higher order cumulants to zero, therefore obtaining a normal
approximation.
Moment closure yields a large speed up in computational time
for the simulation of systems with large populations. For example,
in calcium-induced calcium release in cardiac myocytes [21], an
algebraic relationship (based on the beta-binomial distribution)
is used to express the third moment in terms of lower moments.
However all previous research applies moment closure to mod-
els where the rate laws are of polynomial form. This is the first
0025-5564/$ - see front matter  2011 Elsevier Inc. All rights reserved.
doi:10.1016/j.mbs.2011.02.006
⇑ Corresponding author.
E-mail address: c.gillespie@ncl.ac.uk (C.S. Gillespie).
Mathematical Biosciences 231 (2011) 99–104
Contents lists available at ScienceDirect
Mathematical Biosciences
journal homepage: www.elsevier .com/locate /mbs

Author's personal copy
paper that generates and analyses systems using rational rate laws.
In this paper we consider three specific models, which can be used
as building blocks in other systems. This extension to the moment
closure method now allows us to apply moment closure tech-
niques to models with more complex propensity functions.
2. Example 1: A simple rational rate law
To illustrate our method we will first look at the simplest ra-
tional rate law. In this model, we have a single species and a single
reaction
X ! ; ð1Þ
with propensity function
aðxÞ ¼ Vx
k þ x ; ð2Þ
where k > 0, V > 0 and x is the number of individuals of X at time t.
Thus the probability of (1) occurring in the interval (t, t + dt) is
a(x)dt + o(dt).
Holling [22] proposed (1) with propensity function a(x) as the
Holling type II response function. Holling used this function form
to describe the uptake of substrate. It should be noted that this
reaction system also describes the Michaelis–Menten system see
[23].
The associated forward Kolmogorov equation for (1) with pro-
pensity function a(x) is
dPðx; tÞ
dt
¼ Pðx þ 1; tÞ Vðxþ 1Þ
kþ x þ 1 Pðx; tÞ
Vx
kþ x ; ð3Þ
where P(x; t) is the probability of observing x individuals at time t
see [1]. To find the associated moment equations of (1) we first mul-
tiply (3) by (k + x + 1)(k + x) to give
dPðx; tÞ
dt
ðkþ xþ 1Þðk þ xÞ ¼ Pðx þ 1; tÞVðx þ 1Þðkþ xÞ
 Pðx; tÞVsðkþ x þ 1Þ: ð4Þ
The univariate moment generating function is defined as
Mðh; tÞ ¼
X1
x¼0
Pðx; tÞexh ¼
X1
i¼0
liðtÞh
i
i!
;
where
liðtÞ ¼ E½XðtÞ
i ¼
X1
i¼0
Pðx; tÞxðtÞi:
Multiplying Eq. (4) by exh and summing over all possible values of x
gives
X
i
hi
@ iþ1Mðh; tÞ
@t@hi
¼
X
i
fi
@iMðh; tÞ
@hi
eh 
X
i
gi
@iMðh; tÞ
@hi
; ð5Þ
where
hi ¼
k2 þ k for i ¼ 0;
2kþ 1 for i ¼ 1;
1 for i ¼ 2;
0 otherwise;
8
>>><
>>:
f i ¼
Vðk  1Þ for i ¼ 1;
V for i ¼ 2;
0 otherwise
8
><
>:
and
gi ¼
Vðkþ 1Þ for i ¼ 1;
V for i ¼ 2;
0 otherwise:
8
><
>:
By equating coefficients of h in Eq. (5) we get expressions for each of
the moments. For example, extracting coefficients of h0 gives an
expression for the first moment
ð2kþ 1Þdl1
dt
¼ 2Vl1 
dl2
dt
: ð6Þ
By rewriting Eq. (6) in terms of cumulants we obtain
dj1
dt
¼  Vj1
kþ 1=2þ j1
 dj2
dt
; ð7Þ
where j1 is the mean and j2 is the variance. Setting j2 = 0 in Eq. (7)
gives an approximation to the deterministic Eq. (2). When we have
polynomial rate laws, the deterministic equation is identical to the
moment closure equation with j2 set equal to zero. However, for ra-
tional rate laws this equality does not hold. Instead, Eq. (7) has an
additional term ‘1/2’, in the denominator.
If we equate coefficients of h1 and h2 we retrieve expressions for
the second and third moments respectively,
ð2kþ 1Þdl2
dt
¼ Vðk 1Þl1  3Vl2  kðkþ 1Þ
dl1
dt
 dl3
dt
ð2kþ 1Þdl3
dt
¼ Vðk 1Þl1 þ Vð3 2kÞl2  4Vl3
 kðk þ 1Þ dl2
dt
 dl4
dt
;
where li  li(t).
The ODE for l1(t) (see Eq. (6)) is not closed as it depends l2(t),
similarly the ODE for l2(t) depends on l3(t) etc. Whenever a model
includes non-linear rate laws this dependence structure appears.
Table 1 shows three distributions that could be used to close the
moment equations. Throughout this paper we will use the normal
approximation. Fig. 1 shows the mean with an approximate 95%
confidence interval from the moment closure approximation and
from 100000 Monte Carlo simulations. We can see that the mo-
ment closure approximation gives a good estimate of the true
mean and variance, although between t = 35 and t = 45 the approx-
imation deviates slightly from the true mean and variance.
3. General form
In this section we develop general multivariate results, for for-
ward Kolmogorov equations that contain rational rate coefficients,
allowing the calculation of moment equations up to any order.
Suppose we have N species {y1, . . . ,yN} and l reactions
{R1, . . . ,RL}, where reaction Rl corresponds to
sl1y1 þ    þ slNyN ! sl1y1 þ    þ slNyN;
where sl and sl are the number of reactants and products of each
species involved with reaction l.
Let y be a column vector of species numbers, sl ¼ sl  sl and
sli ¼ sli  sli be the stoichiometric coefficient of species yi in reaction
Rl. Then when reaction Rl occurs yi? yi + sli. Thus we can write the
propensity function for reaction Rl as, al(y)/bl(y), where a and b are
polynomials.
Let p(y)(t) = p(y) be the probability of being in state y at time t,
with initial conditions of y(0). The time evolution of y can be for-
mulated as the forward Kolmogorov equation
Table 1
Table showing some options for closing second order moment for different
distributional assumptions.
Distribution Moments Cumulants
Normal l3 ¼ 3l2l1  2l31 j3 = 0
Poisson l3 ¼ l1 þ 3l21 þ l31 j3 = j2 = j1
Log-normal l3 = (l2/l1)3 j3 ¼ ðj2=j1Þ3 þ 3j22=j1
100 P. Milner et al. /Mathematical Biosciences 231 (2011) 99–104

Author's personal copy
d
dt
PðyÞ ¼
XL
l¼1
Pðy  slÞ
alðy  slÞ
blðy  slÞ
 PðyÞalðyÞ
blðyÞ
; ð8Þ
where bl()– 0 and bl(0,0, . . . ,0)– 0. If we take bl()  1 in (8) then
we have the standard form for the CME. Multiplying (8) by
YL
l¼1
blðy  slÞblðyÞ;
yields
d
dt
PðyÞ
YL
l¼1
blðy  slÞblðyÞ
¼
XL
l¼1
Pðy  slÞalðy  slÞ
YL
r–l
brðy  srÞ
!"

YL
s¼1
bsðyÞ
!#

XL
l¼1
PðyÞalðyÞ
YL
r–l
brðyÞ
! YL
s¼1
bsðy  slÞ
!" #
:
On extracting coefficients of yi and (y  sl)i for each reaction l, yields
d
dt
PðyÞhðyÞ ¼
XL
l¼1
½Pðy  slÞf ðy  slÞ  PðyÞgðyÞ: ð9Þ
where
hðyÞ ¼
YL
l¼1
blðy  slÞbðyÞ ¼
X
i
hiyi; ð10Þ
f ðy  slÞ ¼
XL
l¼1
alðy  slÞ
YL
r–l
brðy  srÞ
YL
s¼1
bsðyÞ
¼
XL
l¼1
X
i
fl;iðy  slÞi; ð11Þ
gðyÞ ¼
XL
l¼1
alðyÞ
YL
r–l
brðyÞ
YL
s¼1
bsðy  slÞ ¼
X
i
giy
i ð12Þ
with hi and gi being the respective coefficients of yi ¼ yi11      y
iN
N ,
and fli being the coefficients of ðy  slÞi ¼ ðy1  sl1Þ
i1
     ðyN  slNÞ
iN for each reaction l. On multiplying (9) by eyh
and summing over all values of y, we get a PDE for the m.g.f, viz.
X
i
hi
@iþ1Mðh; tÞ
@t@hi
¼
XL
l¼1
X
i
fl;i
@iMðh; tÞ
@hi
expðslhÞ
" #

X
i
gi
@iMðh; tÞ
@hi
: ð13Þ
Since h and g are coefficients of P(y), they are independent of l. On
expanding the partial derivatives in (13), we get
@
@t
X
i
hi
X1
j¼i
hji
ðj iÞ!lj ¼
XL
l¼1
X
i
fl;i
X1
j¼i
hji
ðj iÞ!lj
X1
k¼0
ðslhÞk
k!
" #

X
i
gi
X1
j¼i
hji
ðj iÞ!lj; ð14Þ
where
ðslhÞk
k!
¼ ðsl1h1Þ
k1
k1!
     ðslNhNÞ
kN
kN!
and
lj ¼ lj1 ;j2 ;...;jN ðtÞ:
Extracting the coefficients of h from (14) yields
X1
n¼0
hn
X
i
hi
@
@t
liþn ¼
X1
n¼0
hn
XL
l¼1
X
i
fl;i
Xn
k¼0
skl
n
k
 
lnkþi
" #

X
i
gilnþi
( )
;
ð15Þ
where
skl
n
k
 
¼ sk1l1
n1
k1
 
     skNlN
nN
kN
 
and
lnkþi ¼ ln1k1þi1 ;...;nNkNþiN ðtÞ:
To obtain an equation for a particular moment, we extract the coef-
ficients of h1, . . . ,hN from each side of Eq. (15). Hence Eq. (15) is a
general result that can be applied to a wide class of stochastic mod-
els. For example, to obtain the moment equation for the marginal
mean, we extract the coefficients of hj to get an equation for the first
moment of species j. If we required an equation for l1,0,. . .,0 = lj, we
extract coefficients of h1 only, yielding
Time
A
m
ou
nt
0
50
100
150
200
0 10 20 30 40 50
Time
A
bs
ol
ut
e
D
iff
er
en
ce
−6
−4
−2
0
2
4
6
0 10 20 30 40 50
Fig. 1. (a) The mean ± 2 standard deviations obtained from 100,000 simulations (---) and the approximate mean ± 2 standard deviations from the moment closure
approximation (solid), where x(0) = 180 and {V,k} = {4.8,5/3}. (b) The error (approximation-exact) for the mean (solid), the upper 95% (---) and the lower 95% (---).
P. Milner et al. /Mathematical Biosciences 231 (2011) 99–104 101

Author's personal copy
X
i
hi
@
@t
liþn ¼
XL
l¼1
X
i
fl;isnl lnkþi
!

X
i
gilnþi; ð16Þ
where n = 1,0, . . . ,0. On extracting @lj/@t we gain the following
expression
hj
@lj
@t
¼
XL
l¼1
X
i
fl;ili 
X
i
gili 
X
i–j
hi
@li
@t
; ð17Þ
where lj = l1,0,. . .,0 and hj = h1,0,. . .,0 – 0.
3.1. Single species case
If there is only a single species, Eq. (15) simplifies to
X1
n¼0
hn1
X1
i¼0
hi
@
@t
liþn ¼
X1
n¼0
hn1
XL
l¼1
X1
i¼0
fl;i
Xn
k¼0
skl
n
k
 
lnkþi
" #(

X1
i¼0
gilnþi
)
: ð18Þ
This allows us to gain a simple approximation for the first moment,
which cannot be done for multiple species. In the single species case
to find the equation for the first moment we extract the coefficients
of h01, hence the equation for the first moment is
h1
@l1
@t
¼
XL
l¼1
X1
i¼0
fl;ili
" #

X1
i¼0
gili 
X
i–1
hi
@li
@t
ð19Þ
and the equation for the nth moment (nP 2) is
h1
@ln
@t
¼
XL
l¼1
X1
i¼0
fl;i
Xn
k¼0
skl
n
k
 
lnkþi
" #

X1
i¼0
gilnþi 
X1
i–1
hi
@ln1þi
@t
:
ð20Þ
This simplification is not possible for models with more than one
species as this approximation to the first moments from the coeffi-
cients of h01, gives one equation with N (number of species)
unknowns.
4. Further examples
In this section we will apply our method to two further exam-
ples which have rational rate laws. To test the speed and accuracy
of our moment closure approximations, we have simulated ‘exact’
means and variances for each of the examples. The means and vari-
ances are calculated from 100000 simulations of a Gillespie algo-
rithm coded in C. Each of the moment closure approximations
was also coded in C. The moment equations were solved using
the Runge Kutta Cash Karp (rkck) method from the GSL ODE li-
braries. All examples are run on a Intel Xeon E5345 2.33 GHz
processor.
4.1. Example 2: Rate laws involving powers
In this example we consider rational rate laws involving pow-
ers. This model contains a single species Y and a single reaction:
Y ! ;
with rate law
Vy2
kþ y2 :
This is a fairly standard reaction rate and is known as Holling Type
III response function for modelling vertebral predators. The rate has
in a variety of scenarios, such as modelling the outbreak of bud-
worms in forests and predator–prey systems see [22,24,25].
The forward Kolmogorov equation for this reaction is
dPðyÞ
dt
¼ Pðy slÞ
aðy slÞ
bðy slÞ
 PðyÞaðyÞ
bðyÞ ; ð21Þ
where a(y) = Vy2, b(y) = k + y2 and sl = 1. From relationship (19),
we obtain an equation for the first moment is
dl1ðtÞ
dt
¼ 1
h1
X4
i¼2
filiðtÞ 
X4
i¼2
giliðtÞ 
X4
i¼2
hi
dliðtÞ
dt
!
¼ 1
2k
4Vl3ðtÞ  ð2kþ 1Þ
dl2ðtÞ
dt
 2dl3ðtÞ
dt
 dl4ðtÞ
dt
 
:
ð22Þ
Similarly, from relationship (20) we obtain an equation for the sec-
ond moment
2k
dl2ðtÞ
dt
¼ Vððkþ 1Þl2ðtÞ  2l3ðtÞ þ 5l4ðtÞÞ  ðk
2 þ kÞdl1ðtÞ
dt
 ð2kþ 1Þdl3ðtÞ
dt
 2dl4ðtÞ
dt
 dl5ðtÞ
dt
:
This slightly more complicated rate law leads to the ODE for the
mean (Eq. (22)) depending on l2(t), l3(t) and l4(t), whereas in
Example 1 (Section 2), the ODE for l1(t) (Eq. (6)) only depends on
l2(t). We can again apply a normal approximation to close the set
of ODEs, allowing us to solve numerically for given initial conditions
and rate constants.
Fig. 2 shows a plot of the exact mean and variance against the
moment closure approximations of the mean and variance. We
can see that the moment closure approximation gives a good esti-
mate of the true values. As we saw in Example 1 the moment clo-
sure approximation is furthest from the true values as the
population of species Y approaches zero.
4.2. Example 3: Rate laws with multiple species
We now consider a Michaelis–Menten type system where the
product (P) can reform the substrate enzyme complex (SE). We also
consider that the product can be removed from the system and the
substrate immigrates into the system. Thus we have the following
reactions,
R1 : ; ! S; R2 : P ! ;; R3 : S ! P; R4 : P ! S;
with rates k, kP, (V1S)/(Cs + S + CrP) and (CrV2P)/(Cs + S + CrP) respec-
tively. The associated forward Kolmogorov equation is
dPðyÞ
dt
¼
X4
l¼1
Pðy  slÞ
alðy  slÞ
blðy  slÞ
 PðyÞalðyÞ
blðyÞ
; ð23Þ
where y = (y1,y2) are the numbers of molecules of substrate and
product respectively. From Eqs. (10)–(12) we can calculate the coef-
ficients hi, fl,i and gi. Eq. (16) thus allows us to form expressions for
moments of our choice. We have once again used a normal approx-
imation to close the moment equations, yielding five coupled ODEs
(two means, two variances and covariance). These can be solved
numerically to give estimates of the means and variances. We will
consider three approximations to the means, variances and
covariance.
 Full model (Full). All terms for the means, variances and covari-
ance in the moment equations are included.
 Deterministic means (Approximation 1). By assuming the vari-
ance and covariance are zero in the moment equations for the
means, we can gain an approximation to the deterministic solu-
tion. We can combine this approximation with the equations for
the variances and covariance to get a first approximation.
 Assumeall differentials on the LHS equal zero (Approximation 2).
102 P. Milner et al. /Mathematical Biosciences 231 (2011) 99–104

Author's personal copy
We assume that the rate of change for each moment is small
compared to the moment itself, equivalently all hi = 0 " i– 1. That
is we only consider the differential on the LHS that we are inter-
ested in.
For each of these schemes we must solve five moment ODEs.
Fig. 3 shows the mean and variance estimates for each scheme,
along with a Monte Carlo estimate of the true values. We observe
that the full model and both approximations estimate match the
Monte Carlo estimates fairly closely. As we would expect, Approx-
imation 2 performs the poorest of the three schemes.
Table 2 shows the CPU time taken (in seconds) for each of the
three schemes for different simulation times. Overall, using
approximations {1,2} compared to the full moment closure ODES
gives a speed increase of about a factor of {10,100}, respectively.
The CPU time scales linearly with the increase in simulation time.
5. Conclusions
It is becoming increasingly recognised that to accurately cap-
ture the fundamental characteristics of many systems, the associ-
ated mathematical model should contain stochastic features. For
example, Proctor et al. [26] modelled DNA damage-response mech-
anisms in Saccharomyces cerevisiae. This model highlighted the
interaction between DNA repair and checkpoint pathways. Cellular
damage, by its very nature, is a stochastic event. Proctor et al. high-
lights that while an ODE representation of the model captures the
average behaviour, it fails to highlight the variability of the number
of cell divisions in the cdc13-1 mutant strains. Even though this
model was based on mass-action kinetics, it was necessary to make
some limiting assumptions leading to rational laws.
Time
A
m
ou
nt
0
50
100
150
200
0 10 20 30 40 50
Time
A
bs
ol
ut
e
D
iff
er
en
ce
−6
−4
−2
0
2
4
6
0 10 20 30 40 50
Fig. 2. (a) The mean ± 2 standard deviations for Example 2 obtained from 100000 simulations (---) and the approximate mean ± 2 standard deviations from the moment
closure approximation (solid), where y(0) = 180 and {V,k} = {4.8,5/3}. (b) The error (approximation-exact) for the mean (solid), the upper 95% (---) and the lower 95% (---).
Time
A
m
ou
nt
0
50
100
150
200
250
300
Approximation 1
0 2500 5000
Approximation 2
0 2500 5000
Full
0 2500 5000
Fig. 3. The exact mean ± from 100000 simulations 2 standard deviations (---) and the approximate mean ± 2 standard deviations from the moment closure approximations
(---) for Example 3. The initial conditions are y(0,0) = {300,0} and parameters {k,k,V1,V2,Cs,Cr} = {1,0.01,5,5,2,2}.
Table 2
Comparison of simulation times for different moment closure schemes, where
y(0,0) = {300,0} and {k,k,V1,V2,Cs,Cr} = {1,0.01,5,5,2,2}.
Simplification Simulation time CPU time (s)
Full 50 0.57
Approximation 1 50 0.088
Approximation 2 50 0.004
Full 5000 48.75
Approximation 1 5000 7.52
Approximation 2 5000 0.14
P. Milner et al. /Mathematical Biosciences 231 (2011) 99–104 103

Author's personal copy
Waldherr et al. [27] considered stochastic and deterministic
models which linked cellular decisions taking place on a time scale
of years to decades to cellular decisions that were measured in
minutes to hours. By using stochastic bistable switch models
(involving rational rate laws) they demonstrated that population
traits could be predicted from noisy stochastic systems. In their pa-
per they considered only basic models, so future work involving
stochastic models will, by necessity, involve approximations (pos-
sibly using moment closure) of the system.
Stochastic hybrid systems (SHS) are increasingly being used by
modellers see [28] for a review. SHS decompose the state space of
the model into a continuous part and a discrete part. Singh and
Hespanha highlight that moment closure approximations would
be ideal in this setup. The general idea is that fast reactions could
be simulated using a fairly accurate moment closure approxima-
tion, whilst the slow reaction be simulated separately. The key
point is that we are not ignoring the stochasticity of the fast reac-
tions. The method described in this paper provides another tool for
stochastic modellers, allowing moment closure techniques to be
applied to a wider class of models.
References
[1] D.J. Wilkinson, Stochastic Modelling for Systems Biology, Chapman and Hall
CRC Press, 2006.
[2] V. Isham, Assessing the variability of stochastic epidemics, Math. Biosci. 107
(1991) 209.
[3] E. Renshaw, Modelling Biological Populations in Space and Time, Cambridge
University, 1991.
[4] B. Bolker, S.W. Pacala, Using moment equations to understand stochastically
driven spatial pattern formation in ecological systems, Theor. Popul. Biol. 52
(1997) 179.
[5] D.G. Kendall, An artificial realisation of a simple birth and death process, J. R.
Stat. Soc. B 12 (1950) 116.
[6] D.T. Gillespie, Exact stochastic simulation of coupled chemical reactions, J.
Chem. Phys. 81 (1977) 2340.
[7] M.A. Gibson, J. Bruck, Efficient exact stochastic simulation of chemical systems
with many species and many channels, J. Phys. Chem. 104 (2000) 1876.
[8] D.T. Gillespie, A general method for numerically simulating the stochastic time
evolution of coupled chemical reactions, J. Comput. Phys. 22 (1976) 403.
[9] D.T. Gillespie, Approximate accelerated stochastic simulation of chemically
reacting systems, J. Chem. Phys. 115 (2001) 1716.
[10] Y. Cao, D.T. Gillespie, L.R. Petzold, Avoiding negative populations in explicit
poisson tau-leaping, J. Chem. Phys. 123 (2005) 054104.
[11] Y. Cao, D.T. Gillespie, L.R. Petzold, Efficient step size selection for the tau-
leaping simulation method, J. Chem. Phys. 124 (2006) 044109.
[12] A. Gandhi, S. Levin, S. Orszag, Moment closure in spatial ecological models and
moment closure through gaussian approximation, Bull. Math. Biol. 62 (2000)
595.
[13] J.H. Matis, T.R. Kiffe, Effects of immigration on some stochastic logistic models:
a cumulant truncation analysis, Theor. Popul. Biol. 56 (1999) 139.
[14] M.J. Keeling, Multiplicative moments and measures of persistence in ecology, J.
Theor. Biol. 205 (2000) 269.
[15] G. Marion, E. Renshaw, G. Gibson, Stochastic effects in a model of nematode
infection in ruminants, IMA J. Math. Appl. Med. Biol. 15 (1998) 97.
[16] I. Krishnarajah, A. Cook, G. Marion, G. Gibson, Novel moment closure
approximations in stochastic epidemics, Bull. Math. Biol. 67 (2005) 855.
[17] I. Krishnarajah, G. Marion, G. Gibson, Novel bivariate moment-closure
approximations, Math. Biosci. 208 (2007) 621.
[18] C.S. Gillespie, Moment closure approximations for mass-action models, IET
Syst. Biol. 3 (2009) 52.
[19] L.A. Goodman, Population growth of the sexes, Biometrics 9 (1953) 212.
[20] P. Whittle, On the use of the normal approximation in the treatment of
stochastic processes, J. R. Stat. Soc. Ser. B (Methodological) 19 (1957) 268.
[21] G.S. Williams, M.A. Huertas, E.A. Sobie, M.S. Jafri, G.D. Smith, Moment closure
for local control models of calcium-induced calcium release in cardiac
myocytes, Biophys. J. 95 (2008) 1689.
[22] C.S. Holling, The functional response of invertebrate predators to prey density,
Mem. Entomol. Soc. Can. 48 (1966) 1.
[23] L. Michaelis, M.L. Menten, Die kinetik der invertinwirkung, J. Biochem. 49
(1913) 333.
[24] D. Ludwig, D.D. Jones, C.S. Holling, Qualitative analysis of insect out-break
systems: the spruce budworm and forest, J. Anim. Ecol. 47 (1978).
[25] S.B. Hsu, T.W. Huang, Global stability for a class of predator–prey systems,
SIAM J. Appl. Math. 55 (1995) 763.
[26] C.J. Proctor, D.A. Lydall, R.J. Boys, C.S. Gillespie, D.P. Shanley, D.J. Wilkinson,
T.B.L. Kirkwood, Modelling the checkpoint response to telomere uncapping in
budding yeast, J. R. Soc., Interface/R. Soc. 4 (2007) 3.
[27] S. Waldherr, J. Wu, F. Allgöwer, Bridging time scales in cellular decision
making with a stochastic bistable switch, BMC Syst. Biol. 4 (2010) 108.
[28] A. Singh, J.P. Hespanha, Stochastic hybrid systems for studying biochemical
processes, Philos. Trans. R. Soc. A 368 (2010).
104 P. Milner et al. /Mathematical Biosciences 231 (2011) 99–104