Bayesian inference for generalized stochastic popula-
tion growth models with application to aphids
Colin S. Gillespie and Andrew Golightly
Newcastle University, Newcastle Upon Tyne, UK.
Summary. In this paper we analyse the effects of various treatments on cotton aphids (Aphis
gossypii). The standard analysis of count data on cotton aphids determines parameter val-
ues by assuming a deterministic growth model and combines these with the corresponding
stochastic model to make predictions on population sizes, depending on treatment. Here, we
use an integrated stochastic model to capture the intrinsic stochasticity, of both observed aphid
counts and unobserved cumulative population size for all treatment combinations simultane-
ously. Unlike previous approaches, this allows us to explicitly explore and more accurately
assess treatment interactions. Markov chain Monte Carlo methods are used within a Bayesian
framework to integrate over uncertainty associated with the unobserved cumulative population
size and estimate parameters. We restrict attention to data on aphid counts in the Texas High
Plains obtained for three different levels of irrigation water, nitrogen fertiliser and block, but
note that the methods we develop can be applied to a wide range of problems in population
ecology.
Key Words: Cotton aphid, Markov jump process, Moment closure approximation,
Markov chain Monte Carlo.
1. Introduction
The importance of the cotton aphid as an insect pest is well documented (see for example
Leclant and Deguine, 1994). Consequently, accurate modelling of cotton aphid populations
is of economic importance (Rummel et al., 1995; Xia et al., 1999). In the recent litera-
ture, attempts to model cotton aphid dynamics have included the use of dispersion models
(Celini and Vaillant, 2004), computer simulation models (Giarola et al., 2006) and nonlinear
regression models (Matis et al., 2007).
Whilst deterministic population growth models for insect counts are reasonably well
developed (see for example Matis et al., 2006, 2008) little work has taken into account
stochastic effects that are particularly important when population sizes are small (Zheng
and Ross, 1991). Matis et al. (2007) propose a stochastic model for aphid growth, and
use this to make predictions on population sizes using parameter values obtained from the
analogous deterministic model. The model itself assumes a birth rate dependent on current
aphid numbers and a death rate dependent on the latter plus the unobserved cumulative
population size, giving a bivariate Markov jump process with an unobserved component.
We build on the work of Matis et al. (2007) by using the stochastic model to infer unknown
parameters (namely birth and death rates) using discrete time observations on aphid counts
only. The task of inferring parameters governing Markov jump processes with discrete state
space has been examined by Boys et al. (2008) amongst others. It is found therein that
working with the discrete stochastic model can be computationally prohibitive for mod-
els involving more than three reactions and two species. We therefore approximate the

2 Gillespie and Golightly
model and perform exact Bayesian inference for this approximate model. Such approxima-
tions have, in the past, focused on the use of coupled stochastic differential equations to
replace the underlying discrete valued process (see for example Golightly and Wilkinson,
2005; Heron et al., 2007; Purutcuoglu and Wit, 2007). Such an approach, however, is by no
means straightforward, and can be computationally prohibitive for large datasets (Golightly
and Wilkinson, 2008). Our approach also ignores discreteness but assumes that transition
densities for the number of aphids (and the cumulative population size) are Normal, with
parameters determined using moment closure techniques (Milner et al., 2009). Hence, given
complete information on both aphid counts and cumulative population size, Bayesian infer-
ence is straightforward for our approximate model. However, we do not observe cumulative
population size and therefore integrate over the uncertainty associated with it via a compu-
tationally efficient Markov chain Monte Carlo (MCMC) algorithm. In brief, we implement
an MCMC scheme to alternate between draws of the latent data (conditional on the current
parameter value and the observed data), and the parameter value (conditional on the ob-
served and latent data). The methods we develop are general and can in theory be applied
to any stochastic kinetic model.
As well as using synthetic data to validate the methodology, the algorithm is imple-
mented using the data in Matis et al. (2008) (see also Matis et al., 2007), that is, using
measurements on cotton aphid counts in the Texas High Plains. The data were obtained
for several controlled conditions including various irrigation water and nitrogen fertiliser
levels. We apply our method to estimate the effect of these treatments on the stochastic
birth and death rates. In addition, we integrate over the uncertainty associated with the
rate constants to give predictive distributions of the total number of aphids at future time
points. This is particularly effective for comparing the effects of different treatments on
aphid numbers and is therefore of practical importance.
The remainder of this paper is organised as follows: in Section 2, a description of the
model and the data is given. We also outline the moment closure approximation which will
form the basis of the inference algorithm described in Section 3. The algorithm is applied
to real and synthetic data in Section 4 before conclusions are drawn in Section 5.
2. The model
2.1. Description of the Model
Following Matis et al. (2006), we assume a linear aphid population birth rate of λN(t)
where N(t) denotes the size of the aphid population at time t. As discussed in Prajneshu
(1998), aphids excrete honey-dew, forming a cover on the leaf surface and causing aphid
starvation. As the area covered by excretion at time t is proportional to the cumulative
population size at time t, C(t), we assume a death rate of µN(t)C(t) and for simplicity,
assume that there is no removal or decomposition of honey-dew (Matis et al., 2006, 2007,
2008). The model can be represented by the coupled (pseudo) reactions,
N λ−−→ 2N + C (1)
N + C µ−−→ C (2)
since clearly, an occurrence of (1) will lead to a unit increase in both N and C whereas
the reaction in (2) will give a unit decrease in N but leave C unchanged. Such models
are known in the literature as stochastic kinetic models; see Wilkinson (2006) for a concise

Generalized stochastic growth models 3
introduction. We can express the probabilistic laws governing the time evolution of the
process by considering a small interval (t, t+ dt], namely
Pr {N(t+ dt) = n(t) + 1, C(t+ dt) = c(t) + 1 |n(t), c(t)} = λn(t)dt+ o(dt) (3)
Pr {N(t+ dt) = n(t)− 1, C(t+ dt) = c(t) |n(t), c(t)} = µn(t)c(t)dt+ o(dt), (4)
and we adopt the convention that upper case denotes random variable and lower case denotes
the realisation. The expressions in (3) and (4) describe a Markov jump process where each
event occurs at a particular rate dependent on the current state of the system. Note that
ignoring stochasticity leads to a deterministic model where the time course behaviour of
N(t) and C(t) is described by the set of differential equations,








dN(t)
dt = λN(t)− µN(t)C(t)
dC(t)
dt = λN(t)
Let pn,c(t) denote the probability that there are n aphids in the population at time t
and a cumulative population size of c. Then, in principle, inference may proceed by solving
the forward Kolmogorov equation to obtain pn,c(t). This is achieved by writing pn,c(t+∆t)
as the sum of the probabilities of arriving in state (n, c)′ at time t+∆t:
pn,c(t+∆t) = λ(n−1)pn−1,c−1(t)∆t+µc(n+1)pn+1,c(t)∆t+{1− n(λ+ µc)∆t} pn,c(t) . (5)
Plainly, the first two terms on the RHS of (5) give the probability that the system is one
event (either a birth or a death) removed from the state (n, c)′ at time t and then undergoes
such an event in (t, t + ∆t]. The last quantity in (5) is the probability that there are no
events in (t, t+∆t]. Hence, equation (5) leads to the forward Kolmogorov equation
dpn,c(t)
dt = λ(n− 1)pn−1,c−1(t) + µc(n+ 1)pn+1,c(t)− n(λ+ µc)pn,c(t), (6)
which is also known in the context of stochastic chemical kinetic models as the “Chemical
Master Equation”. If equation (6) could be solved to obtain pn,c(t) then Bayesian inference
would be straightforward. Even in the absence of a solution to (6), Boys et al. (2008)
show that it is possible to directly estimate the parameters in Markov jump processes using
Markov chain Monte Carlo algorithms. However, the methods proposed are not suitable
for any reaction system of reasonable size (say, more than three reactions with two species)
or when there is a large amount of data. The dataset we describe in Section 2.3 involves
twenty-seven treatment-block combinations, giving rise to a model with fifty-four reactions
and species. Therefore, it seems unlikely that we can carry out Bayesian inference based on
the exact underlying model. Instead, working with an approximate model but making exact
Bayesian inferences, appears to be a promising approach. Developments on approximating
the underlying model in the recent literature have focused on the diffusion approximation
(see for example Golightly and Wilkinson, 2005, 2006; Heron et al., 2007; Purutcuoglu and
Wit, 2007). The diffusion approximation of the stochastic kinetic model considered here is
obtained by calculating the infinitesimal mean and second moment of the process defined
by (3)–(4) and matching these to the drift and diffusion coefficient of an Itoˆ stochastic
differential equation. We obtain
{
dN(t) = (λN(t)− µN(t)C(t)) dt+
√
λN(t) dW1(t)−
√
µN(t)C(t) dW2(t)
dC(t) = λN(t)dt+
√
λN(t) dW1(t)
(7)

4 Gillespie and Golightly
where W1(t) and W2(t) are uncorrelated Brownian motions. Whilst it is possible to base
inference algorithms around the diffusion approximation, absence of an analytic solution
precludes closed form expressions for the associated transition densities. Consequently,
Bayesian inference is not straightforward. The approach of Eraker (2001) (see also Go-
lightly and Wilkinson, 2005) involves augmenting observations with latent data to allow
a sufficiently accurate Euler-Maruyama approximation of the underlying but unavailable
transition densities. Whilst this approach has been shown to work well, computational cost
scales almost linearly with the number of latent values employed. Therefore, we adopt a dif-
ferent approximation. As with the diffusion approximation, we also ignore the discreteness
associated with the Markov jump process but approximate transition densities as Normal,
with mean and variance determined by moment closure methods. Bayesian inference can
then be performed using these approximate transition densities. Since the moment clo-
sure approximation does not require the need for imputing latent values, the technique is
computationally efficient. The method is outlined in Section 2.2.
2.2. Moment Closure Approximation
For the forward Kolmogorov equation in (6), it is not possible to derive an analytic solution.
One common method of approximating the first few moments of the process is to examine
the moment equations of the system (see for example Krishnarajah et al., 2005; Singh and
Hespanha, 2007; Gillespie, 2009). We first define the bivariate moment generating function
as
M(θ, φ; t) ≡
∞
∑
n,c=0
enθ ecφpn,c(t)
and the associated cumulant generating function as
K(θ, φ; t) ≡ log[M(θ, φ; t)] =
∞
∑
n,c=0
θn
n!
φc
c! κnc .
When dealing with only the first few moments, cumulants are a particularly useful repre-
sentation since κ10 and κ01 are the marginal means of the n(t) and c(t) and {κ20, κ02, κ11}
are the marginal variances and covariance, respectively.
On multiplying equation (6) by enθecφ, summing over {n, c} and converting to the
cumulant generating function, we obtain the partial differential equation
∂K
∂t = λ(e
θ+φ − 1)∂K∂θ + µ(e
−θ − 1)
(
∂2K
∂θ∂φ +
∂K
∂θ
∂K
∂φ
)
(8)
after dropping the notational dependence of K(θ, φ; t) on θ, φ and t. A complete derivation
of equation (8) can be found in Appendix A. An ordinary differential equation for the
mean population can now be obtained by differentiating (8) with respect to θ and setting
θ = φ = 0. Likewise, to obtain an expression for the covariance we differentiate (8) with

Generalized stochastic growth models 5
respect to θ and φ and set θ = φ = 0. Thus we obtain the following cumulant equations
dκ10
dt = λκ10 − µ(κ10κ01 + κ11)
dκ01
dt = λκ10
dκ20
dt = λ(κ10 + 2κ20) + µ(κ11 − 2κ10κ11 − 2κ21 + κ01(κ10 − 2κ20))
dκ11
dt = λ(κ10 + κ20 + κ11)− µ(κ10κ02 + κ01κ11 + κ12)
dκ02
dt = λ(κ10 + 2κ11) . (9)
Since the forward Kolmogorov equation contains a non-linear rate µnc, the ith equation
depends on the (i + 1)th equation. To circumvent this dependency problem we assume
an underlying bivariate Normal distribution, that is we close higher-order moments. This
implies that
κ12(t) = κ21(t) = 0
in equations (9). Matis et al. (2006) extensively explores the validity of assuming an under-
lying Normal distributions for this model, and concludes that the approximation appears
to be suitable. In particular, they conclude
(a) The mean values κ10(t) and κ01(t) are well approximated and the expected final
cumulative count κ01(∞) has an error rate of less than 2.5%
(b) The approximation for the variance is also very accurate, with an error rate of less
than 2.5%.
(c) The skewness for the population and cumulative number of aphids can be substan-
tive, but the limiting distribution is approximated well by a Normal distribution. It
is worth noting that since we are using a truncated Normal distribution, then our
approximating distribution will have non-zero third cumulant.
Overall, the moment closure approximation provides a computationally efficient method for
describing aphid population and final cumulative count distributions (Matis et al., 2006).
2.3. Description of the Data
In this paper we consider the data described in Matis et al. (2008) consisting of five observa-
tions on cotton aphid counts on twenty randomly chosen leaves in each plot, for twenty-seven
treatment-block combinations. The data were recorded in July 2004 in Lamesa, Texas. The
treatments consisted of three nitrogen levels (blanket, variable and none), three irrigation
levels (low, medium and high) and three blocks, each being a distinct area. Irrigation treat-
ments were randomly assigned within each block as whole plots. Nitrogen treatments were
randomly assigned within each whole block as split plots. The data are plotted in Figure 1.
Note that the sampling times are t=0, 1.14, 2.29, 3.57 and 4.57 weeks (i.e. every 7 to 8
days).
Let i, j, k denote the respective level of water, nitrogen and block, with i, j, k ∈ {1, 2, 3}
where for simplicity 1 denotes low water/blanket nitrogen, 2 denotes medium water/variable
nitrogen and 3 denotes high water/zero nitrogen. Now denote the time series obtained for
combination ijk as Dijk, the number of aphids for combination ijk at time t by Nijk(t) and

6 Gillespie and Golightly
Time
Ap
hi
d
Po
pu
la
tio
n
0
500
1000
1500
2000
2500
0 1 2 3 4
Water (H)
Nitrogen (B)
Water (L)
Nitrogen (B)
0 1 2 3 4
Water (M)
Nitrogen (B)
Water (H)
Nitrogen (V)
Water (L)
Nitrogen (V)
0
500
1000
1500
2000
2500
Water (M)
Nitrogen (V)
0
500
1000
1500
2000
2500
Water (H)
Nitrogen (Z)
0 1 2 3 4
Water (L)
Nitrogen (Z)
Water (M)
Nitrogen (Z)
Fig. 1. Aphid counts against time obtained for 27 different scenarios – three levels of water (low,
medium and high) and nitrogen (blanket, variable and zero) and three blocks (with the data from
block 1 given by a circle, block 2 by a triangle and block 3 by a cross).
the cumulative population size at time t by Cijk(t). We then model each dataset Dijk with
a stochastic kinetic model of the form
Nijk
λijk−−−−→ 2Nijk + Cijk
Nijk + Cijk
µijk−−−−→ Cijk
and impose the following nested structure on the birth and death rates
λijk = λ+ αli + βlj + (αβ)lij + ρlk + (αρ)lik + (βρ)ljk (10)
µijk = µ+ αmi + βmj + (αβ)mij + ρmk + (αρ)mik + (βρ)mjk (11)
subject to the constraints
αl1 = αm1 = βl1 = βm1 = ρl1 = ρm1 = 0
and
(αβ)l11 = (αβ)m11 = 0, (αρ)l11 = (αρ)m11 = 0, (βρ)l11 = (βρ)m11 = 0 .
Thus, inference from these data amounts to making statements about thirty-eight parame-
ters. The interpretation of equations (10) and (11) is clear — all (5× 33 =)135 data points
are used to learn about the baseline birth and death rates λ and µ. The (5 × 32 =)45
observations recorded for the medium water scenario inform us about the effect of medium

Generalized stochastic growth models 7
water relative to baseline via αl2, αm2 and so on. Hence, we allow the overall birth and death
rates to vary with each treatment combination but assume underlying (baseline) birth and
death rates of λ and µ for all combinations. Naturally, the structure in equations (10) and
(11) allows us to explore possible interactions between treatments. The analysis of this
dataset is given in Section 4. In the next section we formulate the inference problem for
the stochastic growth model considered here and outline an MCMC approach.
3. Inference
For simplicity, consider first the task of inferring the birth and death rates λ and µ in the
model given by equations (1)–(4) from a single time series. Let X(tu) = (N(tu), C(tu))′ be
the random 2-vector of observed aphid counts and unobserved cumulative population size
at time tu. Let η = (λ, µ)′ be the 2-vector of parameters. Then, using the moment closure
approximation and assuming a Normal distribution, we have
X(tu) |X(tu−1),η ∼ N
(
ψu−1 , Σu−1
)
(12)
where ψu−1 andΣu−1 are calculated via the moment closure approximation outlined in Sec-
tion 2.2 and may depend on the parameters η and state X(tu−1) in a nonlinear way. Note
that N(ψ , Σ) denotes the multivariate Normal distribution with mean vector ψ and covari-
ance matrix Σ. Now suppose that observations on aphid counts n = {n(tu) : u = 0, . . . , T}
are available. Summarising our a priori beliefs about η and the initial unobserved cumula-
tive population size C(t0) via p(η) and p (c(t0)), we obtain the joint posterior distribution
for parameters and unobserved states as
p (η, c |n) ∝ p(η) p (c(t0))
T
∏
u=1
p (x(tu) |x(tu−1),η) (13)
where p (x(tu) |x(tu−1),η) is the transition density associated with (12) and c = {c(tu) :
u = 0, . . . , T} denotes the vector of unobserved states.
Since interest lies largely in the posterior distribution of parameters given aphid counts,
p(η |n), we integrate over our uncertainty for the unobserved states via MCMC. That is, we
sample the joint posterior density in equation (13) via a Gibbs sampler in which we alternate
between draws of η conditional on x = {x(tu) : u = 0, . . . , T} and draws c conditional on
the current η and n. Retaining only those sampled values of η gives a sample from the
target distribution p(η |n) (Tanner and Wong, 1987). Since the required full conditionals
typically preclude analytic tractability, we use a Metropolis-Hastings step where necessary.
Extension of this framework to the case of multiple datasets as in Section 2.3 is straight-
forward. Recall that we have twenty-seven datasets Dijk, for i, j, k ∈ {1, 2, 3}, one for
each water, nitrogen and block combination. Let xijk(tu) = (nijk(tu), cijk(tu))′ be the
vector of aphid counts and cumulative population size at time tu in dataset Dijk. Define
nijk = {nijk(tu) : u = 0, . . . , T} and cijk = {cijk(tu) : u = 0, . . . , T} to be the respec-
tive vectors of aphid counts and cumulative population sizes in dataset Dijk. Redefine
n = {nijk : i, j, k ∈ {1, 2, 3}} and c = {cijk : i, j, k ∈ {1, 2, 3}}. Finally, let η be the
thirty-six vector of parameters given in equations (10) and (11). We then obtain the joint
posterior for parameters and latent data as
p (η, c |n) ∝ p(η) p (c(t0))
3
∏
i=1
3
∏
j=1
3
∏
k=1
4
∏
u=1
p (xijk(tu) |xijk(tu−1),η) (14)

8 Gillespie and Golightly
where c(t0) denotes the vector of initial cumulative population sizes for each dataset. Note
that we take independent prior distributions for each Cijk(t0). Since the cumulative pop-
ulation size should be at least as large as the initial observed population size nijk(t0), and
we do not expect Cijk(t0) to differ hugely from nijk(t0), we specify our prior beliefs via
Cijk(t0) = nijk(t0) + ǫ where ǫ follows a Gamma distribution with shape 1 and scale 0.2.
As for the simpler case of a single dataset, sampling of (14) may proceed by alternating
between parameter draws and latent state draws. In what follows, we describe in detail the
parameter and state updates. For clarity and notational simplicity, we give details of the
MCMC scheme in the context of a single time series. To validate the methodology, we apply
the scheme to synthetic data comprising six time-series (arising from three treatments and
two blocks) in Section 4.1. We analyse the aphid data, for all twenty-seven treatment-block
combinations in Section 4.2.
Parameter updates
Since the birth and death rates λ and µ must be strictly positive, we sample their full
conditional distribution via a random walk Metropolis algorithm (O’Hagan and Forster,
2004) on the log scale. To aid the mixing of the resulting Markov chain, we sample λ and µ
in a single block. Let log(η) = (log(λ), log(µ))′ and we assume independent proper uniform
U(−10, 10) priors for log(λ) and log(µ). If at some iteration of the MCMC algorithm, the
current parameter values are λ and µ, we propose λ∗ and µ∗ jointly, where
log(λ∗) = log(λ) + ǫ1 , ǫ1 ∼ N(0, ω21)
log(µ∗) = log(µ) + ǫ2 , ǫ2 ∼ N(0, ω22)
and ǫ1 is independent of ǫ2. Note that ω1 and ω2 are tuning parameters whose values
determine the mixing properties of the chain. The values of ω1 and ω2 are chosen to achieve
an optimal acceptance probability of around 0.23 (Roberts and Rosenthal, 2001).
A move to log(η∗) is accepted with probability
min











1,
T
∏
u=1
p (x(tu) |x(tu−1),η∗)
T
∏
u=1
p (x(tu) |x(tu−1),η)











provided that the proposed values are consistent with the prior. Note that for the multiple
aphid datasets described in Section 2.3 and analysed in Section 4.2, only the baseline pa-
rameters need be strictly positive. We therefore update the remaining parameters using a
Metropolis random walk update without transforming to the log-scale.
Latent process updates
Since we only have observations on the population level, that is, N(t) and not C(t), we have
to sample the latent process at every iteration of the MCMC algorithm. Here we adopt
the simplest strategy which updates the latent process one value at a time by sampling
C(tu) from its full conditional distribution for u = 0, . . . , T . Using the joint posterior in

Generalized stochastic growth models 9
Table 1. Parameters used in the simulation study. The under-
lying (baseline) birth and death are {λ = 1.75, µ = 0.00095}.
Treatment 2 increases the death rate by 0.0004 (but leaves the
birth rate unchanged) and Treatment 3 increases the birth rate
by 0.35 (but leaves the death rate unchanged). The block ef-
fect reduces the death rate by 0.0003 (but leaves the birth rate
unchanged).
Treatment 1 Treatment 2 Treatment 3
Block 1 {1.75, 0.00095} {1.75, 0.00135} {2.1, 0.00095}
Block 2 {1.75, 0.00065} {1.75, 0.00105} {2.1, 0.00065}
equation (13), we have
p (c(tu) |x(tu−1),x(tu+1), n(tu),η) ∝ p (x(tu+1) |x(tu),η) p (x(tu) |x(tu−1),η)
= N (x(tu+1) ; ψu , Σu)N
(
x(tu) ; ψu−1 , Σu−1
)
(15)
where N(· ; ψ , Σ) denotes the probability density function of a Normal random variable
with mean ψ and covariance matrix Σ. Since ψu and Σu typically depend on C(tu)
in a nonlinear way, the form of (15) precludes analytic tractability. We therefore use a
Metropolis-Hastings step to sample the density in (15). To construct a suitable proposal
distribution we note that the data are equispaced. Hence, smoothing between the neighbours
of N(tu) we arrive at
q (x(tu) |x(tu−1),x(tu+1),η) = N
(
x(tu) ;
1
2
(x(tu−1) + x(tu+1)) ,
1
2
Σu−1
)
and condition this density on N(tu) = n(tu) using standard multivariate Normal condition-
ing arguments. We denote the resulting density by q (c(tu) |x(tu−1),x(tu+1), n(tu),η) and
use this as a proposal density inside a Metropolis-Hastings step. If the current value of the
chain is c(tu), we accept a move to c∗(tu) ∼ q (· |x(tu−1),x(tu+1), n(tu),η) with probability
min
{
1, p (x(tu+1) |x∗(tu),η) p (x∗(tu) |x(tu−1),η)p (x(tu+1) |x(tu),η) p (x(tu) |x(tu−1),η)
×
q (c(tu) |x(tu−1),x(tu+1), n(tu),η)
q (c∗(tu) |x(tu−1),x(tu+1), n(tu),η)
}
where x∗(tu) = (n(tu), c∗(tu))′. We repeat this step for each latent value in turn, thereby
updating the whole latent path. The form of q (· |x(tu−1),x(tu+1), n(tu),η) is given in
Appendix B. Other updating schemes are possible. For example we find that a random
walk also update works well, but this scheme requires the specification of additional tuning
parameters.
4. Results
4.1. Simulation Study
To examine the performance of our inference algorithm, we consider a simple example in
which we have three treatments and two blocks. Let i, j represent the block and treatment
level, i ∈ {1, 2}, j ∈ {1, 2, 3}. Denote the synthetic dataset for each block-treatment
combination by Dij . For each dataset we assume birth and death rates of the form
λij = λ+ αli + βlj and µij = µ+ αmi + βmj

10 Gillespie and Golightly
Time
Ap
hi
d
Po
pu
la
tio
n
0
500
1000
0 1 2 3 4
Treament 1
0 1 2 3 4
Treatment 2
0 1 2 3 4
Treatment 3
Fig. 2. Simulated data sets based on the parameter values in Table 1. The circles and triangle
correspond to Blocks 1 and 2 respectively.
Birth Rate
D
en
si
ty
0
2
4
6
1.6 1.7 1.8 1.9 2.0
Death Rate
D
en
si
ty
0
5000
10000
15000
20000
0.00090 0.00095 0.00100
Fig. 3. Kernel density estimates of the marginal posterior distributions of baseline birth and death
rates (λ and µ) for the synthetic data shown in Figure 2. The true values are the black vertical lines.
with αl1 = αm1 = βl1 = βm1 = 0. Note for simplicity that we assume no interaction between
treatment and block. The six datasets were simulated from the underlying exact discrete
stochastic model (and not our approximate inferential model) using the Gillespie algorithm
(Gillespie, 1977) by taking the values of λij and µij given in Table 1. Note in particular
that λ = 1.75, αl2 = 0, βl2 = 0, βl3 = 0.35, µ = 0.00095, αm2 = −0.0003, βm2 = 0.0004 and
βm3 = 0.
To mirror the real dataset in Section 2.3, the initial population sizes were chosen from a
Poisson distribution with mean 28 and data points were recorded at times 0.0, 1.14, 2.29, 3.57
and 4.57. The resulting dataset is shown in Figure 2.
We ran the MCMC algorithm described in Section 3 for two million iterations and
thinned by a factor of 1000 iterations. This yielded a sample of 2000 iterates with low
correlation to be used as the main monitoring run. The posterior densities for the baseline
birth and death rates (λ and µ) are shown in Figure 3. Plainly, the sampler gives values of
the birth and death rates that are consistent with the true values. In other simulation studies
(not reported here), in general we noticed that the death rate is harder to recover, and we
believe that this is because the parameter is strongly linked with the latent process, that
is, the reaction hazard for aphid death (with rate constant µ) depends on the unobserved
cumulative population size.
The posterior distributions for the treatment and block effects for birth and death are

Generalized stochastic growth models 11
Birth Rate
D
en
si
ty
0
2
4
6
−0.2 0.0 0.2 0.4
X
Block 2
−0.2 0.0 0.2 0.4
X
Treatment 2
−0.2 0.0 0.2 0.4
X
Treatment 3
Death Rates
D
en
si
ty
0
5000
10000
15000
20000
−4e−04 0e+00 4e−04
X
Block 2
−4e−04 0e+00 4e−04
X
Treatment 2
−4e−04 0e+00 4e−04
X
Treatment 3
Fig. 4. Kernel density estimates of the marginal posterior distributions of treatment and block birth
and death rates for the synthetic data shown in Figure 2. Left panel: αl2, βl2 and βl3 respectively. Right
panel: αm2 , βm2 and βm3 respectively. An X corresponds to the true value.
shown in Figure 4. As with the baseline parameters shown in Figure 3, samples of the birth
and death rates are clearly consistent with the true values.
4.2. Application to the Aphid Data Set
The complete aphid dataset described in Section 2.3 comprises twenty-seven treatment
combinations and involves thirty-eight parameters. We believe that this system is the largest
stochastic kinetic model (in terms of numbers of reactions and parameters) that has been
fitted to real data. The inference scheme of Section 3 was run for 2 million iterations with
a thin of 1000. This produced (near) uncorrelated samples from the parameter posterior
distributions. Marginal posterior distributions for the baseline birth and death rates are
shown in Figure 5. In Table 2 we give 95% credible intervals for all parameters. Rather
surprisingly, some of the strongest effects present are in the Block 2 interaction terms.
Further inspection of Table 2 shows that nitrogen plays an important part in the aphid
population dynamics. As levels of nitrogen were decreased (from blanket to variable to
zero), both the birth rate and death rates increased. For example, 95% credible intervals
for βl2 and βl3 (corresponding to the change in birth rate from baseline for variable and zero
nitrogen respectively) are (0.0553, 0.3220) and (0.236, 0.529). This result is in agreement
with Matis et al. (2008).
Since the purpose of this study was to determine how aphid populations are affected
by the different treatments, simply using Table 2 is insufficient. Indeed, it is difficult to
explicitly describe the effect of a particular treatment simply through inspection of Table 2,
since, for example, the birth and the death rate can increase (as with the medium water
effect) and this may have little effect on the total number of aphids. We therefore consider
the effect of different treatment combinations by examining the distributions of the total
number of aphids at future time points in Section 4.2.1.
The sensitivity of the posterior distribution to the choice of prior distribution was as-
sessed by using repeated runs of the inference scheme each with small changes in the prior
specification. We found that such changes resulted in only small changes to the posterior
distributions.
Naturally, the validity of any conclusions that can be drawn from the modelling process
can be questioned if the model does not fit the data sufficiently well. To test model fit,
that is, to see if the model is consistent with the observed data, we compared the predictive

12 Gillespie and Golightly
Table 2. Posterior 95% credible intervals for the parameters in equations (10) and (11), obtained
using the aphid dataset.
Effect Birth Rates Death Rates
Baseline λ: (1.64, 1.86) µ: (0.000904, 0.000987)
Water (M) αl2: (0.501, 0.790) αm2 : (5.73×10−5, 1.64×10−4)
Water (H) αl3: (0.0888, 0.3590) αm3 : (-0.000204, -0.000111)
Nitrogen (V) βl2: (0.0553, 0.3220) βm2 : (1.84×10−5, 1.14×10−4)
Nitrogen (Z) βl3: (0.236, 0.529) βm3 : (0.000109, 0.000208)
Block 2 ρl2: (-0.2420, 0.0112) ρm2 : (-0.000411, -0.000326)
Block 3 ρl3: (-0.169, 0.124) ρm3 : (0.000127, 0.000254)
Water (M) Nitrogen (V) (αβ)l22: (-0.3410, -0.0614) (αβ)m22: (0.000110, 0.000185)
Water (M) Nitrogen (Z) (αβ)l23: (-0.818, -0.525) (αβ)m23: (-1.43×10−4, -5.53×10−5)
Water (H) Nitrogen (V) (αβ)l32: (-0.538, -0.266) (αβ)m32: (0.000021, 0.000102)
Water (H) Nitrogen (Z) (αβ)l33: (-0.550, -0.253) (αβ)m33: (-9.21×10−5, 3.51×10−6)
Water (M) Block 2 (αρ)l22: (-0.697, -0.395) (αρ)m22: (-0.000247, -0.000144)
Water (M) Block 3 (αρ)l23: (-0.613, -0.278) (αρ)m23: (-0.000275, -0.000134)
Water (H) Block 2 (αρ)l32: (-0.2140, 0.0824) (αρ)m32: (0.000183, 0.000276)
Water (H) Block 3 (αρ)l33: (-0.0784, 0.2430) (αρ)m33: (0.000166, 0.000306)
Nitrogen (V) Block 2 (βρ)l22: (0.497, 0.781) (βρ)m22: (-0.000213, -0.000121)
Nitrogen (V) Block 3 (βρ)l23: (-0.111, 0.187) (βρ)m23: (-1.11×10−4, 2.52×10−5)
Nitrogen (Z) Block 2 (βρ)l32: (0.587, 0.884) (βρ)m32: (-1.72×10−6, 9.19×10−5)
Nitrogen (Z) Block 3 (βρ)l33: (0.0907, 0.4030) (βρ)m33: (-1.78×10−4, -4.75×10−5)
Birth Rate
D
en
si
ty
0
2
4
6
1.6 1.7 1.8 1.9 2.0
Death Rate
D
en
si
ty
0
5000
10000
15000
0.00090 0.00095 0.00100
Fig. 5. Kernel density estimates of the marginal posterior distributions of baseline birth and death
rates (λ and µ) for the aphid data shown in Figure 1.
distributions for aphid population size N(t) with the actual observations. The predictive
distribution is obtained by combining the model at a particular time point with the param-
eter posterior and integrating over the uncertainty associated with the parameter value. In
other words, the predictive distribution is taken to be the posterior average of realisations
of the population growth process, where the average is taken over the posterior distribution.
Therefore, obtaining predictive distributions is straightforward, as for each randomly sam-
pled parameter value (λijk , µijk) and the unobserved state Cijk(t0) from the MCMC output,
we can forward simulate at each observation time from the underlying discrete stochastic
model using the Gillespie algorithm. Figure 6 shows box and whisker plots of the predictive
distributions of aphid population size N(t) (for a variety of treatment scenarios) at times
1.14, 2.29, 3.57 and 4.57 weeks, plotted against the actual data. Inspection of Figure 6
suggests that the model captures the underlying characteristics of the data sufficiently well.
Further improvements to the model are discussed in Section 5.

Generalized stochastic growth models 13
Time
Ap
hi
d
Po
pu
la
tio
n
0
500
1000
1500
2000
2500
1.14 2.29 3.57 4.57
X
X
X
X
D 112
1.14 2.29 3.57 4.57
X
X
X
X
D 122
1.14 2.29 3.57 4.57
X
X
X
X
D 113
X
X
X
X
D 123
X
X
X
X
D 121
0
500
1000
1500
2000
2500
X
X
X
X
D131
Fig. 6. Predictive distributions for aphid population sizes (represented by Box and Whisker plots) and
observations for 6 aphid datasets, corresponding to 6 treatment combinations. The data points are
represented by an X.
4.2.1. Additional number of Aphids by treatment
Recall that Cijk(tn) denotes the cumulative population size at time tn in dataset Dijk.
Hence, for some future time point tn and any treatment combination ijk, we can sample
the predictive distribution of Cijk(tn) by averaging over each randomly sampled value of
λijk, µijk and Cijk(t0) (from the MCMC scheme) with fixed Nijk(t0) = nijk(t0) and forward
simulating via the Gillespie algorithm to time tn. Since we wish to compare the effects of
each treatment and their interactions, we fix t0 = 0, Cijk(0) = 1, Nijk(0) = 1 and tn = 6
weeks. Hence, for each treatment effect and interaction, we can compute the distributions
of cumulative population numbers at time 6 in addition to the numbers we would see at
baseline. For example, to obtain samples from the distribution of additional aphid numbers
due to the effect of medium water, we use each MCMC iterated value of the parameters
λ211 = λ+ αl2 and µ211 = µ+ αm2
combined with the Gillespie algorithm, to generate values of C211(6). Denote the sampled
values of C211(6) by {c(i)211(6), i = 1, . . . , N}. We then take
c(i)211(6)− c
(i)
111(6), i = 1, . . . , N
as a sample from the distribution of additional aphids due to medium water. Similar
calculations can be performed to obtain the distributions of additional aphid numbers (over
baseline) for all 18 treatment combinations. These distributions are shown in Figure 7.

14 Gillespie and Golightly
Inspection of Figure 7 suggests that the variable and zero nitrogen treatments have little
effect over baseline with distributions of additional numbers of aphids centred around zero.
Similarly, it appears that water has little effect by itself. This is not entirely unexpected
since in the year that the data were recorded, there was moderately high rainfall. Plainly,
there are strong interactions. In particular, block 2 appears to have a very strong interaction
with nitrogen.
Aphids
D
en
si
ty
0.0000
0.0005
0.0010
0.0015
0.0020
0.0025
0 2000 6000 10000
Block 3 Nitrogen (Z) Block 2 Nitrogen (Z)
0 2000 6000 10000
Block 3 Nitrogen (V) Block 2 Nitrogen (V)
0 2000 6000 10000
Block 3 Water (H) Block 2 Water (H)
Block 3 Water (M) Block 2 Water (M) Water (H) Nitrogen (Z) Water (M) Nitrogen (Z) Water (H) Nitrogen (V)
0.0000
0.0005
0.0010
0.0015
0.0020
0.0025
Water (M) Nitrogen (V)
0.0000
0.0005
0.0010
0.0015
0.0020
0.0025
Block 3
0 2000 6000 10000
Block 2 Nitrogen (Z)
0 2000 6000 10000
Nitrogen (V) Water (H)
0 2000 6000 10000
Water (M)
Fig. 7. Distribution of additional number of aphids due to each treatment over baseline at time 6.
5. Conclusion
In this paper we have considered the task of inference for a stochastic population growth
model of aphid numbers. Whilst such models are well known, interest has usually focused
on their deterministic counterparts (Matis et al. (2006, 2007)). The model considered here
is formulated as a Markov jump process with an unobserved component representing cumu-
lative population size. The associated (but unavailable) transition densities were approxi-
mated using a moment closure approach and an MCMC scheme used to sample the posterior
distribution of parameters and latent data. It is important to note that the methodology
considered here could in theory be applied to any Markov jump process (Milner et al.,
2009). The model was fitted to the data sets of Matis et al. (2008) and consisted of cotton
aphid counts for twenty-seven treatment combinations. To allow the possibility of interac-
tions, and rather than fitting a model to all datasets separately, we used a nested structure
for the birth and death parameters (specific to each treatment combination) resulting in
thirty-eight unknown parameters to be estimated.

Generalized stochastic growth models 15
As well as confirming the results of Matis et al. (2008), we also find strong interaction
effects, in particular in block 2. These effects were quantified probabilistically by calculating
the distribution of additional aphid numbers at time t = 6 over baseline. Note that in
computing these distributions, we take into account parameter uncertainty by integrating
over the unknown parameter vector.
Finally, it is worth noting that although the predictive distributions in Figure 6 suggest
adequate model fit, it may be possible to improve the model by adding in an immigration
term, resulting in
∅ γ−−→ N + C
N λ−−→ 2N + C
N + C µ−−→ C .
This is the subject of ongoing research.
Acknowledgements
We are grateful to the Associate Editor and two referees for their constructive comments and
suggestions which led to improvements in the paper. We also wish to thank Prof Richard
Boys, Dr Malcolm Farrow and Prof Jim Matis for helpful discussions.
A. Constructing the Moment Equations
Recall that the forward Kolmogorov equation governing the probability pn,c(t) is
dpn,c(t)
dt = λ(n− 1)pn−1,c−1(t) + µc(n+ 1)pn+1,c(t)− n(λ+ µc)pn,c(t), (16)
and the bivariate moment generating function is
M(θ, φ; t) ≡
∞
∑
n,c=0
enθ ecφpn,c(t) .
On multiplying equation (16) by enθecφ, we obtain
enθecφ dpn,c(t)dt = λ(n− 1)e
θeφe(n−1)θe(c−1)φpn−1,c−1(t)
+µc(n+ 1)e−θe(n+1)θecφpn+1,c(t)
−n(λ+ µc)enθecφpn,c(t).
Hence, summing both sides of the preceding equation over {n, c} gives
∂M
∂t = λ(e
θ+φ − 1)
∂M
∂θ + µ(e
−θ − 1)
∂2M
∂θ∂φ (17)
where we have dropped the notational dependence of M(θ, φ; t) on θ, φ and t. Now, we
have the cumulant generating function as
K(θ, φ; t) ≡ log[M(θ, φ; t)] =
∞
∑
n,c=0
θn
n!
φc
c! κnc ,

16 Gillespie and Golightly
and arrive at the partial differential governing K(θ, φ; t) by noting that
∂
∂t
(
eK
)
= eK ∂K∂t ,
∂
∂θ
(
eK
)
= eK ∂K∂θ ,
∂
∂θ
{
∂
∂φ
(
eK
)
}
= eK ∂
2K
∂θ∂φ + e
K ∂K
∂θ
∂K
∂φ
where we have dropped the notational dependence of K(θ, φ; t) on θ, φ and t. Hence,
substituting M = eK into equation (17) gives
∂K
∂t = λ(e
θ+φ − 1)∂K∂θ + µ(e
−θ − 1)
(
∂2K
∂θ∂φ +
∂K
∂θ
∂K
∂φ
)
.
This equation can then be used to derive ordinary differential equations for the cumulants
κ10, κ01, κ11, κ02 and κ20.
B. Constructing the Proposal Density for the Latent Data
Recall that X(tu) = (N(tu), C(tu))′ denotes the random 2-vector of population size and
cumulative size, and partition q (x(tu) |x(tu−1),x(tu+1),η) in the following way,
N
{(
n(tu)
c(tu)
)
;
1
2
(
n(tu−1) + n(tu+1)
c(tu−1) + c(tu+1)
)
, 1
2
(
σnnu−1 σncu−1
σcnu−1 σccu−1
)}
.
Then, conditioning on N(tu) = n(tu) gives the proposal density as
q (c(tu) |x(tu−1),x(tu+1), n(tu),η) = N (c(tu) ; ψ , Σ)
where
ψ = 1
2
(c(tu−1) + c(tu+1)) + σcnu−1
(
σnnu−1
)−1
(
n(tu)−
1
2
[n(tu−1) + n(tu+1)]
)
Σ =
1
2
σccu−1 −
1
2
σcnu−1
(
σnnu−1
)−1 σncu−1 .
References
Boys, R. J., D. J. Wilkinson, and T. B. L. Kirkwood (2008). Bayesian inference for a
discretely observed stochastic kinetic model. Statistics and Computing 18, 125–135.
Celini, L. and J. Vaillant (2004). A model of temporal distribution of Aphis gossypii glover
(homoptera: Aphididae) on cotton. Journal of Applied Entomology 128, 133–139.
Eraker, B. (2001). MCMC analysis of diffusion models with application to finance. Journal
of Business and Economic Statistics 19, 177–191.
Giarola, L. T. P., S. G. F. Martins, and M. C. P. Toled Costa (2006). Computer simulation
of Aphis gossypii insects using penna aging model. Physica A: Statistical Mechanics and
its Applications 368, 147–154.
Gillespie, C. S. (2009). Moment-closure approximations for mass-action models. IET Sys-
tems Biology, 3, 52–58.

Generalized stochastic growth models 17
Gillespie, D. T. (1977). Exact stochastic simulation of coupled chemical reactions. Journal
of Physical Chemistry 81, 2340–2361.
Golightly, A. and D. J. Wilkinson (2005). Bayesian inference for stochastic kinetic models
using a diffusion approximation. Biometrics 61 (3), 781–788.
Golightly, A. and D. J. Wilkinson (2006). Bayesian sequential inference for stochastic kinetic
biochemical network models. Journal of Computational Biology 13 (3), 838–851.
Golightly, A. and D. J. Wilkinson (2008). Bayesian inference for nonlinear multivariate
diffusion models observed with error. Computational Statistics & Data Analysis 52,
1674–1693.
Heron, E. A., B. Finkenstadt, and D. A. Rand (2007). Bayesian inference for dynamic
transcriptional regulation; the Hes1 system as a case study. Bioinformatics 23, 2596–
2603.
Krishnarajah, I., A. Cook, G. Marion, and G. Gibson (2005). Novel moment closure ap-
proximations in stochastic epidemics. Bulletin of Mathematical Biology 67, 855–873.
Leclant, F. and J. P. Deguine (1994). Cotton aphids. In G. A. Matthews and J. P. Tunstall
(Eds.), Insect Pests of Cotton, pp. 285–323. U.K.: CAB International.
Matis, J. H., T. R. Kiffe, T. I. Matis, J. A. Jackman, and H. Singh (2007). Population size
models based on cummulative size, with application to aphids. Ecological Modelling 205,
81–92.
Matis, J. H., T. R. Kiffe, T. I. Matis, and D. E. Stevenson (2006). Application of popula-
tion growth models based on cumulative size to Pecan aphids. Journal of Agricultural,
Biological, and Environmental Statistics 11, 425–449.
Matis, J. H., T. R. Kiffe, T. I. Matis, and D. E. Stevenson (2007). Stochastic model-
ing of aphid population growth with nonlinear power-law dynamics. Mathematical Bio-
sciences 208, 469–494.
Matis, J. H., T. R. Zhou, T. R. Kiffe, and T. I. Matis (2007). Fitting cumulative size
mechanistic models to insect population data: A nonlinear random effects model analysis.
Journal of the Indian Society of Agricultural Statistics 61, 147–155.
Matis, T. I., M. N. Parajulee, J. H. Matis, and R. B. Shrestha (2008). A mechanistic
model based analysis of cotton aphid population dynamics data. Agricultural and Forest
Entomology 10, 1–8.
Milner, P., C. S. Gillespie, and D. J. Wilkinson (2009). Parameter inference using moment
closure models. In preparation.
O’Hagan, A. and J. Forster (2004). Kendall’s Advanced Theory of Statistics (2nd ed.),
Volume Volume 2B. Arnold.
Prajneshu (1998). A nonlinear statistical model for aphid population growth. Journal of
the Indian Society of Agricultural Statistics 51, 73–80.

18 Gillespie and Golightly
Purutcuoglu, V. and E. Wit (2007). Bayesian inference of the kinetic parameters of a
realistic MAPK/ERK pathway. BMC Systems Biol. 1, P19.
Roberts, G. O. and J. Rosenthal (2001). Optimal scaling for various Metropolis-Hastings
algorithms. Statistical Science 16 (4), 351–367.
Rummel, D. R., M. D. Arnold, J. E. Slosser, K. C. Neece, and W. E. Pinchak (1995).
Cultural factors influencing the abundance of Aphis gossypii glover in Texas High Plains
cotton. Southwestern Entomologist 20, 395–406.
Singh, A. and J. Hespanha (2007). A derivative matching approach to moment closure for
the stochastic logistic model. Bulletin of Mathematical Biology 69, 1909–1925.
Tanner, M. A. and W. H. Wong (1987). The calculation of posterior distributions by data
augmentation. Journal of the American Statistical Association 82 (398), 528–540.
Wilkinson, D. J. (2006). Stochastic Modelling for Systems Biology. Chapman and Hall/CRC
Press, London.
Xia, J. Y., W. Van der Werf, and R. Rabbinge (1999). Influence of temperature on bionomics
of cotton aphid, Aphis gossypii, on cotton. Entomologia Experimentalis et Applicata 90,
25–35.
Zheng, Q. and J. Ross (1991). Comparison of deterministic and stochastic kinetics for
nonlinear systems. Journal of Chemical Physics 94, 3644–3648.