A new class of highly efficient exact stochastic simulation algorithms for chemical reaction networks

We introduce an alternative formulation of the exact stochastic simulation algorithm (SSA) for sampling trajectories of the chemical master equation for a well-stirred system of coupled chemical reactions. Our formulation is based on factored-out, partial reaction propensities. This novel exact SSA, called the partial propensity direct method (PDM), is highly efficient and has a computational cost that scales at most linearly with the number of chemical species, irrespective of the degree of coupling of the reaction network. In addition, we propose a sorting variant, SPDM, which is especially efficient for multiscale reaction networks.Comment: 23 pages, 3 figures, 4 tables; accepted by J. Chem. Phy

Switching and diffusion models for gene regulation networks

We analyze a hierarchy of three regimes for modeling gene regulation. The most complete model is a continuous time, discrete state space, Markov jump process. An intermediate 'switch plus diffusion' model takes the form of a stochastic differential equation driven by an independent continuous time Markov switch. In the third 'switch plus ODE' model the switch remains but the diffusion is removed. The latter two models allow for multi-scale simulation where, for the sake of computational efficiency, system components are treated differently according to their abundance. The 'switch plus ODE' regime was proposed by Paszek (Modeling stochasticity in gene regulation: characterization in the terms of the underlying distribution function, Bulletin of Mathematical Biology, 2007), who analyzed the steady state behavior, showing that the mean was preserved but the variance only approximated that of the full model. Here, we show that the tools of stochastic calculus can be used to analyze first and second moments for all time. A technical issue to be addressed is that the state space for the discrete-valued switch is infinite. We show that the new 'switch plus diffusion' regime preserves the biologically relevant measures of mean and variance, whereas the 'switch plus ODE' model uniformly underestimates the variance in the protein level. We also show that, for biologically relevant parameters, the transient behaviour can differ significantly from the steady state, justifying our time-dependent analysis. Extra computational results are also given for a protein dimerization model that is beyond the scope of the current analysis

Thermodynamic large fluctuations from uniformized dynamics

Large fluctuations have received considerable attention as they encode information on the fine-scale dynamics. Large deviation relations known as fluctuation theorems also capture crucial nonequilibrium thermodynamical properties. Here we report that, using the technique of uniformization, the thermodynamic large deviation functions of continuous-time Markov processes can be obtained from Markov chains evolving in discrete time. This formulation offers new theoretical and numerical approaches to explore large deviation properties. In particular, the time evolution of autonomous and non-autonomous processes can be expressed in terms of a single Poisson rate. In this way the uniformization procedure leads to a simple and efficient way to simulate stochastic trajectories that reproduce the exact fluxes statistics. We illustrate the formalism for the current fluctuations in a stochastic pump model

Fast stochastic simulation of biochemical reaction systems by\ud alternative formulations of the Chemical Langevin Equation

The Chemical Langevin Equation (CLE), which is a stochastic differential equation (SDE) driven by a multidimensional Wiener process, acts as a bridge between the discrete Stochastic Simulation Algorithm and the deterministic reaction rate equation when simulating (bio)chemical kinetics. The CLE model is valid in the regime where molecular populations are abundant enough to assume their concentrations change continuously, but stochastic fluctuations still play a major role. The contribution of this work is that we observe and explore that the CLE is not a single equation, but a parametric family of equations, all of which give the same finite-dimensional distribution of the variables. On the theoretical side, we prove that as many Wiener processes are sufficient to formulate the CLE as there are independent variables in the equation. On the practical side, we show that in the case where there are m1 pairs of reversible reactions and m2 irreversible reactions only m1+m2 Wiener processes are required in the formulation of the CLE, whereas the standard approach uses 2m1 + m2. We illustrate our findings by considering alternative formulations of the CLE for a\ud HERG ion channel model and the Goldbeter–Koshland switch. We show that there are considerable computational savings when using our insights

Chemical master versus chemical langevin for first-order reaction networks

Markov jump processes are widely used to model interacting species in circumstances where discreteness and stochasticity are relevant. Such models have been particularly successful in computational cell biology, and in this case, the interactions are typically rst-order. The Chemical Langevin Equation is a stochastic dierential equation that can be regarded as an approximation to the underlying jump process. In particular, the Chemical Langevin Equation allows simulations to be performed more eectively. In this work, we obtain expressions for the rst and second moments of the Chemical Langevin Equation for a generic rst-order reaction network. Moreover, we show that these moments exactly match those of the under-lying jump process. Hence, in terms of means, variances and correlations, the Chemical Langevin Equation is an excellent proxy for the Chemical Master Equation. Our work assumes that a unique solution exists for the Chemical Langevin Equation. We also show that the moment matching re- sult extends to the case where a gene regulation model of Raser and O'Shea (Science, 2004) is replaced by a hybrid model that mixes elements of the Master and Langevin equations. We nish with numerical experiments on a dimerization model that involves second order reactions, showing that the two regimes continue to give similar results