The chemical master equation is the well established description of stochastic chemical kinetics in well-mixed conditions. Exact solution of this equation is only known in a handful of cases, namely for systems composed of purely unimolecular reactions and for those obeying detailed balance. However biochemical systems are not in equilibrium (and hence rarely obey detailed balance) and bimolecular interactions abound (for example enzyme – substrate interactions). In the literature, these problems are typically circumvented by using the Gillespie algorithm, a Monte Carlo technique which exactly samples the trajectories of the chemical master equation. This method is feasible for small systems but becomes computationally expensive if one is interested in obtaining detailed statistical information, e.g. mean concentrations and coefficient of variation of the noise, over large regions of parameter space of reasonably sized biochemical networks.
We approach the problem by developing approximate analytical solutions of the chemical master equation. These methods lead to approximate algebraic expressions for the moments of the probability distribution solution of the chemical master equation as a function of the rate constants of the biochemical network; this opens the way to gaining deeper insight into the role played by stochasticity over the whole parameter space than possible via computationally expensive stochastic simulations. A considerable amount of effort in the group is devoted to developing novel approximation methods; two of these are Effective Mesoscopic Rate Equations (EMREs) and the slow-scale linear-noise approximation (ssLNA). These are briefly discussed below.
A. Effective Mesoscopic Rate Equations (EMREs)
EMREs are a new type of rate equations which are accurate over a wider range of molecule numbers than the conventional rate equations, i.e., they provide a more accurate estimate of the mean concentrations as predicted by the chemical master equation than the conventional rate equations. They are derived using the system size expansion of the master equation, in particular by considering terms one order higher than those leading to the classical linear noise approximation (PRL2009, JCP2010). EMREs have been applied to a large class of biochemical systems involving enzyme-substrate interactions (BMC2009, JCP2010), dimerization & trimerization of proteins (JCP2010, NatComm2012), and gene-protein interactions ( NatComm2012, PLoS2012, BIBM2012) and have been shown to lead to accurate results for systems characterized by as little as few tens of molecules. Various phenomena have been elucidated using EMREs including: amplification of substrate concentrations, concentration inversions and downstream blowup of conventional rate equation predictions.
Plots of the mean substrate concentrations versus time for sequential enzyme catalysis of two substrates S1 and S2. Note that the EMRE predictions agree well with the mean concentrations given by the stochastic simulation algorithm (SSA) but the rate equations (RE) predictions are way off. JCP 133, 035101 (2010).
B. Slow-Scale Linear-Noise Approximation (ssLNA)
It is common place to reduce the complexity of a biochemical system by grouping sets of elementary reactions together into an effective reaction and then writing an effective chemical master equation for the dynamics of these lumped reactions. However this method is flawed since the propensities of the lumped effective reactions are guessed based on the effective rates known from deterministic rate equations, e.g. a Michaelis-Menten rate in the rate equations leads to a Michaelis-Menten type propensity in the master equation. Starting from the chemical master equation of the full system of elementary reactions & assuming time scale separation conditions, we have rigorously derived in closed form a Langevin equation for the dynamics of the slowly varying species only (ssLNA) (PRE2012). The ssLNA provides a more accurate description of stochastic kinetics than the conventional effective master equation, in particular revealing how the latter misses out important phenomena such as noise-induced oscillations in genetic networks (BMC2012)
Sample paths of the stochastic simulation algorithm (SSA) for the full network of elementary reactions (gray), of the effective SSA for the reduced network of effective reactions (blue) and of the ssLNA Langevin method for the reduced network of effective reactions (red). Notice how the ssLNA sample paths show the same size of noise as those of the full network but the effective SSA does not. BMC 6:39 (2012)
We have an active interest in developing software, which puts the power of approximation methods at the fingertips of research groups involved in studying the role of stochasticity in biology. To this end we have recently developed iNA (PLoS2012, BIBM2012), an open source software with a user-friendly graphical user interface (GUI) which takes as input an SBML file describing a network (alternatively one can define the network from inside the software) and calculates the mean concentrations and variance / covariance of the noise in all chemical species using the linear-noise approximation (LNA) and higher orders of the system-size expansion (including EMREs). The system-size expansion is handled analytically using the powerful computer algebra Ginac and a just-in-time compiler, which compiles the mathematical model into platform specific machine code at runtime. This strategy yields an accurate fluctuation analysis, which typically outperforms that obtained by means of exact stochastic simulations. This software has been developed together with Hannes Matuschek (Interdisciplinary Center for Dynamics of Complex Systems, Univ. of Potsdam). iNA is free, runs on Mac, Linux and Windows, and can be downloaded from here
Screenshots of iNA's GUI. On the left we have the steady-state analysis panel where the rate equation (RE) and Effective Mesoscopic Rate Eq (EMRE) predictions are shown by blue and red bars respectively and the error bars show the variance calculated by the linear-noise approximation (LNA). On the right we have the time-course analysis where the means are calculated using REs and the standard deviation of intrinsic noise (shaded regions) are calculated using the LNA.
Concentration oscillations are a ubiquitous characteristic of intracellular dynamics. The deterministic oscillatory properties of various biochemical circuits have been extensively studied by means of the rate equation formalism. We are interested in understanding the oscillatory properties of stochastic biochemical systems, where the noise is of internal (molecular noise) or external origin (e.g. noise in the entraining stimulus). The linear-noise approximation, ssLNA (BMC2012) and EMREs formalisms (PLoS2012) provide a systematic means to perform the aforementioned analysis. Some of the predictions are being tested by means of experiments probing single cell circadian rhythms conducted by our collaborators in Edinburgh and abroad.
Stochastic phase space plot for the Brusselator generated using the linear-noise approximation. The Lambda's are non-dimensional parameters. NIO stands for noise-induced oscillations. The regions are as follows: 0 - deterministic limit cycle, 1 - stable node, NIO in species A, 2 - stable focus, NIO in species A, 4 - Stable Focus, NIO in species B, 6 - Stable Focus, NIO in both species A and B, 7 - Stable Node, No NIO, 8 - Stable Focus, no NIO.
The Gillespie algorithm and the equivalent chemical master equation are valid for well-stirred dilute environments. These conditions hold in certain regions of the intracellular environment but not in others. For the most part, the cell is a highly crowded (non-dilute) environment with an excluded fractional volume ranging from 0.1 to 0.4 depending on the location (BioPhysChem2006). We have a general interest in understanding how stochastic biochemical systems operate in crowded environments, in particular how crowding affects the mean concentrations, Fano factors and coefficients of variation of intrinsic noise. A central goal of our investigations is the derivation of effective chemical master equations, which capture the dynamics of Brownian dynamics simulations of stochastic reactions in crowded spaces. We have achieved this for simple systems (JCP2010) and work is in progress for extending these results to any reaction network of interest.
Plot of the normalized coefficient of variation in dimer fluctuations as a function of the occupied area fraction of a membrane. The solid points are obtained from Brownian dynamics of the dimerization reaction A+A ⇔ B whereas the solid line is the prediction of an effective chemical master equation taking into account crowding effects. JCP 132, 185102 (2010)
Due to the difficultly in exactly solving the chemical master equation, various approximation techniques have been developed over the years. A considerable number of these methods are based on ad-hoc assumptions and hence their regime of validity is frequently unclear. Examples of such methods include the non-linear Fokker Planck equation (FPE) which is obtained by a truncation of the Kramers-Moyal expansion to include at most second-order derivatives and moment closure approximations whereby one assumes some information about the underlying probability distribution to be solved for. We have developed a method based on the system-size expansion to estimate the error in the predictions of these approximation methods (applications to the non-linear FPE in JCP2011, and to normal moment-closure methodsin JCP2012). Applications to other approximation techniques and development of other methods of error-estimation are of current interest.
Dependence of the absolute relative error in the mean concentration and variance of intrinsic noise predictions using the nonlinear Fokker-Planck equation versus steady-state number of monomer molecules. This is for the open reaction system: null→ A, A+A→ null where A are the monomers. Red circles show the exact error computed using the exact solutions of the nonlinear FPE and of the chemical master equation while the solid blue line is estimated using a system-size expansion based method. JCP 135, 084103 (2011)