M10: Chemical Fokker-Planck equation and multiscale modelling of (bio)chemical systems

Project completed September 30, 2012

Background

When modelling the biochemical networks that regulate the cell, the systems are traditionally simulated using deterministic (mean-field) models written in terms of ordinary differential equations describing the time evolution of chemical concentrations. Alternative methods include the well-known Stochastic Simulation Algorithm (SSA) attributed to Gillespie. This algorithm is able to take in to account what more traditional models cannot: the low abundancies of some of the key chemical species and the random or sporadic nature with which the reactions involving these species occur. However, this algorithm becomes computationally intractable because these types of systems are highly multiscale – that is, there are some reactions that occur multiple times on a timescale so short that other reactions in the system are unlikely to happen at all. Since the SSA simulates every single reaction that occurs in the system, this can make the analysis of the dynamics of the more slowly-varying quantities in the system computationally intractable.

Researchers at the Oxford Centre for Collaborative Applied Mathematics (OCCAM) have analysed these slowly changing quantities through the development of new multiscale methodologies and the solution of partial differential equations to approximate the invariant density of the system.

Techniques and Challenges

A range of multiscale methods already exist for the simulation of trajectories of the slow variables in a system, following some assumptions regarding the degree of time separation of the different variables. However in this project, the researchers approached this as a parameter estimation problem for an approximating stochastic differential equation. By estimating the effective drift and diffusion of the slow variables on a grid in the slow domain, many different properties of the system can be approximated, for example, the invariant density, average oscillation time, and switching time in multistable systems.

This project carefully considered the problem of approximating the solution of the steady-state chemical Fokker–Planck equation in more than two dimensions. Naive mesh generation with uniform refinement can lead to very inefficient meshes, which may cause either large errors or large computation times. Adaptive finite element methods can lead to more efficient meshes, which are more refined in regions that require it.

Once the invariant densities have been approximated, then post-processing tools need to be developed in order to be able to identify key features, such as multistability and oscillations. These analyses can in turn lead to detailed stochastic bifurcation diagrams for such systems.

Results

The researchers designed multiscale algorithms for finding coefficients of the Effective Fokker–Planck equation and demonstrated how these can be used to construct stochastic bifurcation diagrams for biochemical networks. Furthermore, they developed an adaptive finite element method in which the choice of mesh is assisted by stochastic trajectories of the Gillespie algorithm in order to solve efficiently and accurately Fokker–Planck equations in three dimensions.

A software package, the Stochastic Bifurcation Analysis Toolbox (SBAT), brings together several of the mathematical tools in this project. Currently, the SBAT includes implementations of the algorithms. The ultimate aim of the toolbox is automating a stochastic bifurcation analysis of well-mixed biochemical reaction networks.

The Future

This project represents just the start of the development of a suite of tools that will allow for the automated bifurcation analysis of stochastic descriptions of well-mixed biochemical reactions. SBAT 2.0 will be able to read in files describing generalised systems, and implement numerical multiscale reduction methods, adaptive finite element methods and post-processing tools in order to understand how biochemical networks’ behaviours can change as reaction parameters are altered.

References

[11/65] Cotter S.L., Vejchodský T., Erban R.: Adaptive finite element method assisted by stochastic simulation of chemical systems, submitted to SIAM Journal on Scientific Computing

[11/14] Cotter S.L., Zygalakis K.C., Kevrekidis I.G., Erban R.: A constrained approach to multiscale stochastic simulation of chemically reacting systems, Journal of Chemical Physics, Volume 135, 094102, 2011

[09/18] Erban R., Chapman S.J., Kevrekidis I.G., Vejchodsky T.: Analysis of a stochastic chemical system close to a SNIPER bifurcation of its mean-field model, SIAM Journal on Applied Mathematics, Volume 70, Number 3, pp. 984-1016, 2009

[1] Erban R., Chapman S.J., Maini P.K., A practical guide to stochastic simulations of reaction-diffusion processes, Lecture Notes, 2007

[2] Erban R., Kevrekidis I.G., Adalsteinsson D, Elston T.: Gene regulatory networks: a coarse-grained, equation-free approach to multiscale computation, Journal of Chemical Physics, Volume 124, Issue 8, 084106, 2006

[3] Singer A., Erban R., Kevrekidis I.G., Coifman R., Detecting the slow manifold by anisotropic diffusion maps, Proceedings of the National Academy of Sciences (PNAS), Volume 106, Number 38, pp. 16090-16095, 2009

[4] Erban R., Frewen T., Wang X., Elston T., Coifman R., Nadler B., Kevrekidis I.G., Variable-free exploration of stochastic models: a gene regulatory network example, Journal of Chemical Physics, Volume 126, Issue 15, 155103, 2007

This project is funded by the European Research Council Starting Independent Researcher Grant awarded to Dr Erban.