Simulation and inference algorithms for stochastic biochemical reaction networks: from basic concepts to state-of-the-art

2.2.1. Mono-molecular chain

Consider two chemical species, A and B, and three reactions that form a mono-molecular chain,

with kinetic rate parameters k₁, k₂ and k₃. We adopt the convention that $\mathrm{\varnothing }$ indicates the reactions are part of an open system involving external chemical processes that are not explicitly represented in equation (2.4). Given the state vector, X(t) = [A(t), B(t)]^T, the respective propensity functions areand the stoichiometric vectors areThis mono-molecular chain is interesting since it is a part of a general class of biochemical reaction networks that are analytically tractable, though they are only applicable to relatively simple biochemical processes [36].

2.2.2. Enzyme kinetics

A biologically applicable biochemical reaction network describes the catalytic conversion of a substrate, S, into a product, P, via an enzymatic reaction involving enzyme, E. This is described by Michaelis–Menten enzyme kinetics [37,38],

with kinetic rate parameters, k₁, k₂ and k₃. This particular enzyme kinetic model is a closed system. Here, we have the state vector X(t) = [E(t), S(t), C(t), P(t)]^T, propensity functionsand stoichiometric vectorsSince the first chemical formula involves two reactant molecules, there is significantly less progress that can be made without computational methods.

See MonoMolecularChain.m and MichaelisMenten.m for example code to generate useful data structures for these biochemical reaction networks. These two biochemical reaction networks have been selected to demonstrate two stereotypical problems. In the first instance, the mono-molecular chain model, the network structure enables progress to be made analytically. The enzyme kinetic model, however, represents a more realistic case in which computational methods are required. The focus on these two representative models is done so that the exposition is clear.

3.1.1. Exact stochastic simulation algorithms

Exact SSAs generate sample paths, over some interval t ∈ [0, T], that identically follow the probability laws of the fundamental theory of stochastic chemical kinetics [33]. Take a sufficiently small time interval, [t, t + dt), such that the probability of multiple reactions occurring in this interval is zero. In such a case, the reactions are mutually exclusive events. Hence, based on equations (2.2) and (2.3), the probability of any reaction event occurring in [t, t + dt) is the sum of the individual event probabilities,

where a₀(X(t)) = a₁(X(t)) + · · · + a_M(X(t)) is the total reaction propensity function. Therefore, if we know that the next reaction occurs at time s ∈ [t, t + dt), then we can: (i) randomly select a reaction with probabilities a₁(X(s))/a₀(X(s)), …, a_M(X(s))/a₀(X(s)) and (ii) update the state vector according to the respective stoichiometric vector, ν₁, …, ν_M.

All that remains for an exact simulation method is to determine the time of the next reaction event. Gillespie [21] demonstrates that the time interval between reactions, Δt, may be treated as a random variable that is exponentially distributed with rate a₀(X(t)), that is, Δt ∼ Exp(a₀(X(t))). Therefore, we arrive at the most fundamental exact SSA, the Gillespie direct method:

(1)	initialize, time t = 0 and state vector X = x0;
(2)	calculate propensities, a1(X), …, aM(X), and a0(X) = a1(X) + · · · + aM(X);
(3)	generate random sample of the time to next reaction, Δt ∼ Exp(a0(X));
(4)	if t + Δt > T, then terminate the simulation, otherwise go to step 5;
(5)	randomly select integer j from the set {1, …, M} with $\mathbb{P}(j=1)=a_1(\mathbf{X})/a_0(\mathbf{X})$, $\dots ,\mathbb{P}(j=M)=a_M(\mathbf{X})/a_0(\mathbf{X})$;
(6)	update state, X = X + νj, and time t = t + Δt, then go to step 2.

An example implementation, GillespieDirectMethod.m, and example usage, DemoGillespie.m, are provided. Figure 1 demonstrates sample paths generated by the Gillespie direct method for the mono-molecular chain model (figure 1a) and the enzyme kinetic model (figure 1b).

• Download figure
• Open in new tab
• Download PowerPoint

A mathematically equivalent, but more computationally efficient, exact SSA formulation is derived by Gibson & Bruck [41]. Their method independently tracks the next reaction time of each reaction separately. The next reaction to occur is the one with the smallest next reaction time, therefore no random selection of reaction events is required. It should, however, be noted that the Gillespie direct method may also be improved to yield the optimized direct method [43] with similar performance benefits. Anderson [42] further refines the method of Gibson and Bruck by scaling the times of each reaction so that the scaled times between reactions follow unit-rate exponential random variables. This scaling allows the method to be applied to more complex biochemical reaction networks with time-dependent propensity functions; however, the recently proposed Extrande method [44] is computationally superior. In Anderson’s approach, scaled times are tracked for each reaction independently with t_j being the current time at the natural scale of reaction j. This results in the modified next reaction method:

(1)	initialize, global time t = 0, state vector X = x0 and scaled times t1 = t2 = · · · = tM = 0;
(2)	generate M first reaction times, s1, …, sM ∼ Exp(1);
(3)	calculate propensities, a1(X), …, aM(X);
(4)	rescale time to next reaction, Δtj = (sj − tj)/aj(X) for j = 1, 2, …, M;
(5)	choose reaction k, such that Δtk = min{Δt1, …, ΔtM};
(6)	if t + Δtk > T terminate simulation, otherwise go to step 7;
(7)	update rescaled times, tj = tj + aj(X)Δtk for j = 1, …, M, state X = X + νk and global time t = t + Δtk;
(8)	generate scaled next reaction time for reaction k, Δsk ∼ Exp(1);
(9)	update next scaled reaction time sk = sk + Δsk, and go to step 3.

An example implementation, ModifiedNextReactionMethod.m, and example usage, DemoMNRM.m, are provided. Figure 1 demonstrates sample paths generated by the modified next reaction method for the mono-molecular chain model (figure 1c) and the enzyme kinetic model (figure 1d). Note that the sample paths are different from those generated using the Gillespie direct method, despite the random number generators being initialized the same way. However, both represent exact sample paths, that is, sample paths that exactly follow the dynamics of the biochemical reaction network.

While the Gillespie direct method and the more efficient modified next reaction method and optimized direct method represent the most fundamental examples of exact SSAs, other advanced methods are also available to further improve the computational performance for large and complex biochemical reaction networks. Particular techniques include partial-propensity factorization [45], rejection-based methods [46,47] and composition-rejection [48] methods. We do not discuss these approaches, but we highlight them to indicate that efficient exact SSA method development is still an active area of research.

3.1.2. Approximate stochastic simulation algorithms

Despite some computational improvements provided by the modified next reaction method [41,42], all exact SSAs are computationally intractable for large biochemical populations and with many reactions, since every reaction event is simulated. Several approximate SSAs have been introduced in an attempt to reduce the computational burden while sacrificing accuracy.

The main approximate SSA we consider is also developed by Gillespie [39] almost 25 years after the development of the Gillespie direct method. The key idea is to evolve the system in discrete time steps of length τ, hold the propensity functions constant over the time interval [t, t + τ) and count the number of reaction events that occur. The state vector is then updated based on the net effect of all the reaction events. The number of reaction events within the interval can be shown to be a random variable distributed according to a Poisson distribution with mean a_j(X(t))τ. If Y_j denotes the number of reaction j events in [t, t + τ) then Y_j ∼ Po(a_j(X(t))τ). The result is the tau-leaping method:

(1)	initialize, time t = 0 and state Z = x0;
(2)	if t + τ > T then terminate simulation, otherwise continue;
(3)	calculate propensities, a1(Z), …, aM(Z);
(4)	generate reaction event counts, Yj ∼ Po(aj(Z)τ) for j = 1, …, M;
(5)	update state, Z = Z + Y1ν1 + · · ·YMνM, and time t = t + τ;
(6)	go to step 2.

Note that we use the notation Z(t) to denote an approximation of the true state X(t). An example implementation, TauLeapingMethod.m, and example usage, DemoTauLeap.m, are provided. Figure 1 demonstrates sample paths generated by the tau-leaping method for the mono-molecular chain model (figure 1e) and the enzyme kinetic model (figure 1f). Note that there is a visually obvious difference in the noise patterns of the tau-leaping method sample paths and the exact SSA sample paths (figure 1a–d).

The tau-leaping method is the only approximate SSA that we will explicitly discuss as it captures the essence of what approximations try to achieve; trading accuracy for improved performance. Several variations of the tau-leaping method have been proposed to improve accuracy, such as adaptive τ selection [49,50], implicit schemes [51] and the replacement of Poisson random variables with binomial random variables [52]. Hybrid methods that combine exact SSAs and approximations that split reactions into different time scales are also particularly effective for large scale networks with reactions occurring on very different time scales [53–55]. Other approximate simulation approaches are, for example, based on a continuous state chemical Langevin equation approximation [40,56,57] and employ numerical schemes for stochastic differential equations [35,58].

3.2.1. Using the chemical master equation

From a probability theory perspective, a biochemical reaction network model is a discrete-state, continuous-time Markov process. One key result for discrete-state, continuous-time Markov processes is, given an initial state, X(t) = x₀, one can describe how the probability distribution of states evolves. This is given by the chemical master equation [33],³

where P(x, t|x₀) is the probability that X(t) = x given X(0) = x₀. Solving the chemical master equation provides an explicit means of computing the probability of being in any state at any time. Unfortunately, solutions to the chemical master equation are only known for special cases [36,59].

However, the mean and variance of the biochemical reaction network molecule copy numbers can sometimes be derived without solving the chemical master equation. For example, for the mono-molecular chain model equation (2.4), one may use equation (3.1) to derive the following system of ordinary differential equations (see the electronic supplementary material):

where M_a(t) and V_a(t) (M_b(t) and V_b(t)) are the mean and variance of the copy number A(t) (B(t)) over all possible sample paths. C_a,b(t) is the covariance between A(t) and B(t). Equations (3.2)–(3.6) are linear ordinary differential equations that can be solved analytically and the solution is plotted in figure 2a superimposed with a population of sample paths.

• Download figure
• Open in new tab
• Download PowerPoint

The long time limit behaviour of a biochemical reaction network can also be determined. For the mono-molecular chain, as t → ∞ we have: M_a(t) → k₁/k₂; V_a(t) → k₁/k₂; M_b(t) → k₁/k₃; V_b(t) → k₁/k₃; and C_a,b(t) → 0. We can approximate the long time steady-state behaviour, called the stationary distribution, of the mono-molecular chain model as two independent Gaussian random variables. That is, as t → ∞, A(t) → A_∞ with $A_{\mathrm{\infty }}\sim \mathcal{N}(k_1/k_2,k_1/k_2)$. Similarly, B(t) → B_∞ with $B_{\mathrm{\infty }}\sim \mathcal{N}(k_1/k_3,k_1/k_3)$. This approximation is shown in figure 2b against histograms of sample paths generated using large values of T (see example code, DemoStatDist.m).

In the case of the mono-molecular chain model, the chemical master equation is analytically tractable [36,59]. However, the solution is algebraically complicated and non-trivial to evaluate (see the electronic supplementary material). The full chemical master equation solution in figure 2c–e and the true stationary distribution of two independent Poisson distributions is shown in figure 2f. The true stationary distribution is also compared with the Gaussian approximation in figure 2b; the approximation is reasonably accurate. However, we could have also reasonably surmised the true stationary distribution by noting that, for a Poisson distribution, the mean is equal to the variance.

3.2.2. Monte Carlo methods

The chemical master equation can yield insight for special cases; however, for practical problems, such as the enzyme kinetic model, the chemical master equation is intractable and numerical methods are required. Here, we consider numerical estimation of the mean state vector at a fixed time, T.

In probability theory, the mean of a distribution is defined via an expectation,

where Ω is the set of all possible states. It is important to note that the methods we describe here are equally valid for a more general expectations of the form $\mathbb{E}[f(\mathbf{X}(T))]$ where f is some function that satisfies certain conditions [60].

We usually cannot compute equation (3.7) directly since Ω is typically infinite and the chemical master equation is intractable. However, exact SSAs provide methods for sampling the chemical master equation distribution, X(T) ∼ P(x, T|x₀). This leads to the Monte Carlo estimator

where X(T)⁽¹⁾, …, X(T)⁽ⁿ⁾ are n independent sample paths of the biochemical reaction network of interest (see example implementation MonteCarlo.m).

Unlike equation (3.7), the Monte Carlo estimates, such as equation (3.8), are random variables for finite n. This incurs a probabilistic error. A common measure of the accuracy of a Monte Carlo estimator, $\hat{\mu }(T)$, of some expectation, $\mathbb{E}[\mu (T)]$, is the mean-square error that evaluates the average error behaviour and may be decomposed as follows:

Equation (3.9) highlights that there are two sources of error in a Monte Carlo estimator, the estimator variance and bias, and much of the discussion that follows deals with how to balance both these error sources in a computationally efficient manner.

Through analysis of the mean-square error of an estimator, the rate at which the mean-square error decays as n increases can be determined. Hence, we can determine how large n needs to be to satisfy the condition

where c is a positive constant and h is called the error tolerance.

Since $\mathbb{E}[\hat{\mathbf{X}}(T)]=\mathbb{E}[\mathbf{X}(T)]$, the bias term in equation (3.9) is zero and we call $\hat{\mathbf{X}}(T)$ an unbiased estimator of $\mathbb{E}[\mathbf{X}(T)]$. For an unbiased estimator, the mean-square error is equal to the estimator variance. Furthermore, $\text{Var}[\hat{\mathbf{X}}(T)]=\text{Var}[\mathbf{X}(T)]/n$, so the estimator variance decreases linearly with n, for sufficiently large n. Therefore, $h\propto 1/\sqrt{n}$. That is, to halve h, one must increase n by a factor of four. This may be prohibitive with exact SSAs, especially for biochemical reaction networks with large variance. In the context of the Monte Carlo estimator using an exact SSA, equation (3.8), for n sufficiently large, the central limit theorem (CLT) states that $\hat{\mathbf{X}}(T)\sim \mathcal{N}(\mathbb{E}[\mathbf{X}(T)],\text{Var}[\mathbf{X}(T)]/n)$ (see Wilkinson [22] for a good discussion on the CLT).

Computational improvements can be achieved by using an approximate SSA, such as the tau-leaping method,

where Z(T)⁽¹⁾, …, Z(T)⁽ⁿ⁾ are n independent approximate sample paths of the biochemical reaction network of interest (see the example implementation MonteCarloTauLeap.m).

Note that $\mathbb{E}[\hat{\mathbf{Z}}(T)]=\mathbb{E}[\mathbf{Z}(T)]$. Since $\mathbb{E}[\mathbf{Z}(T)]\ne \mathbb{E}[\mathbf{X}(T)]$ in general, we call $\hat{\mathbf{Z}}(T)$ a biased estimator. Even in the limit of n → ∞, the bias term in equation (3.9) may not be zero, which incurs a lower bound on the best achievable error tolerance, h, for fixed τ. However, it has been shown that the bias of the tau-leaping method decays linearly with τ [61,62]. Therefore, to satisfy the error tolerance condition (equation (3.10)) we not only require n ∝ 1/h² but also τ ∝ h. That is, as h decreases, the performance improvement of Monte Carlo with the tau-leaping method reduces by a factor of τ, because the computational cost of the tau-leaping method is proportional to 1/τ. In figure 3a, the decay of the computational advantage in using the tau-leaping method for Monte Carlo over the Gillespie direct method is evident, and eventually the tau-leaping method will be more computationally burdensome than the Gillespie direct method. By the CLT, for large n, we have $\hat{\mathbf{Z}}(T)\sim \mathcal{N}(\mathbb{E}[\mathbf{Z}(T)],\text{Var}[\mathbf{Z}(T)]/n)$.

• Download figure
• Open in new tab
• Download PowerPoint

The utility of the tau-leaping method for accurate (or exact) Monte Carlo estimation is identified by Anderson & Higham [63] through extending the idea of multilevel Monte Carlo (MLMC) originally proposed by Giles for stochastic differential equations [60,64]. Consider a sequence of L + 1 tau-leaping method time steps τ₀, τ₁, …, τ_L, with τ_ℓ < τ_ℓ−1 for ℓ = 1, …, L. Let Z_ℓ(T) denote the state vector of a tau-leaping method approximation using time step τ_ℓ. Assume τ_L is small enough that $\mathbb{E}[{\mathbf{Z}}_L(T)]$ is a good approximation of $\mathbb{E}[\mathbf{X}(T)]$. Note that for large τ_ℓ (small ℓ), sample paths are cheap to generate, but inaccurate; conversely, small τ_ℓ (large ℓ) results in computationally expensive, but accurate sample paths.

We can write

Importantly, the final estimator on the right of equation (3.12), called the multilevel telescoping summation, is equivalent in bias to $\mathbb{E}[{\mathbf{Z}}_L(T)]$. At first glance, equation (3.12) looks to have complicated the computational problem and inevitably decreased performance of the Monte Carlo estimator. The insight of Giles [64], in the context of stochastic differential equation models for financial applications, is that the bias correction terms may be computed using Monte Carlo approaches that involve generating highly correlated sample paths in estimation of each of the correction terms, thus reducing the variance in the bias corrections. If the correlation is strong enough, then the variance decays such that few of the most accurate tau-leaping method sample paths are required; this can result in significant computational savings.

A contribution of Anderson & Higham [63] is an efficient method of generating correlated tau-leaping method sample path pairs (Z_ℓ(T), Z_ℓ−1(T)) in the case when τ_ℓ = τ_ℓ−1/δ for some positive integer scale factor δ. The algorithm is based on the property that the sum of two Poisson random variables is also a Poisson random variable. This enables the sample path with τ_ℓ−1 to be constructed as an approximation to the sample path with τ_ℓ directly. Figure 3b demonstrates tau-leaping method sample paths of B(t) in the mono-molecular chain model with τ = 2, 4, 8 generated directly from a tau-leaping method sample path with τ = 1. The algorithm can be thought of as generating multiple approximations of the same exact biochemical reaction network sample path. The algorithm is the correlated tau-leaping method:

(1)	initialize time t = 0, and states Zℓ, Zℓ−1 corresponding to sample paths with τ = τℓ and τ = τℓ−1, respectively;
(2)	if t + τℓ > T, then terminate simulation, otherwise continue;
(3)	calculate propensities for path Zℓ, a1(Zℓ), …, aM(Zℓ);
(4)	if t/τℓ is not an integer multiple of δ, then go to step 6, otherwise continue;
(5)	calculate propensities for path Zℓ−1, a1(Zℓ−1), …, aM(Zℓ−1), initialize intermediate state $\overline{\mathbf{Z}}={\mathbf{Z}}_{\ell -1}$;
(6)	for each reaction j = 1, …, M; 6.1 calculate virtual propensities, bj,1 = min{aj(Zℓ), aj(Zℓ−1)}, bj,2 = aj(Zℓ) − bj,1 and bj,3 = aj(Zℓ−1) − bj,1; 6.2 generate virtual reaction event counts, Yj,1 ∼ Po(bj,1τℓ), Yj,2 ∼ Po(bj,2τℓ) and Yj,3 ∼ Po(bj,3τℓ);
(1)	set, Zℓ = Zℓ + (Y1,1 + Y1,2)ν1 + · · · + (YM,1 + YM,2)νM;
(1)	set, $\overline{\mathbf{Z}}=\overline{\mathbf{Z}}+(Y_{1,1}+Y_{1,3}){\mathit{\nu }}_1+\cdots +(Y_{M,1}+Y_{M,3}){\mathit{\nu }}_M$;
(1)	if (t + τℓ)/τℓ is an integer multiple of δ, then set ${\mathbf{Z}}_{\ell -1}=\overline{\mathbf{Z}}$;
(1)	update time t = t + τℓ, and go to step 2.

See example implementation, CorTauLeapingMethod.m, and example usage, DemoCorTauLeap.m.

Given the correlated tau-leaping method, Monte Carlo estimation can be applied to each term in equation (3.12) to give

where ${\mathbf{Z}}_0{(T)}^{(1)},\dots ,{\mathbf{Z}}_0{(T)}^{(n_0)}$ are n₀ independent tau-leaping method sample paths with τ = τ₀, and (Z_ℓ(T)⁽¹⁾, ${\mathbf{Z}}_{\ell -1}{(T)}^{(1)}),\dots ,({\mathbf{Z}}_{\ell }{(T)}^{(n_{\ell })},{\mathbf{Z}}_{\ell -1}{(T)}^{(n_{\ell })})$ are n_ℓ paired correlated tau-leaping method sample paths with time steps τ = τ_ℓ, τ = τ_ℓ−1 and τ_ℓ−1 = δτ_ℓ for each bias correction ℓ = 1, 2, …, L. Given an error tolerance, h, it is possible to calculate an optimal sequence of sample path numbers n₀, n₁, …, n_L such that the total computation time is optimized [63–65]. The results are shown in figure 3a for a more computationally challenging parameter set of the mono-molecular chain model. See the example implementation, MultilevelMonteCarlo.m and DemoMonteCarlo.m, for the full performance comparison.

As formulated here, MLMC results in a biased estimator, though it is significantly more efficient to reduce the bias of this estimator than by direct use of the tau-leaping method. If an unbiased estimator is required, then this can be achieved by correlating exact SSA sample paths with approximate SSA sample paths. Anderson & Higham [63] demonstrate a method for correlating tau-leaping method sample paths and modified next reaction method sample paths, and Lester et al. [65] demonstrate correlating tau-leaping method sample paths and Gillespie direct method sample paths. Further refinements such as adaptive and non-nested τ_ℓ steps are also considered by Lester et al. [66], a multilevel hybrid scheme is developed by Moraes et al. [67] and Wilson & Baker [68] use MLMC and maximum entropy methods to generate approximations to the chemical master equation.

4.1.1. Example data

The computation examples given in this review are based on two synthetically generated data sets, corresponding to the biochemical reaction network models given in §2. This enables the comparison between inference methods and the accuracy of inference.

The data for inference on the mono-molecular chain model (equation (2.4)) are taken as perfect observations, that is, K = N, A = I and $\mathbb{P}(\mathit{\xi }=\mathbf{0})=1$. A single sample path is generated for the mono-molecular chain model with true parameters, θ_true = [1.0, 0.1, 0.05], and initial condition A(0) = 100 and B(0) = 0 using the Gillespie direct method. Observations are made using equation (4.1) applied at n_t = 4 discrete times, t₁ = 25, t₂ = 50, t₃ = 75 and t₄ = 100. The resulting data are given in the electronic supplementary material.

The data for inference on the enzyme kinetic model equation (2.7) assumes incomplete and noisy observations. Specifically, only the product is observed, so K = 1, A = [0, 0, 0, 1]. Furthermore, we assume that there is some measurement error, $\xi \sim \mathcal{N}(0,4)$; that is, the error standard deviation is two product molecules. The data are generated using the Gillespie direct method with true parameters θ_true = [0.001, 0.005, 0.01] and initial condition E(0) = 100, S(0) = 100, C(0) = 0 and P(0) = 0. Equation (4.1) is evaluated at n_t = 5 discrete times, t₁ = 0, t₂ = 20, t₃ = 40, t₄ = 60 and t₅ = 80 (including an observation of the initial state), yielding the data in the electronic supplementary material.

4.3.1. Exact Bayesian sampling

Assuming the likelihood can be evaluated, that is, the chemical master equation is tractable, a naive method of generating m samples from the posterior is the exact rejection sampler:

(1)	initialize index i = 0;
(2)	generate a prior sample θ* ∼ p(θ);
(3)	calculate acceptance probability α = p(Yobs|θ*);
(4)	with probability α, accept θ(i+1) = θ* and i = i + 1;
(5)	if i = m, terminate sampling, otherwise go to step 2.

Unsurprisingly, this approach is almost never viable as the likelihood probabilities are often extremely small. In the code example, DemoExactBayesRejection.m, the acceptance probability is never more than 9 × 10⁻¹⁵.

The most common solution is to construct a type of stochastic process (a Markov chain) in the parameter space to generate m_n steps, ${\mathit{\theta }}^{(0)},\dots ,{\mathit{\theta }}^{(m_n)}$. An essential property of the Markov chain used is that its stationary distribution is the target posterior distribution. This approach is called Markov chain Monte Carlo (MCMC), and a common method is the Metropolis–Hastings method [95,96]:

(1)	initialize n = 0 and select starting point θ(0);
(2)	generate a proposal sample, θ* ∼ q(θ|θ(n));
(3)	calculate acceptance probability $$\alpha =min(1,\frac{\hspace{0.17em}p({\mathbf{Y}}_{\mathrm{obs}}\mid {\mathit{\theta }}^{\ast })p({\mathit{\theta }}^{\ast })q({\mathit{\theta }}^{(n)}\mid {\mathit{\theta }}^{\ast })}{p({\mathbf{Y}}_{\mathrm{obs}}\mid {\mathit{\theta }}^{(n)})p({\mathit{\theta }}^{(n)})q({\mathit{\theta }}^{\ast }\mid {\mathit{\theta }}^{(n)})});$$
(4)	with probability α, set θ(n+1) = θ*, and with probability 1 − α, set θ(n+1) = θ(n);
(5)	update time, n = n + 1;
(6)	if n > mn, terminate simulation, otherwise go to step 2.

The acceptance probability is based on the relative likelihood between two parameter configurations, the current configuration θ⁽ⁿ⁾ and a proposed new configuration θ*, as generated by the user-defined proposal kernel q(θ|θ⁽ⁿ⁾). It is essential to understand that since the samples, ${\mathit{\theta }}^{(0)},\dots ,{\mathit{\theta }}^{(m_n)}$, are produced by a Markov chain we cannot treat the m_n steps as m_n independent posterior samples. Rather, we need to take m_n to be large enough that the m_n steps are effectively equivalent to m independent samples.

One challenge in MCMC sampling is the selection of a proposal kernel such that the size of m_n required for the Markov chain to reach the stationary distribution, called the mixing time, is small [97]. If the variance of the proposal is too low, then the acceptance rate is high, but the mixing is poor since only small steps are ever taken. On the other hand, a proposal variance that is too high will almost always make proposals that are rejected, resulting in many wasted likelihood evaluations before a successful move event. Selecting good proposal kernels is a non-trivial exercise and we refer the reader to Cotter et al. [98], Green et al. [99] and Roberts & Rosenthal [100] for detailed discussions on the wide variety of MCMC samplers including proposal techniques. Other techniques used to reduce correlations in MCMC samples include discarding the first m_b steps, called burn-in iterations, or sub-sampling the chain by using only every m_hth step, also called thinning. However, in general, the use of thinning decreases the statistical efficiency of the MCMC estimator [101,102].

Alternative approaches for exact Bayesian sampling can also be based on importance sampling [22,103]. Consider a random variable that cannot be simulated, X ∼ p(x), but suppose that it is possible to simulate another random variable, Y ∼ q(y). If X, Y ∈ Ω and p(x) = 0 whenever q(y) = 0, then

where Y⁽¹⁾, …, Y^(m) are independent samples of q(Y). Using equation (4.8), one can show that if the distributions of the target, p(x), and the proposal, q(y), are similar, then a collection of samples, Y⁽¹⁾, …, Y^(m), can be used to generate m approximate samples from p(x). This is called importance resampling:

(1)	generate samples Y(1), …, Y(m) ∼ q(y);
(2)	compute weights w(1) = p(Y(1))/q(Y(1)), …, w(m) = p(Y(m))/q(Y(m)); for i = 1, 2, …, m;
(3)	generate samples {X(1), …, X(m)} by drawing from {Y(1), …, Y(m)} with replacement using probabilities $\mathbb{P}(X=Y^{(i)})=w^{(i)}/\sum _{i=1}^m{w}^{(i)}$.

In Bayesian applications, the prior is often very different from the posterior. In such a case, importance resampling may be applied using a sequence of intermediate distributions. This is called sequential importance resampling and is one approach from the family of sequential Monte Carlo (SMC) [104] samplers. However, like the Metropolis–Hastings method, these methods also require explicit calculations of the likelihood function in order to compute the weights, thus all of these approaches are infeasible for practical biochemical reaction networks. Therefore, we only present more practical forms of these methods later.

More recently, it has been shown that the MLMC telescoping summation can accelerate the computation of posterior expectations when using MCMC [105,106] or SMC [107]. The key challenges in these applications is the development of appropriate coupling strategies. We do not cover these technical details in this review.

4.3.2. Likelihood-free methods

Since exact Bayesian sampling is rarely tractable, due to the intractability of the chemical master equation, alternative, likelihood-free sampling methods are required. Two main classes of likelihood-free methods exist: (i) so-called pseudo-marginal MCMC and (ii) approximate Bayesian computation (ABC). The focus of this review is ABC methods, though we first briefly discuss the pseudo-marginal approach.

The basis of the pseudo-marginal approach is to use a MCMC sampler (e.g. Metropolis–Hastings method), but replace the explicit likelihood evaluation with a likelihood estimate obtained through Monte Carlo simulation of the forwards problem [108]. For example, a direct unbiased Monte Carlo estimator is

where ${[\mathbf{X}{(t_1)}^{(j)},\dots ,\mathbf{X}{(t_{n_t})}^{(j)}]}^T$ for j = 1, 2, …, n are independent sample paths of the biochemical reaction network of interest observed discretely at times $t_1,t_2,\dots ,t_{n_t}$. The most successful class of pseudo-marginal techniques, particle MCMC [109], apply a SMC sampler to approximate the likelihood and inform the MCMC proposal mechanism. The particle marginal Metropolis–Hastings method is a popular variant [30,109]. However, the recent model based proposals variant also holds promise for biochemical reaction networks specifically [110].

A particularly nice feature of pseudo-marginal methods is that they are unbiased,⁶ that is, the Markov chain will still converge to the exact target posterior distribution. This property is sometimes called an ‘exact approximation’ [22,30]. Unfortunately, the Markov chains in these methods typically converge more slowly than their exact counterparts. However, computational improvements have been obtained through application of MLMC [111].

Another popular likelihood-free approach is ABC [73,112–114]. ABC methods have enabled very complex models to be investigated, that are otherwise intractable [32]. Furthermore, ABC methods are very intuitive, leading to wide adoption within the scientific community [32]. Applications are particularly prevalent in the life sciences, especially in evolutionary biology [112,113,115–117], cell biology [118–120], epidemiology [121], ecology [122,123] and systems biology [31,90].

The basis of ABC is a discrepancy metric ρ(Y_obs, S_obs) that provides a measure of closeness between the data, Y_obs, and simulated data S_obs generated through stochastic simulation of the biochemical reaction network with simulated measurement error. Thus, acceptance probabilities are replaced with a deterministic acceptance condition, ρ(Y_obs, S_obs) ≤ ε, where ε is the discrepancy threshold. This yields an approximation to Bayes’ theorem,

The key insight here is that the ability to produce sample paths of a biochemical reaction network model, that is, the forwards problem, enables an approximate algorithm for inference, that is, the inverse problem, regardless of the intractability or complexity of the likelihood. In fact, a formula for the likelihood need not even be known.

The discrepancy threshold determines the level of approximation; however, under the assumption of model and observation error, equation (4.9) can be treated as exact [124]. As ε → 0 then p(θ|ρ(Y_obs, S_obs) ≤ ε) → p(θ|Y_obs) [125,126]. Using data for the mono-molecular chain model, we can demonstrate this convergence, as shown in figure 4. The marginal posterior distributions are shown for each parameter, for various values of ε and compared with the exact marginal posteriors. The discrepancy metric used is

where S(t) is simulated data generated using the Gillespie direct method and equation (4.1). The example code, DemoABCConvergence.m, is used to generate these marginal distributions.

• Download figure
• Open in new tab
• Download PowerPoint

For any ε > 0, ABC methods are biased, just as the tau-leaping method is biased for the forwards problem. Therefore, a Monte Carlo estimate of a posterior summary statistic, such as equation (4.7), needs to take this bias into account. Just like the tau-leaping method, the rate of convergence in mean-square of ABC-based Monte Carlo is degraded because of this bias [125,126]. Furthermore, as the dimensionality of the data increases, small values of ɛ are not computationally feasible. In such cases, the data dimensionality may be reduced by computing lower dimensional summary statistics [32]; however, these must be sufficient statistics in order to ensure the ABC posterior still converges to the correct posterior as ε → 0. We refer the reader to Fearnhead & Prangle [126] for more detail on this topic.

4.3.3. Samplers for approximate Bayesian computation

We now focus on computational methods for generating m samples, ${\mathit{\theta }}_{ϵ}^{(1)},\dots ,{\mathit{\theta }}_{ϵ}^{(m)}$, from the ABC posterior equation (4.9) with ε > 0 and discrepancy metric as given in equation (4.10). Throughout, we denote s(S_obs; θ) as the process for generating simulated data given a parameter vector; this process is identical to the processes used to generate our synthetic example data.

In general, the ABC samplers are only computationally viable if the data simulation process is not computationally expensive, in the sense that it is feasible to generate millions of sample paths. However, this is not always realistic, and many extensions exist to standard ABC in an attempt to deal with these more challenging cases. The lazy ABC method [127] applies an additional probability rule that terminates simulations early if it is likely that ρ(Y_obs, S_obs) > ε. The approximate ABC method [128,129] uses a small set of data simulations to construct an approximation to the data simulation process.

The most notable early applications of ABC samplers are Beaumont et al. [115], Pritchard et al. [112] and Taveré et al. [113]. The essential development of this early work is that of the ABC rejection sampler:

(1)	initialize index i = 0;
(2)	generate a prior sample θ* ∼ p(θ);
(3)	generate simulated data, ${\mathbf{S}}_{\mathrm{obs}}^{\ast }\sim s({\mathbf{S}}_{\mathrm{obs}};{\mathit{\theta }}^{\ast })$;
(4)	if $\rho ({\mathbf{Y}}_{\mathrm{obs}},{\mathbf{S}}_{\mathrm{obs}}^{\ast })\le ϵ$, accept ${\mathit{\theta }}_{ϵ}^{(i+1)}={\mathit{\theta }}^{\ast }$ and set i = i + 1, otherwise continue;
(5)	if i = m, terminate sampling, otherwise go to step 2.

There is a clear connection with the exact rejection sampler. Note that every accepted sample of the ABC posterior corresponds to at least one simulation of the forwards problem as shown in figure 5. As a result, the computational burden of the inverse problem is significantly higher than the forwards problem, especially for small ε. The example code, ABCRejectionSampler.m, provides an implementation of the ABC rejection sampler.

Unfortunately, for small ε, the computational burden of the ABC rejection sampler may be prohibitive as the acceptance rate is very low (this is especially an issue for biochemical reaction networks with highly variable dynamics). For example in figure 4, m = 100 posterior samples takes approximately one minute for ε = 25, but nearly ten hours for ε = 12.5. However, for ε = 12.5 the marginal ABC posterior for k₃ has not yet converged.

• Download figure
• Open in new tab
• Download PowerPoint

Marjoram et al. [130] provide a solution via an ABC modification to the Metropolis–Hastings method. This results in a Markov chain with the ABC posterior as the stationary distribution. This method is called ABC Markov chain Monte Carlo (ABCMCMC):

(1)	initialize n = 0 and select starting point ${\mathit{\theta }}_{ϵ}^{(0)}$;
(2)	generate a proposal sample, ${\mathit{\theta }}^{\ast }\sim q(\mathit{\theta }\mid {\mathit{\theta }}_{ϵ}^{(n)})$;
(3)	generate simulated data, ${\mathbf{S}}_{\mathrm{obs}}^{\ast }\sim s({\mathbf{S}}_{\mathrm{obs}};{\mathit{\theta }}^{\ast })$;
(4)	if $\rho ({\mathbf{Y}}_{\mathrm{obs}},{\mathbf{S}}_{\mathrm{obs}}^{\ast })>ϵ$, then set ${\mathit{\theta }}_{ϵ}^{(n+1)}={\mathit{\theta }}_{ϵ}^{(n)}$ and go to step 7, otherwise continue;
(5)	calculate acceptance probability $$\alpha =min(1,\frac{\hspace{0.17em}p({\mathit{\theta }}^{\ast })q({\mathit{\theta }}_{ϵ}^{(n)}\mid {\mathit{\theta }}^{\ast })}{p({\mathit{\theta }}_{ϵ}^{(n)})q({\mathit{\theta }}^{\ast }\mid {\mathit{\theta }}_{ϵ}^{(n)})});$$
(6)	with probability α, set ${\mathit{\theta }}_{ϵ}^{(n+1)}={\mathit{\theta }}^{\ast }$, and with probability 1 − α, set ${\mathit{\theta }}_{ϵ}^{(n+1)}={\mathit{\theta }}_{ϵ}^{(n)}$;
(7)	update time, n = n + 1;
(8)	if n > mn, terminate simulation, otherwise go to step 2.

An example implementation, ABCMCMCSampler.m, is provided.

Just as with the Metropolis–Hastings method, the efficacy of the ABCMCMC rests upon the non-trivial choice of the proposal kernel, $q(\mathit{\theta }\mid {\mathit{\theta }}_{ϵ}^{(n)})$. The challenge of constructing effective proposal kernels is equally non-trivial for ABCMCMC as for the Metropolis–Hastings method. Figure 6 highlights different Markov chain behaviours for heuristically chosen proposal kernels based on Gaussian random walks with variances that we alter until the Markov chain seems to be mixing well.

• Download figure
• Open in new tab
• Download PowerPoint

For the mono-molecular chain model (figure 6a–c), ABCMCMC seems to be performing well with the Markov chain state moving largely in the regions of high posterior probability. However, for the enzyme kinetics model (figure 6d–f), the Markov chain has taken a long excursion into a region of low posterior probability (figure 6f) where it is a long way from the true parameter value of k₃ = 0.01 for many steps. Such an extended path into the low probability region results in significant bias in this region and many more steps of the algorithm are required to compensate for this [131,132]. Just as with exact MCMC, burn-in samples and thinning may be used to reduce correlations but may not improve the final Monte Carlo estimate.

To avoid the difficulties in ensuring efficient ABCMCMC convergence, Sisson et al. [132] developed an ABC variant of SMC. Alternative versions of the method are also designed by Beaumont et al. [131] and Toni et al. [31]. The fundamental idea is to use sequential importance resampling to propagate m_p samples, called particles, through a sequence of R + 1 ABC posterior distributions defined through a sequence of discrepancy thresholds, ε₀, ε₁, …, ε_R with ε_r > ε_r+1 for r = 0, 1, …, R − 1 and ε₀ = ∞ (that is, $p({\mathit{\theta }}_{{ϵ}_0}\mid {\mathbf{Y}}_{\mathrm{obs}})$ is the prior). This results in ABCSMC:

(1)	initialize r = 0 and weights $w_r^{(i)}=1/m_{\hspace{0.17em}p}$ for i = 1, …, mp;
(2)	generate mp particles from the prior, ${\mathit{\theta }}_{{ϵ}_r}^{(i)}\sim p(\mathit{\theta })$, for i = 1, 2, …, mp;
(3)	set index i = 0;
(4)	randomly select integer, j, from the set {1, …, mp} with $\mathbb{P}(\hspace{0.17em}j=1)=w_r^{(1)},\dots ,\mathbb{P}(\hspace{0.17em}j=m_{\hspace{0.17em}p})=w_r^{(m_{\hspace{0.17em}p})}$;
(5)	generate proposal, ${\mathit{\theta }}^{\ast }\sim q(\mathit{\theta }\mid {\mathit{\theta }}_{{ϵ}_r}^{(j)})$;
(6)	generate simulated data, ${\mathbf{S}}_{\mathrm{obs}}^{\ast }\sim s({\mathbf{S}}_{\mathrm{obs}};{\mathit{\theta }}^{\ast })$;
(7)	if $\rho ({\mathbf{Y}}_{\mathrm{obs}},{\mathbf{S}}_{\mathrm{obs}}^{\ast })>{ϵ}_{r+1}$, go to step 4, otherwise continue;
(8)	set ${\mathit{\theta }}_{{ϵ}_{r+1}}^{(i+1)\ast }={\mathit{\theta }}^{\ast }$, $w_{r+1}^{(i+1)\ast }=p({\mathit{\theta }}^{\ast })/[\sum _{\hspace{0.17em}j=1}^{m_{\hspace{0.17em}p}}{w}_r^{(j)}q({\mathit{\theta }}^{\ast }\mid {\mathit{\theta }}_{{ϵ}_r}^{(j)})]$ and i = i + 1;
(9)	if i < mp, go to step 4, otherwise continue;
(10)	set index i = 0;
(11)	randomly select integer, j, from the set {1, …, mp} with $\mathbb{P}(\hspace{0.17em}j=1)=w_{r+1}^{(1)\ast },\dots ,\mathbb{P}(j=m_{\hspace{0.17em}p})=w_{r+1}^{(m_{\hspace{0.17em}p})\ast }$;
(12)	set ${\mathit{\theta }}_{r+1}^{i+1}={\mathit{\theta }}_{r+1}^{(j)\ast }$ and i = i + 1;
(13)	if i < mp, go to step 11, otherwise continue;
(14)	set $w_{r+1}^{(i)}=1/[\sum _{\hspace{0.17em}j=1}^{m_{\hspace{0.17em}p}}{w}_{r+1}^{(j)\ast }]$ and r = r + 1;
(15)	if r < R, go to step 3, otherwise terminate simulation.

An implementation of the ABCSMC method is provided in ABCSMCSampler.m.

The ABCSMC method avoids some issues inherent in ABCMCMC such as strong sample correlation and long excursions into low probability regions. However, it is still the case that the number of particles, m_p, often needs to be larger than the desired number of independent samples, m, from the ABC posterior with ε_R. This is especially true when there is a larger discrepancy between two successive ABC posteriors in the sequence. Although techniques exist to adaptively generate the sequence of acceptance thresholds [133]. Also note that ABCSMC still requires a proposal kernel to mutate the particles at each step and the efficiency of the method is affected by the choice of proposal kernel. Just as with ABCMCMC, the task of selecting an optimal proposal kernel is non-trivial [134].

The formulation of ABCSMC using a sequence of discrepancy thresholds hints that MLMC ideas could also be applicable. Recently, a variety of MLMC methods for ABC have been proposed [135–137]. All of these approaches are similar in their application of the multilevel telescoping summation to compute expectations with respect to ABC posterior distributions,

where ε₀, ε₁, …, ε_L is a sequence of discrepancy thresholds with ε_ℓ > ε_ℓ+1 for all ℓ = 0, 1, …, L − 1.

Unlike the forwards problem, no exact solution has been found for the generation of correlated ABC posterior pairs, $({\mathit{\theta }}_{{ϵ}_{\ell }},{\mathit{\theta }}_{{ϵ}_{\ell -1}})$, with ℓ = 1, …, L. Several approaches to this problem have been proposed: Guha & Tan [135] use a sequence of correlated ABCMCMC samplers in a similar manner to Dodwell et al. [105] and Efendiev et al. [106]; Jasra et al. [136] apply the MLSMC estimator of Beskos et al. [107] directly to equation (4.11); and Warne et al. [137] develop a MLMC variant of ABC rejection sampler through the introduction of an approximation based on the empirical marginal posterior distributions.

The coupling scheme of Warne et al. [137] is the simplest to present succinctly. The scheme is based on the inverse transform method for sampling univariate distributions. Given random variable X ∼ p(x), the cumulative distribution function is defined as

Note that 0 ≤ F(z) ≤ 1. If F(z) has an inverse function, F⁻¹(u), and U⁽¹⁾, U⁽²⁾, …, U^(m) are independent samples of $\mathcal{U}(0,1)$, then independent samples from p(x) can be generated using X⁽¹⁾ = F⁻¹(U⁽¹⁾), X⁽²⁾ = F⁻¹(U⁽²⁾), …, X^(m) = F⁻¹(U^(m)).

Given samples, ${\mathit{\theta }}_{{ϵ}_{\ell }}^1,{\mathit{\theta }}_{{ϵ}_{\ell }}^2,\dots ,\hspace{0.17em}{\mathit{\theta }}_{{ϵ}_{\ell }}^{m_{\ell }}$, an estimate for the posterior marginal cumulative distribution functions can be obtained using

where $k_{{ϵ}_{\ell },j}^{(i)}$ is the jth component of ${\mathit{\theta }}_{{ϵ}_{\ell }}^{(i)}$ for all j = 1, …, M, and ${\mathbb{1}}_z(x)=1$ if x ≤ z and ${\mathbb{1}}_z(x)=0$ otherwise. Assuming we already have an estimate of the posterior marginal cumulative distribution functions for acceptance threshold ε_ℓ−1, then we apply the transform $k_{{ϵ}_{\ell -1},j}^{(i)}={\hat{F}}_{\ell -1,j}^{-1}({\hat{F}}_{\ell ,j}(k_{{ϵ}_{\ell },j}^{(i)}))$. This results in a coupled marginal pair ${(k_{{ϵ}_{\ell },j},k_{{ϵ}_{\ell -1},j})}^{(i)}$ since it is equivalent to using the inverse transform method to generate $k_{{ϵ}_{\ell },j}^{(i)}$ and $k_{{ϵ}_{\ell -1},j}^{(i)}$ using the same uniform sample $U^{(i)}\sim \mathcal{U}(0,1)$. These approximately correlated samples can be obtained using:

(1)	generate sample from level ℓ ABC posterior ${\mathit{\theta }}_{{ϵ}_{\ell }}\sim p({\mathit{\theta }}_{{ϵ}_{\ell }}\mid {\mathbf{Y}}_{\mathrm{obs}})$ using the ABC rejection sampler;
(2)	set $k_{{ϵ}_{\ell -1},j}={\hat{F}}_{\ell -1,j}^{-1}({\hat{F}}_{\ell ,j}(k_{{ϵ}_{\ell },j}))$ for j = 1, 2, …, M;
(3)	set ${\mathit{\theta }}_{{ϵ}_{\ell -1}}=[k_{{ϵ}_{\ell -1},1},k_{{ϵ}_{\ell -1},2},\dots ,\hspace{0.17em}{k}_{{ϵ}_{\ell -1},M}]$.

Thus, the ABC multilevel Monte Carlo (ABCMLMC) [137] method proceeds though computing the Monte Carlo estimate to equation (4.11),

where ${\mathit{\theta }}_{{ϵ}_0}^{(1)},\dots ,\hspace{0.17em}{\mathit{\theta }}_{{ϵ}_0}^{(m_0)}$ are m₀ independent samples using ABC rejection sampler and $({\mathit{\theta }}_{{ϵ}_{\ell }}^{(1)},{\mathit{\theta }}_{{ϵ}_{\ell -1}}^{(1)})$, $\dots ,({\mathit{\theta }}_{{ϵ}_{\ell }}^{(m_{\ell })},{\mathit{\theta }}_{{ϵ}_{\ell -1}}^{(m_{\ell })})$ are m_ℓ independent sample pairs using the approximately correlated ABC rejection sampler for ℓ = 1, 2, …, L. It is essential that each bias correction term in equation (4.12) be estimated in ascending order to ensure the marginal CDF estimates for ε_ℓ−1 are available for use in the coupling scheme for the ℓth bias correction term. The complete algorithm is provided in the electronic supplementary material and the example code is provided in ABCMLMC.m. We refer the reader to Warne et al. [137] for further details.

Provided successive ABC posteriors are similar in correlation structure, then the ABCMLMC method accelerates ABC rejection sampler through reduction of the number of m_L samples required for a given error tolerance. However, additional bias is introduced by the approximately correlated ABC rejection sampler since the telescoping summation is not strictly satisfied. In practice, this affects the choice of the acceptance threshold sequence [137].

We provide a comparison of the ABC rejection sampler, ABCMCMC, ABCSMC and ABCMLMC methods as presented here. Posterior means are computed for both the mono-molecular model and enzyme kinetics model inference problems along with the $95\text{\%}$ CIs for each estimate⁷ (see electronic supplementary material). Computation times and parameter estimates are provided in tables 1 and 2. Example code, DemoABCMethodsMonoMol.m and DemoABCMethodsMichMent.m, is provided. Algorithm configurations are supplied in the electronic supplementary material.

Note that no extensive tuning of the ABCMCMC, ABCSMC or ABCMLMC algorithms is performed, thus these outcomes do not reflect a detailed and optimized implementation of the methods. Such a benchmark would be significantly more computationally intensive than the comparison of Monte Carlo methods for the forwards problem (figure 3). The computational statistics literature demonstrates that ABCMCMC [130], ABCSMC [131,132] and ABCMLMC [135–137] can be tuned to provide very competitive results on a given inference problem. However, a large number of trial computations are often required to achieve this tuning, or more complex adaptive schemes need to be exploited [100,133]. Instead, our comparisons represent a more practical guide in which computations are kept short and almost no tuning is performed, thus giving a fair comparison for a fixed, short computational budget.

For both inference problems, ABCSMC and ABCMLMC produce more accurate parameter estimates in less computation time than ABC rejection sampler and ABCMCMC. ABCSMC performs better on the mono-molecular chain model and ABCMLMC performs better on the enzyme kinetic model. This is not to say that ABCMCMC cannot be tuned to perform more efficiently, but it does indicate that ABCMCMC is harder to tune without more extensive testing. An advantage of ABCMLMC is that it has fewer components that require tuning, since the sample numbers may be optimally chosen, and the user need only provide a sequence of acceptance thresholds [137]. However, the ability to tune the proposal kernel in ABCSMC can be a significant advantage, despite the challenge of determining a good choice for such a proposal [133,134].