Stochastic modelling and simulation of biochemical reaction networks

The time-evolution of interacting types of molecules (e.g. genes, RNAs, proteins, metabolites) within a cell can be described by biochemical reaction networks. The population (no of molecules) of each variable (molecular specie) evolves in time by interacting with other (or same) variables-species of molecules through a set of coupled reactions.

​

Traditionally, the dynamics of these interactions are described by dynamical systems of ordinary differential equations which are deterministic (non-random) models. However, the recent biotechnological advances made clear that stochastic fluctuations in molecular populations are significant and can produce non-ignorable variability between genetically identical cells placed in seemingly identical experimental conditions. 

​

Stochastic processes are therefore more appropriate to describe the dynamics of these networks and continuous-time Markov-jump processes, governed by the Kolmogorov (Chemical Master) equation are typically considered a fine description of the underlying process. However, their Kolmogorov equations have analytic solutions in the simplest cases, and while it is fairly easy to construct algorithms to simulate them exactly (see Gillespie's SSA), these simulations are often very slow for any bio-molecular network of practical interest.  

​

Therefore, in order to predict, analyse, design experiments, and infer the underlying process one needs to rely on models that approximate these Markov Jump processes. The target is to develop fast, and accurate approximations. Some methods focus on developing simulation algorithms (see tau-leaping, Chemical Langevin (diffusion) equation) that speed up simulation, while another method, namely the linear noise approximation (LNA), also provide analytical expressions for the probability distributions of the state of the process and therefore not only speeds up simulation, but also provide means of formal analysis of the system. 

 

While the LNA is accurate in some fairly simple situations, we have observed that it quickly becomes inaccurate in oscillatory-periodic systems. It is also inappropriate for bi-stable and other types of systems. 

​

We recently developed a method that applies a phase correction to the LNA (pcLNA) to provide very fast simulations (100s to 1000s of times faster than SSA), and analytical expressions for probability distributions of the state of the system at given phases, which are both shown to be indistinguishable from SSA empirical distributions.  

​

We currently work on developing models and simulation algorithms for other types of dynamics that are commonly observed in cellular biology and elsewhere.