Fast Bayesian parameter estimation for stochastic logistic growth models

Jonathan Heydari Conor Lawless David A. Lydall Darren J. Wilkinson
Newcastle University, UK

Abstract

The transition density of a stochastic, logistic population growth model with multiplicative intrinsic noise is analytically intractable. Inferring model parameter values by fitting such stochastic differential equation (SDE) models to data therefore requires relatively slow numerical simulation. Where such simulation is prohibitively slow, an alternative is to use model approximations which do have an analytically tractable transition density, enabling fast inference. We introduce two such approximations, with either multiplicative or additive intrinsic noise, each derived from the linear noise approximation of the logistic growth SDE. After Bayesian inference we find that our fast LNA models, using Kalman filter recursion for computation of marginal likelihoods, give similar posterior distributions to slow arbitrarily exact models. We also demonstrate that simulations from our LNA models better describe the characteristics of the stochastic logistic growth models than a related approach. Finally, we demonstrate that our LNA model with additive intrinsic noise and measurement error best describes an example set of longitudinal observations of microbial population size taken from a typical, genome-wide screening experiment.
Keywords: Kalman Filter; Linear Noise Approximation; Logistic; Population Growth; Stochastic Modelling;

1 Introduction

Stochastic models simultaneously describe dynamics and noise or heterogeneity in real systems (Chen et al., 2010). For example, stochastic models are increasingly recognised as necessary tools for understanding the behaviour of complex biological systems (Wilkinson, 2012, 2009) and are also used to capture uncertainty in financial market behaviour (Kijima, 2013; Koller, 2012). Many such models are written as continuous stochastic differential equations (SDEs) which often do not have analytical solutions and are slow to evaluate numerically compared to their deterministic counterparts. Simulation speed is often a particularly critical issue when inferring model parameter values by comparing simulated output with observed data (Hurn et al., 2007).

For SDE models where no explicit expression for the transition density is available, it is possible to infer parameter values by simulating a latent process using a data augmentation approach (Golightly and Wilkinson, 2005). However, this method is computationally intensive and not practical for all applications. When fast inference for SDEs is important, for example for real-time analysis as part of decision support systems or for big data inference problems where we simultaneously fit models to many thousands of datasets (e.g. Heydari et al. (2012)), we need an alternative approach. Here we demonstrate one such approach: developing an analytically tractable approximation to the original SDE, by making linear noise approximations (LNAs). We apply this approach to a SDE describing logistic population growth for the first time.

The logistic model of population growth, an ordinary differential equation (ODE) describing the self-limiting growth of a population of size $X_{t}$ at time $t$ , was developed by Verhulst (1845)

\displaystyle\frac{dX_{t}}{dt}

\displaystyle=rX_{t}-\frac{r}{K}X_{t}^{2}.

(1)

The ODE has the following analytic solution:

\displaystyle X_{t}

\displaystyle=\frac{K}{1+Qe^{-r(t-t_{0})}},

(2)

where $Q=\left(\frac{K}{P}-1\right)e^{rt_{0}}$ and $P=X_{t_{0}}$ . The model describes a population growing from an initial size $P$ with an intrinsic growth rate $r$ , undergoing approximately exponential growth which slows as the availability of some critical resource (e.g. nutrients or space) becomes limiting (Turner Jr. et al., 1976). Ultimately, population density saturates at the carrying capacity (maximum achievable population density) $K$ , once the critical resource is exhausted. Where further flexibility is required, generalized forms of the logistic growth process (Tsoularis and Wallace, 2002; Peleg et al., 2007) may be used instead.

To account for uncertainty about processes affecting population growth which are not explicitly described by the deterministic logistic model, we can include a term describing intrinsic noise and consider a SDE version of the model (Li et al., 2011; Román-Román and Torres-Ruiz, 2012). Here we extend the ODE in (1) by adding a term representing multiplicative intrinsic noise (3) to give a model which we refer to as the stochastic logistic growth model (SLGM).

\displaystyle dX_{t}

\displaystyle=\left[rX_{t}-\frac{r}{K}X_{t}^{2}\right]dt+\sigma X_{t}dW_{t},

(3)

where $X_{t_{0}}=P$ and is independent of $W_{t}$ , $t\geq t_{0}$ ,

Alternative stochastic formulations of the logistic ODE can be generated (Campillo et al., 2013). The Kolmogorov forward equation has not been solved for (3) (or for any similar formulation of a logistic SDE) and so no explicit expression for the transition density is available.

Román-Román and Torres-Ruiz (2012) introduce a diffusion process approximating the SLGM (which we label RRTR) with a transition density that can be derived explicitly. The Bayesian approach can be applied in a natural way to carry out parameter inference for state space models with tractable transition densities (West and Harrison, 1997). A state space model describes the probabilistic dependence between an observation process variable $X_{t}$ and state process $S_{t}$ . The transition density is used to describe the state process $S_{t}$ and a measurement error structure is chosen to describe the relationship between $X_{t}$ and $S_{t}$ .

The Kalman filter (Kalman, 1960) is typically used to infer the hidden state process of interest $S_{t}$ and is an optimal estimator, minimising the mean square error of estimated parameters when all noise in the system can be assumed to be Gaussian. We use the Kalman filter to reduce computational time in a parameter inference algorithm by recursively computing the marginal likelihood (West and Harrison, 1997).

The RRTR can be fit to data within an acceptable time frame by assuming multiplicative measurement error to give a linear Gaussian structure, allowing us to use a Kalman filter for inference.

We introduce two new first order linear noise approximations (LNAs) (Wallace, 2010; Komorowski et al., 2009) of (3), one with multiplicative and one with additive intrinsic noise, which we label LNAM and LNAA respectively. The LNA reduces a SDE to a linear SDE with additive noise, which can be solved to give an explicit expression for the transition density. We derive transition densities for the two approximate models and construct a Kalman filter by choosing measurement noise to be either multiplicative or additive to construct a linear Gaussian structure. Exact simulations from the SGLM are compared with each of the three approximate models. We compare the utility of each of the approximate models during parameter inference by comparing simulations with both synthetic and real datasets.

2 The Román-Román and Torres-Ruiz (2012) diffusion process

Román-Román and Torres-Ruiz (2012) present a logistic growth diffusion process (RRTR) which has a transition density that can be written explicitly, allowing inference of model parameter values from discrete sampling trajectories.
The RRTR is derived from the following ODE:

\displaystyle\frac{dx_{t}}{dt}

\displaystyle=\frac{Qr}{e^{rt}+Q}x_{t}.

(4)

The solution to (4) is given in (2) (it has the same solution as (1)).

Román-Román and Torres-Ruiz (2012) see (4) as a generalisation of the Malthusian growth model with a deterministic, time-dependent fertility $h(t)=\frac{Qr}{e^{rt}+Q}$ , and replace this with $\frac{Qr}{e^{rt}+Q}+\sigma W_{t}$ to obtain the following approximation to the SLGM:

\displaystyle dX_{t}

\displaystyle=\frac{Qr}{e^{r{t}}+Q}X_{t}d{t}+{\sigma}X_{t}dW_{t},

(5)

where $Q=\left(\frac{K}{P}-1\right)e^{rt_{0}}$ , $P=X_{t_{0}}$ and is independent of $W_{t}$ , $t\geq t_{0}$ . The process described in (5) is a particular case of the lognomal process with exogenous factors, therefore an exact transition density is available (Gutiérrez et al., 2006). The transition density for $Y_{t}$ , where $Y_{t}=\log(X_{t})$ , can be written:

\displaystyle\begin{split}(Y_{t_{i}}|Y_{t_{i-1}}&=y_{t_{i-1}})\sim\operatorname{N}\left(\mu_{t_{i}},\Xi_{t_{i}}\right),\\ \text{where }a&=r,\quad b=\frac{r}{K},\\ \mu_{t_{i}}=&\log(y_{t_{i-1}})+\log(\frac{1+be^{-at_{i}}}{1+be^{-at_{i-1}}})-\frac{\sigma^{2}}{2}(t_{i}-t_{i-1})\text{ and}\\ \Xi_{t_{i}}&=\sigma^{2}(t_{i}-t_{i-1}).\end{split}

(6)

3 Linear noise approximation with multiplicative noise

We now take a different approach to approximating the SLGM (3), which will turn out to be closer to the exact solution of the SLGM than the RRTR (5). Starting from the original model (3), we apply Itô’s lemma with the transformation $f(t,X_{t})\equiv Y_{t}=\log X_{t}$ to obtain the following Itô drift-diffusion process:

\displaystyle dY_{t}=\left(r-\frac{1}{2}\sigma^{2}-\frac{r}{K}e^{Y_{t}}\right)dt+\sigma dW_{t}.

(7)

The log transformation from multiplicative to additive noise, gives a constant diffusion term, so that the LNA will give a good approximation to (3). We now separate the process $Y_{t}$ into a deterministic part $V_{t}$ and a stochastic part $Z_{t}$ so that $Y_{t}=V_{t}+Z_{t}$ and consequently $dY_{t}=dV_{t}+dZ_{t}$ . We choose $V_{t}$ to be the solution of the deterministic part

dV_{t}=\left(r-\frac{1}{2}\sigma^{2}-\frac{r}{K}e^{V_{t}}\right)dt.

After redefining our notation as follows: $a=r-\frac{\sigma^{2}}{2}$ and $b=\frac{r}{K}$ , we solve (7) for $V_{t}$ :

V_{t}=\log\left(\frac{aPe^{aT}}{bP(e^{aT}-1)+a}\right).

(8)

Differentiating w.r.t. t, we obtain:

\frac{dV_{t}}{dt}=\frac{a(a-bP)}{bP(e^{aT}-1)+a}.

Writing down an expression for $dZ_{t}$ , where $dZ_{t}=dY_{t}-dV_{t}$ .

dZ_{t}=\left(r-\frac{1}{2}\sigma^{2}-\frac{r}{K}e^{Y_{t}}\right)dt+\sigma dW_{t}-\frac{a(a-bP)}{bP(e^{aT}-1)+a}dt

and simplifying:

dZ_{t}=\left(\frac{baPe^{aT}}{bP(e^{aT}-1)+a}-be^{Y_{t}}\right)dt+\sigma dW_{t}.

Recognizing that $e^{V_{t}}=\frac{aPe^{aT}}{bP(e^{aT}-1)+a}$ :

dZ_{t}=b\left(e^{V_{t}}-e^{Y_{t}}\right)dt+\sigma dW_{t}.

We substitute in $Y_{t}=V_{t}+Z_{t}$ to give

dZ_{t}=b\left(e^{V_{t}}-e^{V_{t}+Z_{t}}\right)dt+\sigma dW_{t}.

We now apply the LNA, by making a first-order approximation of $e^{Z_{t}}\approx 1+Z_{t}$ and then simplify to give $dZ_{t}=-be^{V_{t}}Z_{t}dt+\sigma dW_{t}$ . Finally, substitute $e^{V_{t}}=\frac{aPe^{aT}}{bP(e^{aT}-1)+a}$ to obtain

dZ_{t}=-\frac{baPe^{aT}}{bP(e^{aT}-1)+a}Z_{t}dt+\sigma dW_{t}.

(9)

This process is a particular case of the time-varying Ornstein-Uhlenbeck process, which can be solved explicitly. The transition density for $Y_{t}$ (derivation in Appendix A) is then:

\displaystyle\begin{split}(Y_{t_{i}}|Y_{t_{i-1}}&=y_{t_{i-1}})\sim\operatorname{N}\left(\mu_{t_{i}},\Xi_{t_{i}}\right),\\ \text{redefine }y_{t_{i-1}}&=v_{t_{i-1}}+z_{t_{i-1}},Q=\left(\frac{\frac{a}{b}}{P}-1\right)e^{at_{0}},\\ \mu_{t_{i}}&=v_{t_{i-1}}+\log\left(\frac{1+Qe^{-at_{i-1}}}{1+Qe^{-at_{i}}}\right)+e^{-a(t_{i}-t_{i-1})}\frac{1+Qe^{-at_{i-1}}}{1+Qe^{-at_{i}}}z_{t_{i-1}}\text{ and}\\ \Xi_{t_{i}}&=\sigma^{2}\left[\frac{4Q(e^{at_{i}}-e^{at_{i-1}})+e^{2at_{i}}-e^{2at_{i-1}}+2aQ^{2}(t_{i}-t_{i-1})}{2a(Q+e^{at_{i}})^{2}}\right].\end{split}

(10)

The LNA of the SLGM with multiplicative intrinsic noise (LNAM) can then be written as

\displaystyle d\log X_{t}=\left[dV_{t}+be^{V_{t}}V_{t}-be^{V_{t}}\log X_{t}\right]dt+\sigma dW_{t},

where $P=X_{t_{0}}$ and is independent of $W_{t}$ , $t\geq t_{0}$ .
Note that the RRTR given in (5) can be similarly derived using a zero-order noise approximation ( $e^{Z_{t}}\approx 1$ ) instead of the LNA.

4 Linear noise approximation with additive noise

As in Section 3, we start from the SLGM, given in (3). Without first log transforming the process, the LNA will lead to a worse approximation to the diffusion term of the SLGM, but we will see in the coming sections that there are nevertheless advantages. We separate the process $X_{t}$ into a deterministic part $V_{t}$ and a stochastic part $Z_{t}$ so that $dX_{t}=dV_{t}+dZ_{t}$ and consequently $X_{t}=V_{t}+Z_{t}$ . We choose $V_{t}$ to be the solution of the deterministic part

dV_{t}=\left(rV_{t}-\frac{r}{K}V_{t}^{2}\right)dt.

After redefining our previous notation as follows: $a=r$ and $b=\frac{r}{K}$ , we solve $dV_{t}$ to give:

V_{t}=\frac{aPe^{aT}}{bP(e^{aT}-1)+a}.

Differentiating w.r.t. $t$ , we obtain:

\frac{dV_{t}}{dt}=\frac{a^{2}Pe^{aT}(a-bP)}{\left(bP(e^{aT}-1)+a\right)^{2}}.

We now solve $dZ_{t}$ , where $dZ_{t}=dX_{t}-dV_{t}$ . Expressions for both $dX_{t}$ and $dV_{t}$ are known:

dZ_{t}=\left(rX_{t}-\frac{r}{K}X_{t}^{2}\right)dt+\sigma X_{t}dW_{t}-\frac{a^{2}Pe^{aT}(a-bP)}{\left(bP(e^{aT}-1)+a\right)^{2}}dt.

Simplifying:

dZ_{t}=\left(aX_{t}-bX_{t}^{2}-\frac{a^{2}Pe^{aT}(a-bP)}{\left(bP(e^{aT}-1)+a\right)^{2}}\right)dt+\sigma X_{t}dW_{t}.

We then substitute in $X_{t}=V_{t}+Z_{t}$ and rearrange to give

	$\displaystyle dZ_{t}=$	$\displaystyle\left(aV_{t}-bV_{t}^{2}-\frac{a^{2}Pe^{aT}(a-bP)}{\left(bP(e^{aT}-1)+a\right)^{2}}+(a-2bV_{t})Z_{t}-bZ_{t}^{2}\right)dt$
		$\displaystyle+\left(\sigma V_{t}+\sigma Z_{t}\right)dW_{t}.$

We now apply the LNA, by setting second-order terms $-bZ^{2}dt=0$ and $\sigma Z_{t}dW_{t}=0$ to obtain

dZ_{t}=\left(aV_{t}-bV_{t}^{2}-\frac{a^{2}Pe^{aT}(a-bP)}{\left(bP(e^{aT}-1)+a\right)^{2}}+(a-2bV_{t})Z_{t}\right)dt+\sigma V_{t}dW_{t}.

(11)

This process is a particular case of the Ornstein-Uhlenbeck process, which can be solved. The transition density for $X_{t}$ (derivation in Appendix B) is then

\displaystyle\begin{split}(X_{t_{i}}|X_{t_{i-1}}&=x_{t_{i-1}})\sim N(\mu_{t_{i}},\Xi_{t_{i}}),\\ \text{where }x_{t_{i-1}}&=v_{t_{i-1}}+z_{t_{i-1}},\\ \mu_{t_{i}}&=v_{t_{i-1}}+\left(\frac{aPe^{aT_{i}}}{bP(e^{aT_{i}}-1)+a}\right)-\left(\frac{aPe^{aT_{i-1}}}{bP(e^{aT_{i-1}}-1)+a}\right)\\ &+e^{a(t_{i}-t_{i-1})}\left(\frac{bP(e^{aT_{i-1}}-1)+a}{bP(e^{aT_{i}}-1)+a}\right)^{2}Z_{t_{i-1}}\text{ and}\\ \Xi_{t}&=\frac{1}{2}\sigma^{2}aP^{2}e^{2aT_{i}}\left(\frac{1}{bP(e^{aT_{i}}-1)+a}\right)^{4}\\ &\times[b^{2}P^{2}(e^{2aT_{i}}-e^{2aT_{i-1}})+4bP(a-bP)(e^{aT_{i}}-e^{aT_{i-1}})\\ &\;\>\>\>\>+2a(t_{i}-t_{i-1})(a-bP)^{2}].\end{split}

(12)

The LNA of the SLGM, with additive intrinsic noise (LNAA) can then be written as

\displaystyle dX_{t}=\left[b{V_{t}}^{2}+\left(a-2bV_{t}\right)X_{t}\right]dt+\sigma V_{t}dW_{t},

where $P=X_{t_{0}}$ and is independent of $W_{t}$ , $t\geq t_{0}$ .

5 Simulation and Bayesian inference for logistic SDE and approximations

To test which of the three approximate models best represent the SLGM, we first compare simulated forward trajectories from the RRTR, LNAM and LNAA with simulated forward trajectories from the SLGM (Figure 1). We use the Euler-Maruyama method (Carletti, 2006) with very fine discretisation to give arbitrarily exact simulated trajectories from each SDE.

The LNAA and LNAM trajectories are visually indistinguishable from the SLGM (Figure 1 A,C & D). On the other hand, population sizes simulated with the RRTR display large deviations from the mean as the population approaches stationary phase (Figure 1A & B). Figure 1E further highlights the increases in variation as the population approaches stationary phase for simulated trajectories of the RRTR, in contrast to the SLGM and LNA models.

Refer to caption — Figure 1: Forward trajectories (No. of simulations=100) for logistic SDE and approximations. See Table 1 for parameter values. A) The stochastic logistic growth model (SLGM). B) The Román-Román and Torres-Ruiz (2012) (RRTR) approximation. C) The linear noise approximation with multiplicative intrinsic noise (LNAM). D) The linear noise approximation with additive intrinsic noise (LNAA). E) Standard deviations of simulated trajectories over time for the SLGM (black), RRTR (red), LNAM (green) and LNAA (blue).

5.1 Bayesian parameter inference with approximate models

To compare the quality of parameter inference using each of these approximations we simulated synthetic time-course data from the SLGM and combined this with either lognormal or normal measurement error. Carrying out Bayesian inference with broad priors (see (13) and (14)) we compared the parameters recovered using each approximation with those used to generate the synthetic dataset. The synthetic time-course datasets consist of 27 time points generated using the Euler-Maruyama method with very fine intervals (Kloeden and Platen, 1992).

We formulate our inference problem as a dynamic linear state space model. To allow fast parameter inference we take advantage of a linear Gaussian structure and construct a Kalman filter recursion for marginal likelihood computation (Appendix D). We therefore assume lognormal (multiplicative) error for the RRTR and LNAM, and for the LNAA we assume normal (additive) measurement error. Dependent variable $y_{t_{i}}$ and independent variable $\{t_{i},i=1,...,N\}$ are data input to the model (where $t_{i}$ is the time at point $i$ and $N$ is the number of time points). $X_{t}$ is the state process, describing the population size. See Table C.1 for prior hyper-parameter values.

For the RRTR and LNAM,

	$\displaystyle\log(y_{t_{i}})$	$\displaystyle\sim\operatorname{N}(X_{t_{i}},{\nu}^{2}),$
	$\displaystyle(X_{t_{i}}\|X_{t_{i-1}}=x_{t_{i-1}})$	$\displaystyle\sim\operatorname{N}\left(\mu_{t_{i}},\Xi_{t_{i}}\right),\text{ where }x_{t_{i}}=v_{t_{i}}+z_{t_{i}},$		(13)

$\mu_{t_{i}}$ and $\Xi_{t_{i}}$ are given by (6) and (10) for the RRTR and LNAM respectively. Priors are as follows:

	$\displaystyle\log X_{0}\equiv\log~P$	$\displaystyle\sim\operatorname{N}({\mu}_{P},{\tau_{P}}^{-1}),$	$\displaystyle\quad\quad\log~K$	$\displaystyle\sim\operatorname{N}({\mu}_{K},{\tau_{K}}^{-1}),$	$\displaystyle\quad\quad\log~r$	$\displaystyle\sim\operatorname{N}({\mu}_{r},{\tau_{r}}^{-1}),$
	$\displaystyle\log~\nu^{-2}$	$\displaystyle\sim\operatorname{N}({\mu}_{\nu},{\tau_{\nu}}^{-1}),$	$\displaystyle\log~{\sigma^{-2}}$	$\displaystyle\sim\operatorname{N}({\mu}_{\sigma},{\tau_{\sigma}}^{-1})I_{[1,\infty]}.$

For the LNAA,

	$\displaystyle y_{t_{i}}$	$\displaystyle\sim\operatorname{N}(X_{t_{i}},{\nu}^{2}),$
	$\displaystyle(X_{t_{i}}\|X_{t_{i-1}}=x_{t_{i-1}})$	$\displaystyle\sim\operatorname{N}\left(\mu_{t_{i}},\Xi_{t_{i}}\right),\text{ where }x_{t_{i}}=v_{t_{i}}+z_{t_{i}},$		(14)

$\mu_{t_{i}}$ and $\Xi_{t_{i}}$ are given by (12). Priors are as in (13). Our prior for $\log~{\sigma^{-2}}$ is truncated below 1 to avoid unnecessary exploration of extremely low probability regions, which could be caused when there are problems identifying $\nu$ , for example when $\log~\nu^{-2}$ takes large values, and to ensure that intrinsic noise does not dominate the process. The truncation limit was chosen by visual inspection.

To see how the inference from our approximate models compares with slower “exact” models, we consider Euler-Maruyama approximations (Kloeden and Platen, 1992) of (3) and of the log transformed process, using fine intervals. Given these approximations we can construct a state space model for an “exact” SLGM with lognormal measurement error (SLGM+L) and similarly for the SLGM with normal measurement error (SLGM+N), priors are as in (13).

Model fitting is carried out using standard MCMC techniques (the Gibbs sampler) (Gamerman and Lopes, 2006). Posterior means are used to obtain point estimates and standard deviations for describing variation of inferred parameters. The Heidelberger and Welch convergence diagnostic (Heidelberger and Welch, 1981) and effective sample size diagnostic (Plummer et al., 2006) are used to determine whether convergence has been achieved for all parameters.

Computational times for convergence of our MCMC schemes (code is available at https://github.com/jhncl/LNA.git) can be compared using estimates for the minimum effective sample size per second (ESS_min/sec). The average ESS_min/sec of our approximate model (coded in C) is $\sim$ 100 and “exact” model $\sim$ 1 (coded in JAGS (Plummer, 2010) with 15 imputed states between time points, chosen to maximise ESS_min/sec). We find that our C code is typically twice as fast as the simple MCMC scheme used by JAGS, indicating that our inference is ${\sim}50\times$ faster than an “exact” approach. A more efficient “exact” approach could speed up further, say by another factor of 5, but our approximate approach will at least be an order of magnitude faster. We use a burn-in of 600,000 and a thinning of 4,000 to obtain a final posterior sample size of 1,000 for MCMC convergence of all our models.

To compare the approximate models ability to recover parameters from the SLGM with simulated lognormal measurement error, we simulate data and carry out Bayesian inference. Figure 2 shows that all three approximate models can capture the synthetic time-course well, but that the RRTR model is the least representative with the largest amount of drift occurring at the saturation stage, a property not found in the SLGM or the two new LNA models. Comparing forwards trajectories with measurement error (Figure 2), the “exact” model is visually similar to all our approximate models, but least similar to the RRTR. Further, Table 1 demonstrates that parameter posterior means are close to the true values and that standard deviations are small for all models and each parameter set. By comparing posterior means and standard deviations to the true values, Table 1 shows that all our models are able to recover the three different parameter sets considered.

To compare the approximations to the SLGM with simulated normal measurement error, we simulate data and carry out Bayesian inference. Figure 3 shows that of our approximate models, only the LNAA model can appropriately represent the simulated time-course as both our models with lognormal measurement error, the RRTR and LNAM do not closely bound the data. Comparing forwards trajectories with measurement error (Figure 3), the “exact” model is most visually similar to the LNAA, which shares the same measurement error structure. Further, Table 1 demonstrates that only our models with normal measurement error have posterior means close to the true values and that standard deviations are larger in the models with lognormal measurement error. Observing the posterior means for $K$ for each parameter set (Table 1), we can see that the RRTR has the largest standard deviations and that, of the approximate models, its posterior means are furthest from both the true values and the “exact” model posterior means. Comparing LNA models to the “exact” models with matching measurement error, we can see in Table 1 that they share similar posterior means and only slightly larger standard deviations. Example posterior diagnostics given in Appendix E, demonstrate that posteriors are distributed tightly around true values for our LNAA and data from the SLGM with Normal measurement error.

Table 1: Bayesian state space model parameter posterior means, standard deviations and true values for Figure 2, 3 and 4. True values for the simulated data used for Figure 1, 2 and 3 are also given.

Panel	Model	$\hat{K}$		$\hat{r}$		$\hat{P}$		$\hat{\nu}$		$\hat{\zeta}$
Figure 2, SLGM with lognormal error
A	SGLM+L	$\numprint{0.1495369660}$	( $\numprint{0.001}$ )	$\numprint{2.9822427445}$	( $\numprint{0.01350133}$ )	$\numprint{1.001567e-04}$	( $\numprint{1.111598e-06}$ )	$\numprint{3.8597809e-03}$	( $\numprint{2.126679e-03}$ )	$\numprint{0.0173228498}$	( $\numprint{0.005037563}$ )
B	RRTR	$\numprint{0.149999335}$	( $\numprint{0.00341}$ )	$\numprint{2.99015285}$	( $\numprint{0.01092}$ )	$\numprint{9.93119e-05}$	( $\numprint{1.06863e-06}$ )	$\numprint{5.6843212256e-03}$	( $\numprint{2.36e-03}$ )	$\numprint{0.011547265557}$	( $\numprint{0.00618}$ )
C	LNAM	$\numprint{0.149614787}$	( $\numprint{0.001}$ )	$\numprint{2.98773783}$	( $\numprint{0.01318}$ )	$\numprint{9.97960e-05}$	( $\numprint{1.124219e-06}$ )	$\numprint{4.1396496436e-03}$	( $\numprint{2.179554e-03}$ )	$\numprint{0.016364920984}$	( $\numprint{0.00517}$ )
D	LNAA	$\numprint{0.14953488}$	( $\numprint{0.0005}$ )	$\numprint{3.00514737}$	( $\numprint{0.01995}$ )	$\numprint{9.64713e-05}$	( $\numprint{2.945503e-06}$ )	$\numprint{3.098832e-05}$	( $\numprint{2.53443e-05}$ )	$\numprint{0.0194997429}$	( $\numprint{0.00344}$ )
E	SGLM+L	$\numprint{0.1096583}$	( $\numprint{0.0007301528}$ )	$\numprint{3.974696}$	( $\numprint{0.04724553}$ )	$\numprint{5.054189e-05}$	( $\numprint{1.568196e-06}$ )	$\numprint{6.158735e-03}$	( $\numprint{5.527275e-03}$ )	$\numprint{0.05052818}$	( $\numprint{0.01393793}$ )
F	RRTR	$\numprint{0.1086907}$	( $\numprint{0.006533453}$ )	$\numprint{3.983695}$	( $\numprint{0.03462259}$ )	$\numprint{5.046032e-05}$	( $\numprint{1.137058e-06}$ )	$\numprint{5.928172e-03}$	( $\numprint{4.596061e-03}$ )	$\numprint{0.03704644}$	( $\numprint{0.00901516}$ )
G	LNAM	$\numprint{0.1097371}$	( $\numprint{0.0007490187}$ )	$\numprint{3.985208}$	( $\numprint{0.04606133}$ )	$\numprint{5.043313e-05}$	( $\numprint{1.580017e-06}$ )	$\numprint{6.187716e-03}$	( $\numprint{5.190529e-03}$ )	$\numprint{0.05156039}$	( $\numprint{0.0126579}$ )
H	LNAA	$\numprint{0.1097912}$	( $\numprint{0.0007867848}$ )	$\numprint{3.959346}$	( $\numprint{0.06704208}$ )	$\numprint{5.207434e-05}$	( $\numprint{4.309876e-06}$ )	$\numprint{4.539877e-05}$	( $\numprint{4.394698e-05}$ )	$\numprint{0.05912248}$	( $\numprint{0.00978808}$ )
I	SGLM+L	$\numprint{0.3003784}$	( $\numprint{0.0009748994}$ )	$\numprint{5.996904}$	( $\numprint{0.02851697}$ )	$\numprint{1.962434e-05}$	( $\numprint{4.0414e-07}$ )	$\numprint{9.543053e-03}$	( $\numprint{4.03492e-03}$ )	$\numprint{0.02406215}$	( $\numprint{0.01543261}$ )
J	RRTR	$\numprint{0.3005997}$	( $\numprint{0.0044928}$ )	$\numprint{6.015265}$	( $\numprint{0.01697153}$ )	$\numprint{1.942819e-05}$	( $\numprint{2.835413e-07}$ )	$\numprint{1.240917e-02}$	( $\numprint{2.306528e-03}$ )	$\numprint{0.007745952}$	( $\numprint{0.006357571}$ )
K	LNAM	$\numprint{0.3004755}$	( $\numprint{0.0009411441}$ )	$\numprint{6.014510}$	( $\numprint{0.03065974}$ )	$\numprint{1.953061e-05}$	( $\numprint{4.202024e-07}$ )	$\numprint{8.942722e-03}$	( $\numprint{4.251644e-03}$ )	$\numprint{0.02690633}$	( $\numprint{0.01575695}$ )
L	LNAA	$\numprint{0.3004586}$	( $\numprint{0.001146215}$ )	$\numprint{6.037006}$	( $\numprint{0.06747846}$ )	$\numprint{1.895225e-05}$	( $\numprint{1.502188e-06}$ )	$\numprint{8.122229e-05}$	( $\numprint{1.595971e-04}$ )	$\numprint{0.04733552}$	( $\numprint{0.007512705}$ )
Figure 3, SLGM with normal error
A	SLGM+N	$\numprint{0.1502746}$	( $\numprint{0.002206958}$ )	$\numprint{3.099127}$	( $\numprint{0.08534504}$ )	$\numprint{9.298825e-05}$	( $\numprint{7.304899e-06}$ )	$\numprint{5.326174e-03}$	( $\numprint{1.008844e-03}$ )	$\numprint{0.05948689}$	( $\numprint{0.02986999}$ )
B	RRTR	$\numprint{0.2127014398}$	( $\numprint{0.12323}$ )	$\numprint{1.367943557}$	( $\numprint{0.26339}$ )	$\numprint{4.55242737e-03}$	( $\numprint{2.1184933e-03}$ )	$\numprint{2.5392518626e-01}$	( $\numprint{1.0969e-01}$ )	$\numprint{0.41885683478}$	( $\numprint{0.12943}$ )
C	LNAM	$\numprint{0.1714519199}$	( $\numprint{0.03280}$ )	$\numprint{1.579950397}$	( $\numprint{0.27107}$ )	$\numprint{5.2409643e-03}$	( $\numprint{2.0479954e-03}$ )	$\numprint{2.05421122203e-01}$	( $\numprint{7.805e-02}$ )	$\numprint{0.4730700594}$	( $\numprint{0.05100}$ )
D	LNAA	$\numprint{0.150452393}$	( $\numprint{0.00222}$ )	$\numprint{2.9899791}$	( $\numprint{0.26239}$ )	$\numprint{1.189e-04}$	( $\numprint{7.098716e-05}$ )	$\numprint{5.49023564e-03}$	( $\numprint{1.06e-03}$ )	$\numprint{0.052829431599}$	( $\numprint{0.03326}$ )
E	SLGM+N	$\numprint{0.1094272}$	( $\numprint{0.0007586429}$ )	$\numprint{4.183277}$	( $\numprint{0.07356718}$ )	$\numprint{4.389757e-05}$	( $\numprint{4.128667e-06}$ )	$\numprint{9.679088e-04}$	( $\numprint{2.806267e-04}$ )	$\numprint{0.05662970}$	( $\numprint{0.01165919}$ )
F	RRTR	$\numprint{0.1574587528}$	( $\numprint{0.08749706}$ )	$\numprint{2.6313040010}$	( $\numprint{0.3367076}$ )	$\numprint{4.398197e-04}$	( $\numprint{1.677808e-04}$ )	$\numprint{1.040098786e-01}$	( $\numprint{1.008875e-01}$ )	$\numprint{0.3735663206}$	( $\numprint{0.1616813}$ )
G	LNAM	$\numprint{0.1159608439}$	( $\numprint{0.009422523}$ )	$\numprint{3.0185447100}$	( $\numprint{0.373545}$ )	$\numprint{4.967272e-04}$	( $\numprint{1.396969e-04}$ )	$\numprint{3.34648648e-02}$	( $\numprint{4.308955e-02}$ )	$\numprint{0.4753192595}$	( $\numprint{0.04357369}$ )
H	LNAA	$\numprint{0.1095810}$	( $\numprint{0.0007913119}$ )	$\numprint{4.009893}$	( $\numprint{0.1581109}$ )	$\numprint{5.011655e-05}$	( $\numprint{1.442577e-05}$ )	$\numprint{1.093103e-03}$	( $\numprint{3.638245e-04}$ )	$\numprint{0.05275491}$	( $\numprint{0.01322628}$ )
I	SLGM+N	$\numprint{0.3052235554}$	( $\numprint{0.00282216}$ )	$\numprint{5.2670916442}$	( $\numprint{0.1249732}$ )	$\numprint{3.263387e-04}$	( $\numprint{3.407446e-05}$ )	$\numprint{1.11938047e-02}$	( $\numprint{1.974384e-03}$ )	$\numprint{0.0446024422}$	( $\numprint{0.0310588}$ )
J	RRTR	$\numprint{0.313643513}$	( $\numprint{0.05677519}$ )	$\numprint{3.029806490}$	( $\numprint{0.2325609}$ )	$\numprint{1.307166e-03}$	( $\numprint{2.897121e-04}$ )	$\numprint{2.22843445e-01}$	( $\numprint{3.707941e-02}$ )	$\numprint{0.074919987}$	( $\numprint{0.08628192}$ )
K	LNAM	$\numprint{0.312629722}$	( $\numprint{0.02045301}$ )	$\numprint{3.391999960}$	( $\numprint{0.4300475}$ )	$\numprint{1.118432e-03}$	( $\numprint{3.2685e-04}$ )	$\numprint{1.17641652e-01}$	( $\numprint{8.435387e-02}$ )	$\numprint{0.360004398}$	( $\numprint{0.1647072}$ )
L	LNAA	$\numprint{0.3022532}$	( $\numprint{0.002361439}$ )	$\numprint{5.862218}$	( $\numprint{0.5228873}$ )	$\numprint{2.889783e-05}$	( $\numprint{2.598836e-05}$ )	$\numprint{8.773735e-03}$	( $\numprint{1.466338e-03}$ )	$\numprint{0.04084204}$	( $\numprint{0.02844386}$ )
Figure 4, observed yeast data
A	SLGM+L	$\numprint{0.1096916}$	( $\numprint{0.007456264}$ )	$\numprint{4.098380}$	( $\numprint{0.2994734}$ )	$\numprint{7.603410e-06}$	( $\numprint{3.205816e-06}$ )	$\numprint{3.456518e-01}$	( $\numprint{5.318567e-02}$ )	$\numprint{0.1128021}$	( $\numprint{0.109331}$ )
B	SLGM+N	$\numprint{0.1098586}$	( $\numprint{0.003273399}$ )	$\numprint{3.905401}$	( $\numprint{0.1725181}$ )	$\numprint{1.043704e-05}$	( $\numprint{3.085561e-06}$ )	$\numprint{1.851678e-04}$	( $\numprint{7.460388e-05}$ )	$\numprint{0.1666784}$	( $\numprint{0.02814871}$ )
C	RRTR	$\numprint{0.1143094}$	( $\numprint{0.02602507}$ )	$\numprint{3.763728}$	( $\numprint{0.2006525}$ )	$\numprint{1.079002e-05}$	( $\numprint{3.154911e-06}$ )	$\numprint{3.378525e-01}$	( $\numprint{4.839998e-02}$ )	$\numprint{0.07771955}$	( $\numprint{0.07732242}$ )
D	LNAM	$\numprint{0.1103702}$	( $\numprint{0.01113706}$ )	$\numprint{3.776766}$	( $\numprint{0.216272}$ )	$\numprint{1.077129e-05}$	( $\numprint{3.277218e-06}$ )	$\numprint{3.36179e-01}$	( $\numprint{5.137415e-02}$ )	$\numprint{0.1036413}$	( $\numprint{0.1079923}$ )
E	LNAA	$\numprint{0.109176}$	( $\numprint{0.003306315}$ )	$\numprint{3.832318}$	( $\numprint{0.1984635}$ )	$\numprint{1.068951e-05}$	( $\numprint{3.680477e-06}$ )	$\numprint{1.768693e-04}$	( $\numprint{6.606697e-05}$ )	$\numprint{0.1637506}$	( $\numprint{0.03259627}$ )
True values		K		r		P		$\nu$		$\sigma$
Figure 1, panels A, B, C and D		0.11		4		0.00005		N/A		0.05
Figure 2 and 3, panels A, B, C & D		0.15		3		0.0001		0.005		0.01
Figure 2 and 3, panels E, F, G and H		0.11		4		0.00005		0.001		0.05
Figure 2 and 3, panels I, J, K and L		0.3		6		0.0002		0.01		0.02

5.2 Application to observed yeast data

We now consider which diffusion equation model can best represent observed microbial population growth curves taken from a Quantitative Fitness Analysis (QFA) experiment (Addinall et al., 2011; Banks et al., 2012), see Figure 4. The data consists of scaled cell density estimates over time for budding yeast Saccharomyces cerevisiae. Independent replicate cultures are inoculated on plates and photographed over a period of 5 days. The images captured are then converted into estimates of integrated optical density (IOD, which we assume are proportional to cell population size), by the software package Colonyzer (Lawless et al., 2010). The dataset chosen for our model fitting is a representative set of 10 time-courses, each with 27 time points.

As in Figure 3, we see that the LNAA model is the only approximation that can appropriately represent the time-course and that both the RRTR and LNAM fail to bound the data as tightly as the LNAA (Figure 4). Our two “exact” models are visually similar to our approximate models with the same measurement error, with the SLGM+N most similar to the LNAA and the SLGM+L to the RRTR and LNAM. This is as expected due to matching measurement error structures. Table 1 summarises parameter estimates for the observed yeast data using each model. The variation in the the LNAA model parameter posteriors is much smaller than the RRTR and LNAM, indicating a more appropriate model fit. Comparing the LNA models and “exact” models with matching measurement error, we can see in Table 1 that they share similar posterior means and standard deviations for all parameters and in particular, they are very similar for both $K$ and $r$ , which are important phenotypes for calculating fitness (Addinall et al., 2011).

In Table 2, to compare quality of parameter inference for 10 observed yeast time-courses with each approximate model. Mean square error (MSE) for 1000 posterior sample forward simulations are calculated for each yeast time course and summed to give a Total MSE for each model. It is clear that the RRTR is the worst overall representation of the 10 yeast time courses, with the highest total MSE and a much larger total MSE than the “exact” SLGM+L. It is interesting to see there is a very similar total MSE for the SLGM+L and LNAM, and similarly for the SLGM+N and LNAA, demonstrating that our approximations perform well.

Table 2: Total mean squared error (MSE) for 10 observed yeast growth time courses, each with 1000 forward simulated time-courses with measurement error. Parameter values are taken from posterior samples. Standard Deviations give the variation between the sub-total MSEs for each yeast time course fit (n=10).

Model	SLGM+N	SLGM+L	RRTR	LNAM	LNAA
Total MSE	$\numprint{29.84698}$	$\numprint{100.165}$	$\numprint{600.6007}$	$\numprint{99.39689}$	$\numprint{30.95872}$
Standard Deviation	$\numprint{1.688649}$	$\numprint{8.390729}$	$\numprint{55.72009}$	$\numprint{9.263094}$	$\numprint{2.029686}$

6 Conclusion

We have presented two new diffusion processes for modelling logistic growth data where fast inference is required: the linear noise approximation (LNA) of the stochastic logistic growth model (SLGM) with multiplicative noise and the LNA of the SLGM with additive intrinsic noise. Both the LNAM and LNAA are derived from the linear noise approximation of the stochastic logistic growth model (SLGM). The new diffusion processes approximate the SLGM more closely than an alternative approximation (RRTR) proposed by Román-Román and Torres-Ruiz (2012). The RRTR lacks a mean reverting property that is found in the SLGM, LNAM and LNAA, resulting in increasing variance during the stationary phase of population growth (see Figure 1).

We compared the ability of each of the three approximate models and the SLGM to recover parameter values from simulated datasets using standard MCMC techniques. When modelling stochastic logistic growth with lognormal measurement error we find that our approximate models are able to represent data simulated from the original process and that the RRTR is least representative, with large variation over the stationary phase (see Figure 2). When modelling stochastic logistic growth with normal measurement error we find that only our models with normal measurement error can appropriately bound data simulated from the original process (see Figure 3). We also compared parameter posterior distribution summaries with parameter values used to generate simulated data after inference using both approximate and “exact” models (see Table 1). We find that, when using the RRTR model, posterior distributions for the carrying capacity parameter $K$ are less precise than for the LNAM and LNAA approximations. We also note that it is not possible to model additive measurement error while maintaining a linear Gaussian structure (which allows fast inference with the Kalman filter) when carrying out inference with the RRTR. We conclude that when measurement error is additive, the LNAA model is the most appropriate approximate model.

To test model performance during inference with real population data, we fitted our approximate models and the “exact” SLGM to microbial population growth curves generated by quantitative fitness analysis (QFA) (see Figure 4). We found that the LNAA model was the most appropriate for modelling experimental data. It seems likely that this is because a normal error structure best describes this particular dataset, placing the LNAM and RRTR models at a disadvantage. We demonstrate that arbitrarily exact methods and our fast approximations perform similarly during inference for 10 diverse, experimentally observed, microbial population growth curves (see Table 2) which show that, in practise, our fast approximations are as good as “exact” methods. We conclude that our LNA models are preferable to the RRTR for modelling QFA data.

It is interesting to note that, although the LNAA is not a better approximation of the original SGLM process than the LNAM, it is still quite reasonable. Figure 1A and 1D shows that the SLGM and LNAA processes are visually similar. Figure 1E demonstrates that forward trajectories of the LNAA also share similar levels of variation over time with the SLGM and LNAM.

Fast inference with the LNAA gives us the potential to develop large hierarchical Bayesian models which simultaneously describe thousands of independent time-courses from QFA with a diffusion equation, allowing us to infer the existence of genetic interactions on a genome-wide scale using realistic computational resources.

Here, we have concentrated on a biological model of population growth. However, we expect that the approach we have demonstrated: generating linear noise approximations of stochastic processes to allow fast Bayesian inference with Kalman filtering for marginal likelihood computation, will be useful in a wide range of other applications where simulation is prohibitively slow.

Appendix A LNAM Solution

First we look to solve $dZ_{t}$ , given in equation (9). We define $f(t)=-\frac{baPe^{aT}}{bP(e^{aT}-1)+a}$ to obtain the following,

dZ_{t}=f(t)Z_{t}dt+\sigma dW_{t}.

Define a new process $U_{t}=e^{-\int^{t}_{t_{0}}f(s)ds}Z_{t}$ and solve the integral,

\int^{t}_{t_{0}}f(s)ds=\int^{t}_{t_{0}}-\frac{baPe^{aS}}{bP(e^{aS}-1)+a}ds=\log\left(\frac{a}{bP(e^{aT}-1)+a}\right),

where, $S=s-{t_{0}}$ and $T=t-{t_{0}}$ . Apply the chain rule to $U_{t}$ ,

U_{t}=e^{-\int^{t}_{t_{0}}f(s)ds}dZ_{t}-f(t)e^{-\int^{t}_{t_{0}}f(s)ds}Z_{t}dt.

Now substitute in $dZ_{t}=f(t)Z_{t}dt+\sigma dW_{t}$ and simplify to give

U_{t}=e^{-\int^{t}_{t_{0}}f(s)ds}\sigma dW_{t}.

Apply the following notation $\phi(t)=e^{\int^{t}_{t_{0}}f(s)ds}=\frac{a}{bP(e^{aT}-1)+a}$ and $\psi(t)=\sigma$ to give

U_{t}=\phi(t)^{-1}\psi(t)dW_{t}.

$U_{t}$ , has the following solution,

U_{t}=U_{0}+\int^{t}_{t_{0}}\phi(s)^{-1}\psi(s)dW_{s}.

As $U_{t}=\phi(t)^{-1}Z_{t}$ , $Z_{t}$ then has the following solution,

Z_{t}=\phi(t)\left[Z_{0}+\int^{t}_{t_{0}}\phi(s)^{-1}\psi(s)dW_{s}\right].

Finally, the distribution at time t is $Z_{t}|Z_{0}\sim N(M_{t},E_{t})$ , where
$M_{t}=\phi(t)Z_{0}$ and $E_{t}=\phi(t)\int^{t}_{t_{0}}\left[{\phi(s)}^{-1}\psi(s)\right]\left[{\phi(s)}^{-1}\psi(s)\right]^{T}ds\>\phi(t)^{T}$ .
Further, $M_{t}=\frac{a}{bP(e^{aT}-1)+a}Z_{0}$ and $E_{t}=\sigma^{2}\left[\frac{a}{bP(e^{aT}-1)+a}\right]^{2}\int^{t}_{t_{0}}\left[\frac{a}{bP(e^{aS}-1)+a}\right]^{-2}$ .
As $\int^{t}_{t_{0}}\left[\frac{a}{bP(e^{aS}-1)+a}\right]^{-2}ds=\frac{b^{2}P^{2}e^{2aT}+4bP(a-bP)e^{aT}+2at(a-bP)^{2}}{2a^{3}}-\frac{b^{2}P^{2}+4bP(a-bP)+2at_{0}(a-bP)^{2}}{2a^{3}}$ ,

	$\displaystyle E_{t}=$	$\displaystyle\sigma^{2}\left[\frac{a}{bP(e^{aT}-1)+a}\right]^{2}\left[\frac{b^{2}P^{2}(e^{2aT}-1)+4bP(a-bP)(e^{aT}-1)+2aT(a-bP)^{2}}{2a^{3}}\right]$
	$\displaystyle=$	$\displaystyle\sigma^{2}\left[\frac{b^{2}P^{2}(e^{2aT}-1)+4bP(a-bP)(e^{aT}-1)+2aT(a-bP)^{2}}{2a\left(bP(e^{aT}-1)+a\right)^{2}}\right].$

Taking our solutions for $V_{t}$ (8) and $Z_{t}$ , we can now write our solution for the LNA to the log of the logistic growth process (7).
As $Y_{t}=V_{t}+Z_{t}$ ,

Y_{t}|Y_{0}\sim\mathcal{N}\left(\log\left[\frac{aPe^{aT}}{bP(e^{aT}-1)+a}\right]+M_{t},E_{t}\right).

Note: $\frac{aPe^{aT}}{bP(e^{aT}-1)+a}$ has the same functional form as the solution to the deterministic part of the logistic growth process (3) and is equivalent when $\sigma=0$ (such that $a=r-\frac{\sigma^{2}}{2}=r$ ).

Further, as $Y_{t}$ is normally distributed, we know $X_{t}=e^{Y_{t}}$ will be log normally distributed and

X_{t}|X_{0}\sim\log\>\mathcal{N}(\log\left(\frac{aPe^{aT}}{bP(e^{aT}-1)+a}\right)+M_{t},E_{t}).

Alternatively set $Q=\left(\frac{\frac{a}{b}}{P}-1\right)e^{at_{0}}$ ,

X_{t}|X_{0}\sim\log\>\mathcal{N}(\log\left(\frac{\frac{a}{b}}{1+Qe^{-at}}\right)+M_{t},E_{t}).

From our solution to the log process we can obtain the following transition density

\displaystyle\begin{split}(Y_{t_{i}}|Y_{t_{i-1}}&=y_{t_{i-1}})\sim\operatorname{N}\left(\mu_{t_{i}},\Xi_{t_{i}}\right),\\ \text{where }y_{t_{i-1}}&=v_{t_{i-1}}+z_{t_{i-1}},Q=\left(\frac{\frac{a}{b}}{P}-1\right)e^{at_{0}},\\ \mu_{t_{i}}&=v_{t_{i-1}}+\log\left(\frac{1+Qe^{-at_{i-1}}}{1+Qe^{-at_{i}}}\right)+e^{-a(t_{i}-t_{i-1})}\frac{1+Qe^{-at_{i-1}}}{1+Qe^{-at_{i}}}z_{t_{i-1}}\text{ and}\\ \Xi_{t_{i}}&=\sigma^{2}\left[\frac{4Q(e^{at_{i}}-e^{at_{i-1}})+e^{2at_{i}}-e^{2at_{i-1}}+2aQ^{2}(t_{i}-t_{i-1})}{2a(Q+e^{at_{i}})^{2}}\right].\end{split}

Appendix B LNAA Solution

First we look to solve $dZ_{t}$ , given in (11). We define $f(t)=a-2bV_{t}$ and $g(t)=aV_{t}-bV_{t}^{2}-\frac{a^{2}Pe^{aT}(a-bP)}{\left(bP(e^{aT}-1)+a\right)^{2}}$ to obtain the following,

dZ_{t}=\left(g(t)+f(t)Z_{t}\right)dt+\sigma V_{t}dW_{t}.

Define a new process $U_{t}=e^{-\int^{t}_{t_{0}}f(s)ds}Z_{t}$ and solve the integral,

\int^{t}_{t_{0}}f(s)ds=\int^{t}_{t_{0}}(a-2bV_{s})ds=aT-2\log\left(\frac{bP(e^{aT}-1)+a}{a}\right),

as $\int^{t}_{t_{0}}V_{s}ds=\frac{1}{b}\log\left(\frac{bP(e^{aT}-1)+a}{a}\right)$ , where $S=s-{t_{0}}$ and $T=t-{t_{0}}$ . Apply the chain rule to $U_{t}$ ,

U_{t}=e^{-\int^{t}_{t_{0}}f(s)ds}dZ_{t}-f(t)e^{-\int^{t}_{t_{0}}f(s)ds}Z_{t}dt.

Now substitute in $dZ_{t}=\left(g(t)+f(t)Z_{t}\right)dt+\sigma V_{t}dW_{t}$ and simplify to give,

U_{t}=e^{-\int^{t}_{t_{0}}f(s)ds}g(t)dt+e^{-\int^{t}_{t_{0}}f(s)ds}\sigma V_{t}dW_{t}.

Apply the following notation $\phi(t)=e^{\int^{t}_{t_{0}}f(s)ds}=e^{aT}\left(\frac{a}{bP(e^{aT}-1)+a}\right)^{2}$ and $\psi(t)=\sigma V_{t}$ to give,

U_{t}=\phi(t)^{-1}g(t)+\phi(t)^{-1}\psi(t)dW_{t}.

$U_{t}$ has the following solution,

U_{t}=U_{0}+\int^{t}_{t_{0}}\phi(s)^{-1}g(s)ds+\int^{t}_{t_{0}}\phi(s)^{-1}\psi(s)dW_{s}.

As $U_{t}=\phi(t)^{-1}Z_{t}$ , $Z_{t}$ has the following solution,

Z_{t}=\phi(t)\left[Z_{0}+\int^{t}_{t_{0}}\phi(s)^{-1}g(s)ds+\int^{t}_{t_{0}}\phi(s)^{-1}\psi(s)dW_{s}\right].

Finally the distribution at time t is $Z_{t}|Z_{0}\sim N(M_{t},E_{t})$ , where
$M_{t}=\phi(t)\left(Z_{0}+\int^{t}_{t_{0}}\phi(s)^{-1}g(s)ds\right)$ and $E_{t}=\phi(t)\int^{t}_{t_{0}}\left[{\phi(t)}^{-1}\psi(t)\right]\left[{\phi(t)}^{-1}\psi(t)\right]^{T}ds\>\phi(t)^{T}$ .

	$\displaystyle M_{t}$	$\displaystyle=\phi(t)\left(Z_{0}+\int^{t}_{t_{0}}\left[\left(\phi(t)\right)^{-1}\left(aV_{s}-bV_{s}^{2}-\frac{a^{2}Pe^{aS}(a-bP)}{\left(bP(e^{aS}-1)+a\right)^{2}}\right)\right]ds\right)$
		$\displaystyle=\phi(t)\left(Z_{0}+\int^{t}_{t_{0}}\left[\left(\phi(t)\right)^{-1}\left(0\right)\right]ds\right)$
		$\displaystyle=e^{aT}\left(\frac{a}{bP(e^{aT}-1)+a}\right)^{2}Z_{0}$

and

E_{t}=\left(e^{aT}\left(\frac{a}{bP(e^{aT}-1)+a}\right)^{2}\right)^{2}\int^{t}_{t_{0}}\left[e^{aS}\left(\frac{a}{bP(e^{aS}-1)+a}\right)^{2}\right]^{-2}\sigma^{2}V_{s}^{2}ds

=\sigma^{2}\left(e^{aT}\left(\frac{a}{bP(e^{aT}-1)+a}\right)^{2}\right)^{2}\int^{t}_{t_{0}}\left[e^{aS}\left(\frac{a}{bP(e^{aS}-1)+a}\right)^{2}\right]^{-2}\left[\frac{aPe^{aS}}{bP(e^{aS}-1)+a}\right]^{2}ds.

=\sigma^{2}\left(e^{aT}\left(\frac{a}{bP(e^{aT}-1)+a}\right)^{2}\right)^{2}\int^{t}_{t_{0}}\left[e^{-2aS}\left(\frac{a}{bP(e^{aS}-1)+a}\right)^{-4}\right]\left[\frac{aPe^{aS}}{bP(e^{aS}-1)+a}\right]^{2}ds

=\sigma^{2}\left(e^{aT}\left(\frac{1}{bP(e^{aT}-1)+a}\right)^{2}\right)^{2}\int^{t}_{t_{0}}\left[a^{2}P^{2}\left(\frac{1}{bP(e^{aS}-1)+a}\right)^{-2}\right]ds,

as $\int^{t_{i}}_{t_{0}}\left(\frac{1}{bP(e^{aS}-1)+a}\right)^{-2}ds=\frac{b^{2}P^{2}e^{2aT}+4bP(a-bP)e^{aT}+2at(a-bP)^{2}}{2a}-\frac{b^{2}P^{2}+4bP(a-bP)+2at_{0}(a-bP)^{2}}{2a}$ ,

	$\displaystyle E_{t}=$	$\displaystyle\frac{1}{2}\sigma^{2}aP^{2}e^{2aT}\left(\frac{1}{bP(e^{aT}-1)+a}\right)^{4}$
		$\displaystyle\times\left[b^{2}P^{2}(e^{2aT}-1)+4bP(a-bP)(e^{aT}-1)+2aT(a-bP)^{2}\right].$

From our solution to the process we can obtain the following transition density

\displaystyle\begin{split}(X_{t_{i}}|X_{t_{i-1}}=&x_{t_{i-1}})\sim N(\mu_{t_{i}},\Xi_{t_{i}}),\\ \text{where }x_{t_{i-1}}=&v_{t_{i-1}}+z_{t_{i-1}},\\ \mu_{t_{i}}=&v_{t_{i-1}}+\left(\frac{aPe^{aT_{i}}}{bP(e^{aT_{i}}-1)+a}\right)-\left(\frac{aPe^{aT_{i-1}}}{bP(e^{aT_{i-1}}-1)+a}\right)\\ &+e^{a(t_{i}-t_{i-1})}\left(\frac{bP(e^{aT_{i-1}}-1)+a}{bP(e^{aT_{i}}-1)+a}\right)^{2}Z_{t_{i-1}}\text{ and}\\ \Xi_{t}=&\frac{1}{2}\sigma^{2}aP^{2}e^{2aT_{i}}\left(\frac{1}{bP(e^{aT_{i}}-1)+a}\right)^{4}\\ &\times[b^{2}P^{2}(e^{2aT_{i}}-e^{2aT_{i-1}})+4bP(a-bP)(e^{aT_{i}}-e^{aT_{i-1}})\\ &\;\>\>\>\>+2a(t_{i}-t_{i-1})(a-bP)^{2}].\end{split}

Appendix C Prior Hyper-parameters for Bayesian State Space Models

Table C.1: Prior hyper-parameters for Bayesian sate space models, Lognormal with mean (

\mu

) and precision (

\tau

)

Parameter Name	Value
${\mu}_{K}$	$\log(0.1)$
$\tau_{K}$	2
${\mu}_{r}$	$\log(3)$
$\tau_{r}$	5
${\mu}_{P}$	$\log(0.0001)$
$\tau_{P}$	0.1
${\mu}_{\sigma}$	$\log(100)$
$\tau_{\sigma}$	0.1
${\mu}_{\nu}$	$\log(10000)$
$\tau_{\nu}$	0.1

Appendix D Kalman Filter

To find $\pi(y_{t_{1:N}})$ for the LNAA with normal measurement error we can use the following Kalman Filter algorithm. First we assume the following:

	$\displaystyle\theta_{{t_{i}}}\|y_{1:{{t_{i}}}}$	$\displaystyle\sim\operatorname{N}(m_{t_{i}},C_{t_{i}}),$
	$\displaystyle m_{t_{i}}$	$\displaystyle=a_{t_{i}}+R_{t_{i}}F(F^{T}R_{{t_{i}}}F+U)^{-1}[y_{t_{i}}-F^{T}a_{t_{i}}],$
	$\displaystyle C_{t_{i}}$	$\displaystyle=R_{t_{i}}-R_{t_{i}}F(F^{T}R_{t_{i}}F+U)^{-1}F^{T}R_{t_{i}}$

and initialize with $m_{0}=P$ and $C_{0}=0$ . Now suppose that,

	$\displaystyle\theta_{t_{i}}\|y_{1:{t_{i-1}}}$	$\displaystyle\sim\operatorname{N}(a_{t_{i}},R_{t_{i}}),$
	$\displaystyle a_{t_{i}}$	$\displaystyle=G_{{t_{i}}}m_{{t_{i-1}}}$
	$\displaystyle\text{and }R_{t_{i}}$	$\displaystyle=G_{{t_{i}}}C_{{t_{i-1}}}G_{t_{i}}^{T}+W_{t_{i}}.$

The transition density distribution, see (12) is as follows:

	$\displaystyle\theta_{{t_{i}}}\|\theta_{{t_{i-1}}}$	$\displaystyle\sim\operatorname{N}(G_{{t_{i}}}\theta_{{t_{i-1}}},W_{t_{i}})$
	$\displaystyle\text{or equivalently }(X_{t_{i}}\|X_{t_{i-1}}=x_{t_{i-1}})$	$\displaystyle\sim\operatorname{N}\left(\mu_{t_{i}},\Xi_{t_{i}}\right),\text{ where }x_{t_{i-1}}=v_{t_{i-1}}+z_{t_{i-1}},$
	$\displaystyle\theta_{t}$	$\displaystyle=\begin{pmatrix}1\\ X_{t_{i}}\end{pmatrix}=\begin{pmatrix}1&0\\ H_{\alpha,t_{i}}&H_{\beta,t_{i}}\end{pmatrix}\begin{pmatrix}1\\ X_{t_{i-1}}\end{pmatrix}$
		$\displaystyle=G_{t_{i}}\theta_{{t_{i-1}}},$
	$\displaystyle G_{t_{i}}$	$\displaystyle=\begin{pmatrix}1&0\\ H_{\alpha,t_{i}}&H_{\beta,t_{i}}\end{pmatrix},\quad W_{t_{i}}=\begin{pmatrix}0&0\\ 0&\Xi_{t_{i}}\end{pmatrix}$
	$\displaystyle\text{where }H_{\alpha,t_{i}}=H_{\alpha}({t_{i}},{t_{i-1}})=$	$\displaystyle V_{t}-V_{t-1}e^{a(t_{i}-t_{i-1})}\left(\frac{bP(e^{aT_{i-1}}-1)+a}{bP(e^{aT_{i}}-1)+a}\right)^{2}$
	$\displaystyle\text{and }H_{\beta,t_{i}}=$	$\displaystyle H_{\beta}({t_{i}},{t_{i-1}})=e^{a(t_{i}-t_{i-1})}\left(\frac{bP(e^{aT_{i-1}}-1)+a}{bP(e^{aT_{i}}-1)+a}\right)^{2}.$

The measurement error distribution is as follows:

	$\displaystyle y_{t_{i}}\|\theta_{{t_{i}}}{\sim}$	$\displaystyle\operatorname{N}(F^{T}\theta_{{t_{i}}},U)$
	$\displaystyle\text{or equivalently }y_{t_{i}}\|\theta_{{t_{i}}}{\sim}$	$\displaystyle\operatorname{N}(X_{{t_{i}}},\sigma_{\nu}^{2}),$
	$\displaystyle\text{where }F=$	$\displaystyle\begin{pmatrix}0\\ 1\end{pmatrix}\text{ and }U=\sigma_{\nu}^{2}.$

Matrix Algebra:

	$\displaystyle a_{t_{i}}=$	$\displaystyle G_{{t_{i}}}m_{{t_{i-1}}}$
	$\displaystyle=$	$\displaystyle\begin{pmatrix}1&0\\ H_{\alpha,t_{i}}&H_{\beta,t_{i}}\end{pmatrix}\begin{pmatrix}1\\ m_{{t_{i-1}}}\end{pmatrix}=\begin{pmatrix}1\\ H_{\alpha,t_{i}}+H_{\beta,t_{i}}m_{{t_{i-1}}}\end{pmatrix}$

	$\displaystyle R_{t_{i}}$	$\displaystyle=G_{{t_{i}}}C_{{t_{i-1}}}G_{t_{i}}^{T}+W_{t_{i}}$
		$\displaystyle=\begin{pmatrix}0&0\\ 0&{H_{\beta,t_{i}}}^{2}c_{{t_{i-1}}}^{2}\end{pmatrix}+\begin{pmatrix}0&0\\ 0&\Xi_{t_{i}}\end{pmatrix}=\begin{pmatrix}0&0\\ 0&{H_{\beta,t_{i}}}^{2}c_{{t_{i-1}}}^{2}+\Xi_{t_{i}}\end{pmatrix}$

\displaystyle C_{t_{i-1}}

\displaystyle=\begin{pmatrix}0&0\\ 0&c_{{t_{i-1}}}^{2}\end{pmatrix}

	$\displaystyle R_{t_{i}}F(F^{T}R_{{t_{i}}}F+U)^{-1}=$	$\displaystyle\begin{pmatrix}0&0\\ 0&{H_{\beta,t_{i}}}^{2}c_{{t_{i-1}}}^{2}+\Xi_{t_{i}}\end{pmatrix}\begin{pmatrix}0\\ 1\end{pmatrix}$
		$\displaystyle\times\left[\begin{pmatrix}0&1\\ \end{pmatrix}\begin{pmatrix}0&0\\ 0&{H_{\beta,t_{i}}}^{2}c_{{t_{i-1}}}^{2}+\Xi_{t_{i}}\end{pmatrix}\begin{pmatrix}0\\ 1\end{pmatrix}+\sigma_{\nu}^{2}\right]^{-1}$
	$\displaystyle=$	$\displaystyle\left[\begin{pmatrix}{H_{\beta,t_{i}}}^{2}c_{{t_{i-1}}}^{2}+\Xi_{t_{i}}+\sigma_{\nu}^{2}\end{pmatrix}\right]^{-1}\begin{pmatrix}0\\ {H_{\beta,t_{i}}}^{2}c_{{t_{i-1}}}^{2}+\Xi_{t_{i}}\end{pmatrix}$

	$\displaystyle m_{t_{i}}=$	$\displaystyle a_{t_{i}}+R_{t_{i}}F(F^{T}R_{{t_{i}}}F+U)^{-1}[y_{t_{i}}-F^{T}a_{t_{i}}]$
	$\displaystyle=$	$\displaystyle\begin{pmatrix}1\\ H_{\alpha,t_{i}}+H_{\beta,t_{i}}m_{{t_{i-1}}}\end{pmatrix}$
		$\displaystyle+\left[\begin{pmatrix}{H_{\beta,t_{i}}}^{2}c_{{t_{i-1}}}^{2}+\Xi_{t_{i}}+\sigma_{\nu}^{2}\end{pmatrix}\right]^{-1}\begin{pmatrix}0\\ {H_{\beta,t_{i}}}^{2}c_{{t_{i-1}}}^{2}+\Xi_{t_{i}}\end{pmatrix}\left[y_{t_{i}}-\begin{pmatrix}0&1\\ \end{pmatrix}\begin{pmatrix}1\\ H_{\alpha,t_{i}}+H_{\beta,t_{i}}m_{{t_{i-1}}}\end{pmatrix}\right]$
	$\displaystyle=$	$\displaystyle\begin{pmatrix}0\\ H_{\alpha,t_{i}}+H_{\beta,t_{i}}m_{{t_{i-1}}}+\frac{{H_{\beta,t_{i}}}^{2}c_{{t_{i-1}}}^{2}+\Xi_{t_{i}}}{{H_{\beta,t_{i}}}^{2}c_{{t_{i-1}}}^{2}+\Xi_{t_{i}}+\sigma_{\nu}^{2}}\left[y_{t_{i}}-H_{\alpha,t_{i}}-H_{\beta,t_{i}}m_{{t_{i-1}}}\right]\end{pmatrix}$

	$\displaystyle C_{t_{i}}=$	$\displaystyle R_{t_{i}}-R_{t_{i}}F(F^{T}R_{t_{i}}F+U)^{-1}F^{T}R_{t_{i}}$
	$\displaystyle=$	$\displaystyle\begin{pmatrix}0&0\\ 0&{H_{\beta,t_{i}}}^{2}c_{{t_{i-1}}}^{2}+\Xi_{t_{i}}\end{pmatrix}$
		$\displaystyle-\left[\begin{pmatrix}{H_{\beta,t_{i}}}^{2}c_{{t_{i-1}}}^{2}+\Xi_{t_{i}}+\sigma_{\nu}^{2}\end{pmatrix}\right]^{-1}\begin{pmatrix}0\\ {H_{\beta,t_{i}}}^{2}c_{{t_{i-1}}}^{2}+\Xi_{t_{i}}\end{pmatrix}\left[\begin{pmatrix}0&1\end{pmatrix}\begin{pmatrix}0&0\\ 0&{H{\beta,t_{i}}}^{2}c_{{t_{i-1}}}^{2}+\Xi_{t_{i}}\end{pmatrix}\right]$
	$\displaystyle=$	$\displaystyle\begin{pmatrix}0&0\\ 0&{H_{\beta,t_{i}}}^{2}c_{{t_{i-1}}}^{2}+\Xi_{t_{i}}-\frac{\left({H_{\beta,t_{i}}}^{2}c_{{t_{i-1}}}^{2}+\Xi_{t_{i}}\right)^{2}}{{H_{\beta,t_{i}}}^{2}c_{{t_{i-1}}}^{2}+\Xi_{t_{i}}+\sigma_{\nu}^{2}}\end{pmatrix}$

With $m_{t_{i}}$ and $C_{t_{i}}$ for $i=1:T$ , we can evaluate $a_{t_{i}}$ , $R_{t_{i}}$ and $\pi(x_{t_{i}}|y_{t_{1:(i-1)}})$ for $i=1:T$ . We are interested in $\pi(y_{t_{1:i}})=\prod^{T}_{i=1}\pi(y_{t_{i}}|y_{t_{1:(i-1)}})$ , where $\pi(y_{t_{i}}|y_{t_{1:(i-1)}})=\int_{x}\pi(y_{t_{i}}|x_{t_{i}})\pi(x_{t_{i}}|y_{t_{1:(i-1)}})dx_{t_{i}}$ gives a tractable Gaussian integral. Finally,

	$\displaystyle\log\pi(y_{{t_{1:(i-1)}}})$	$\displaystyle=\sum^{T}_{i=1}\log\pi(y_{t_{i}}\|y_{{t_{1:(i-1)}}})$
		$\displaystyle=\sum^{T}_{i=1}\left[-\log\left({\sqrt{2\pi(\sigma_{f}^{2}+\sigma_{g}^{2})}}\right){-\frac{(\mu_{f}-\mu_{g})^{2}}{2(\sigma_{f}^{2}+\sigma_{g}^{2})}}\right],$

	$\displaystyle\text{where }\mu_{f}-\mu_{g}=$	$\displaystyle y_{t_{i}}-a_{t_{i}}=y_{t_{i}}-H_{\alpha,t_{i}}-H_{\beta,t_{i}}m_{t_{i-1}}$
	$\displaystyle\text{and }\sigma_{f}^{2}+\sigma_{g}^{2}=$	$\displaystyle\sigma_{\nu}^{2}+R_{t_{i}}=\sigma_{\nu}^{2}+{H_{\beta,t_{i}}}^{2}c_{{t_{i-1}}}^{2}+\Xi_{t_{i}}.$

Procedure

1. Set $i=1$ . Initialize $m_{0}=\log(P)$ and $C_{0}=0$ .

2. Evaluate and store the following log likelihood term:

	$\displaystyle\log\pi(y_{t_{i}}\|y_{t_{1:(i-1)}})=$	$\displaystyle\left[-\log\left({\sqrt{2\pi(\sigma_{f}^{2}+\sigma_{g}^{2})}}\right){-\frac{(\mu_{f}-\mu_{g})^{2}}{2(\sigma_{f}^{2}+\sigma_{g}^{2})}}\right],$
	$\displaystyle\text{where }\mu_{f}-\mu_{g}=$	$\displaystyle y_{t_{i}}-H_{\alpha,t_{i}}-H_{\beta,t_{i}}m_{{t_{i-1}}}\text{ and }\sigma_{f}^{2}+\sigma_{g}^{2}=\sigma_{\nu}^{2}+{H_{\beta,t_{i}}}^{2}c_{{t_{i-1}}}^{2}+\Xi_{t_{i}}.$

3. Create and store both $m_{t_{i}}$ , and $C_{t_{i}}$ ,

	$\displaystyle\text{where }m_{t_{i}}=$	$\displaystyle H_{\alpha,t_{i}}+H_{\beta,t_{i}}m_{t_{i-1}}+\frac{{H_{\beta,t_{i}}}^{2}c_{{t_{i-1}}}^{2}+\Xi_{t_{i}}}{{H_{\beta,t_{i}}}^{2}c_{{t_{i-1}}}^{2}+\Xi_{t_{i}}+\sigma_{\nu}^{2}}\left[y_{t_{i}}-H_{\alpha,t_{i}}-H_{\beta,t_{i}}m_{{t_{i-1}}}\right]$
	$\displaystyle\text{and }c_{{t_{i}}}^{2}=$	$\displaystyle{H_{\beta,t_{i}}}^{2}c_{{t_{i-1}}}^{2}+\Xi_{t_{i}}-\frac{\left({H_{\beta,t_{i}}}^{2}c_{{t_{i-1}}}^{2}+\Xi_{t_{i}}\right)^{2}}{{H_{\beta,t_{i}}}^{2}c_{{t_{i-1}}}^{2}+\Xi_{t_{i}}+\sigma_{\nu}^{2}}.$

4. Increment $i$ , $i$ = $(i+1)$ and repeat steps 2-3 till $\log\pi(y_{t_{N}}|y_{t_{1:(N-1)}})$ is evaluated.

5. Calculate the sum:

\log\pi(y_{{t_{1:N}}})=\sum^{N}_{i=1}\log\pi(y_{t_{i}}|y_{{t_{1:(i-1)}}}).

Appendix E Parameter posterior means

Acknowledgments

This research was supported by grants from both the Biotechnology and Biological Sciences Research Council and the Medical Research Council.

References

Addinall et al. (2011) Addinall, S. G., E.-M. Holstein, C. Lawless, M. Yu, K. Chapman, A. P. Banks, H.-P. Ngo, L. Maringele, M. Taschuk, A. Young, A. Ciesiolka, A. L. Lister, A. Wipat, D. J. Wilkinson, and D. Lydall (2011). Quantitative Fitness Analysis Shows That NMD Proteins and Many Other Protein Complexes Suppress or Enhance Distinct Telomere Cap Defects. PLoS Genet 7(4), e1001362.
Banks et al. (2012) Banks, A., C. Lawless, and D. Lydall (2012). A Quantitative Fitness Analysis Workflow. J. Vis. Exp 66, e4018.
Campillo et al. (2013) Campillo, F., M. Joannides, and I. Larramendy-Valverde (2013). Estimation of the parameters of a stochastic logistic growth model. In submission, (available at, http://arxiv.org/abs/1307.2217).
Carletti (2006) Carletti, M. (2006). Numerical solution of stochastic differential problems in the biosciences. J. Comput. Appl. Math. 185(2), 422–440.
Chen et al. (2010) Chen, Y., C. Lawless, C. S. Gillespie, J. Wu, R. J. Boys, and D. J. Wilkinson (2010). Calibayes and basis: integrated tools for the calibration, simulation and storage of biological simulation models. Briefings in Bioinformatics 11(3), 278–289.
Gamerman and Lopes (2006) Gamerman, D. and H. Lopes (2006). Markov Chain Monte Carlo: Stochastic Simulation for Bayesian Inference. Chapman and Hall/CRC Texts in Statistical Science Series. Taylor & Francis.
Golightly and Wilkinson (2005) Golightly, A. and D. J. Wilkinson (2005). Bayesian inference for stochastic kinetic models using a diffusion approximation. Biometrics 61(3), 781–788.
Gutiérrez et al. (2006) Gutiérrez, R., N. Rico, P. Román-Román, and F. Torres-Ruiz (2006). Approximate and generalized confidence bands for some parametric functions of the lognormal diffusion process with exogenous factors. Scientiae Mathematicae Japonicae 64(2), 313–330.
Heidelberger and Welch (1981) Heidelberger, P. and P. D. Welch (1981). A spectral method for confidence interval generation and run length control in simulations. Commun. ACM 24(4), 233–245.
Heydari et al. (2012) Heydari, J. J., C. Lawless, D. A. Lydall, and D. J. Wilkinson (2012). Bayesian hierarchical modelling for inferring genetic interactions in yeast. In submission.
Hurn et al. (2007) Hurn, A. S., J. I. Jeisman, and K. A. Lindsay (2007). Seeing the Wood for the Trees: A Critical Evaluation of Methods to Estimate the Parameters of Stochastic Differential Equations. Journal of Financial Econometrics 5(3), 390–455.
Kalman (1960) Kalman, R. E. (1960). A new approach to linear filtering and prediction problems. Transactions of the ASME–Journal of Basic Engineering 82(Series D), 35–45.
Kijima (2013) Kijima, M. (2013). Stochastic Processes with Applications to Finance, Second Edition. Chapman & Hall/CRC financial mathematics series. CRC Press.
Kloeden and Platen (1992) Kloeden, P. and E. Platen (1992). Numerical Solution of Stochastic Differential Equations. Applications of mathematics : stochastic modelling and applied probability. Springer.
Koller (2012) Koller, M. (2012). Stochastic Models in Life Insurance. EAA Series. Springer.
Komorowski et al. (2009) Komorowski, M., B. Finkenstadt, C. V. Harper, and D. A. Rand (2009). Bayesian inference of biochemical kinetic parameters using the linear noise approximation. BMC Bioinformatics 10, 343.
Lawless et al. (2010) Lawless, C., D. J. Wilkinson, A. Young, S. G. Addinall, and D. A. Lydall (2010). Colonyzer: automated quantification of micro-organism growth characteristics on solid agar. BMC Bioinformatics 11, 287.
Li et al. (2011) Li, W., K. Wang, and H. Su (2011). Optimal harvesting policy for stochastic logistic population model. Applied Mathematics and Computation 218(1), 157 – 162.
Peleg et al. (2007) Peleg, M., M. G. Corradini, and M. D. Normand (2007). The logistic (Verhulst) model for sigmoid microbial growth curves revisited. Food Research International 40(7), 808 – 818.
Plummer (2010) Plummer, M. (2010). rjags: Bayesian graphical models using MCMC. R package version 2.1.0-10.
Plummer et al. (2006) Plummer, M., N. Best, K. Cowles, and K. Vines (2006). CODA: Convergence Diagnosis and Output Analysis for MCMC. R News 6(1), 7–11.
Román-Román and Torres-Ruiz (2012) Román-Román, P. and F. Torres-Ruiz (2012). Modelling logistic growth by a new diffusion process: Application to biological systems. Biosystems 110(1), 9–21.
Tsoularis and Wallace (2002) Tsoularis, A. and J. Wallace (2002). Analysis of logistic growth models. Mathematical Biosciences 179(1), 21 – 55.
Turner Jr. et al. (1976) Turner Jr., M. E., E. L. Bradley Jr., K. A. Kirk, and K. M. Pruitt (1976). A theory of growth. Mathematical Biosciences 29(3–4), 367 – 373.
Verhulst (1845) Verhulst, P. F. (1845). Recherches mathématiques sur la loi d’accroissement de la population. Nouveaux mémoires de l’Academie Royale des Science et Belles-Lettres de Bruxelles 18, 1–41.
Wallace (2010) Wallace, E. W. J. (2010). A simplified derivation of the Linear Noise Approximation. In submission,( available at, http://arxiv.org/abs/1004.4280).
West and Harrison (1997) West, M. and J. Harrison (1997). Bayesian Forecasting and Dynamic Models (Second ed.). Springer Series in Statistics. New York: Springer-Verlag.
Wilkinson (2009) Wilkinson, D. J. (2009). Stochastic modelling for quantitative description of heterogeneous biological systems. Nature Reviews Genetics 10(2), 122–133.
Wilkinson (2012) Wilkinson, D. J. (2012). Stochastic Modelling for Systems Biology. Chapman & Hall / CRC Mathematical and computational biology series. CRC Press.