Approximate Bayesian computation scheme for parameter inference and model selection in dynamical systems

Abstract

Approximate Bayesian computation (ABC) methods can be used to evaluate posterior distributions without having to calculate likelihoods. In this paper, we discuss and apply an ABC method based on sequential Monte Carlo (SMC) to estimate parameters of dynamical models. We show that ABC SMC provides information about the inferability of parameters and model sensitivity to changes in parameters, and tends to perform better than other ABC approaches. The algorithm is applied to several well-known biological systems, for which parameters and their credible intervals are inferred. Moreover, we develop ABC SMC as a tool for model selection; given a range of different mathematical descriptions, ABC SMC is able to choose the best model using the standard Bayesian model selection apparatus.

Figure 1

(a) Trajectories of prey (solid curve) and predator (dashed curve) populations of the deterministic LV system and the data points (circles, prey data; triangles, predator data). (b) Parameters inferred by the ABC rejection sampler.

Figure 2

Histograms of populations (a) 0 (uniform prior distribution), (b) 1, (c) 2, (d) 3, (e) 4 and (f) 5 (approximation of posterior distribution) of parameters (i) a and (ii) b of the LV system.

Figure 3

Histograms of the approximated posterior distributions of parameters (a) c₁, (b) c₂ and (c) c₃ of the stochastic LV dynamics.

Figure 4

(a) Histograms of the approximate marginal posterior distributions of parameters α₀, n, β and α of the deterministic repressilator model. (b) The normalized 95% interquantile ranges (qr) of each population. The narrower the interval for a given tolerance ϵ_t, the more sensitive the model is to the corresponding parameter. The interquantile range reached in population 9 is determined by the added experimental noise. As ϵ₉ was chosen accordingly, one cannot proceed by lowering the tolerance further. The sharp change in the interquantile ranges, which occurs, for example, for parameter α₀ between populations 1 and 2, can be explained by the steep gradient of the likelihood surface along α₀. (c) The output (i.e. the accepted particles) of the ABC SMC algorithm as two-dimensional scatterplots. The particles from population 1 are in yellow, particles from population 4 in black, particles from population 7 in blue and those from the last population in red. Islands of particles are observed in population 4 and they can be explained by the multimodality of the fourth intermediate distribution. (d) PCA of the set of accepted particles (population 9). Owing to the dependence of the PCA on the scaling of original variables, the PCA was performed on the correlation matrix. The first PC explains 70.0% of the total variance, the second 24.6%, the third 5.3% and the fourth 0.1% of the variance. Pie charts show the fraction of the length of PCs explained by individual parameters.

Figure 5

(a) Histograms of the approximated posterior distributions of parameters α₀, n, β and α of the stochastic repressilator model. (b) The output of the ABC SMC algorithm as two-dimensional scatterplots. The particles from population 1 are in yellow, particles from population 2 in black, particles from population 3 in blue and those from the last population in red. We note that the projection to parameter n in population 2 is a sample from a multimodal intermediate distribution, while this distribution becomes unimodal from population 3 onwards.

Figure 6

Histograms show populations (a–i) 1–9 of the model parameter m. Population 9 represents the final posterior marginal estimate of m. Population 0 (discrete uniform prior distribution) is not shown.

Figure 7

Populations of the marginal posterior distribution of m. Models 1–4 correspond to equations (3.10)–(3.13), respectively. An interesting phenomenon is observed in populations 1–12, where model 2 has the highest probability, in contrast to model 3 having the highest inferred probability in the last population. The most probable explanation for this is that a local maximum favouring model 2 has been passed on route to a global maximum of the posterior probability favouring model 3. Populations (a–o) 1–15. Population 0 (discrete uniform prior distribution) is not shown.

Figure 8

Histograms of the approximated posterior distributions of parameters (a) γ, (b) v, (c) δ and (d) S(0) of the SIR model with a latent phase of infection (3.12).

Figure 9

(a) The probability density function of a mixture of two normal distributions, $(1/2)\varphi (0,1/100)+(1/2)\varphi (0,1)$, taken as a toy example used for a comparison of ABC SMC with ABC PRC. (b,c) Plots show how the variance of approximated intermediate distributions π_t changes with populations (t=1, …, 10 on the x-axis). The red curves plot the variance of the ABC SMC population and the blue curves the variance of the non-weighted ABC PRC populations. The perturbation kernel in both algorithms is uniform, $K_t=\sigma U(-1,1)$. In (b), σ=1.5 and this results in poor estimation of the posterior variance with the ABC PRC algorithm. In (c), σ is updated in each population so that it expands over the whole population range. Such σ is big enough for non-weighted ABC PRC to perform equally well as ABC SMC.