Inference for biochemical reaction systems based on a weak noise approximation
∂x . The marginal likelihood function �(x, t) is defined as the probability that a process with drift vector f(x) starting at state x and time t generated all observations yk with tk ≥ t. In theory, �(x, t) can be calculated exactly by solving the backward Fokker-Planck equation of the diffusion model with the end condition �(x, tn) = p(yn|x(tn)) and the jump conditions �(x, t − ) = �(x, t + ) p(yk|x(tk)). Here �(x, t − ) and �(x, t + ) denote the values of �(x, t) before and after the k-th observation. Instead we obtain an approximate solution in the weak noise limit by setting D → ǫD and assuming a scaling form �(x, t) ≈ z(t)exp[− 1 2ǫ (x −b(t)) ⊤ S −1 (t)(x −b(t))] for ǫ → 0. To leading order, we find ODEs for b(t), S(t), and z(t). By solving them backwards in time, we obtain approximations to the drift g(x, t) of the posterior process and the free energy F(θ) ≡ −ln p(y1, . . . , yn|θ), i.e. the negative log-likelihood of the data. The posterior statistics can now be computed either by applying the weak noise limit again (leading to an approximate Gaussian posterior process), or more simply by sampling. Additionally, type II maximum likelihood estimates of the parameters θ can be obtained by minimising F(θ) ≈ (2π) d/2 det(S(t1)) 1/2 z(t1). We have applied our method to the stochastic Lotka-Volterra model and an autoregulatory gene network [1]. In both cases the original process, which generated the observations, lies inside the inferred confidence interval most of the time. Thus our algorithm gives reliable information about the state of the reaction system, even if the internal noise is rather large. Parameter estimation also works well, as the inferred values are reasonable close to the real values of the reaction constants. We are currently working on full Bayesian inference by using a Hybrid Monte Carlo method, which samples from the posterior p(θ|y1, . . . , yn) ≈ 1 Z exp(−F(θ))p(θ).