Stochastic simulation

Stochastic simulation algorithms and methods were initially developed to analyse chemical reactions involving large numbers of species with complex reaction kinetics. The first algorithm, the Gillespie algorithm was proposed by Dan Gillespie in 1977. It is an exact procedure for numerically simulating the time evolution of a well-stirred chemically reacting system.

The algorithm is a Monte Carlo type method.

Discrete, exact variants

In order to determine the next event in a stochastic simulation, the rates of all possible changes to the state of the model are computed, and then ordered in an array. Next, the cumulative sum of the array is taken, and the final cell contains the number R, where R is the total event rate. This cumulative array is now a discrete cumulative distribution, and can be used to choose the next event by picking a random number z~U(0,R) and choosing the first event, such that z is less than the rate associated with that event.

In historical order

Direct and first reaction methods

Published by Dan Gillespie in 1977, and is a linear search on the cumulative array. See Gillespie algorithm.

Next reaction method

Published 2000. This is an improvement over the first reaction method where the unused reaction times are reused. To make the sampling of reactions more efficient, an indexed priority queue is used to store the reaction times. On the other hand, to make the recomputation of propensities more efficient, a dependency graph is used. This dependency graph tells which reaction propensities to update after a particular reaction has fired.

Optimised and sorting direct methods

Published 2004 and 2005. These methods sort the cumulative array to reduce the average search depth of the algorithm. The former runs a presimulation to estimate the firing frequency of reactions, whereas the latter sorts the cumulative array on-the-fly.

Logarithmic direct method

Published in 2006. This is a binary search on the cumulative array, thus reducing the worst-case time complexity of reaction sampling to O(log M).

Partial-propensity methods

Published in 2009, 2010, and 2011 (Ramaswamy 2009, 2010, 2011). Use factored-out, partial reaction propensities to reduce the computational cost to scale with the number of species in the network, rather than the (larger) number of reactions. Four variants exist:

• PDM, the partial-propensity direct method. Has a computational cost that scales linearly with the number of different species in the reaction network, independent of the coupling class of the network (Ramaswamy 2009).
• SPDM, the sorting partial-propensity direct method. Uses dynamic bubble sort to reduce the pre-factor of the computational cost in multi-scale reaction networks where the reaction rates span several orders of magnitude (Ramaswamy 2009).
• PSSA-CR, the partial-propensity SSA with composition-rejection sampling. Reduces the computational cost to constant time (i.e., independent of network size) for weakly coupled networks (Ramaswamy 2010) using composition-rejection sampling (Slepoy 2008).
• dPDM, the delay partial-propensity direct method. Extends PDM to reaction networks that incur time delays (Ramaswamy 2011) by providing a partial-propensity variant of the delay-SSA method (Bratsun 2005, Cai 2007).

The use of partial-propensity methods is limited to elementary chemical reactions, i.e., reactions with at most two different reactants. Every non-elementary chemical reaction can be equivalently decomposed into a set of elementary ones, at the expense of a linear (in the order of the reaction) increase in network size.

Continuous, approximate variants

τ leap and modified Poisson τ leap methods

First published in 2001; modified in 2005.

See also

• Gillespie algorithm
• Network simulation
• Network traffic simulation
• Simulation language
• Queueing theory

References

• (Slepoy 2008): A. Slepoy, A. P. Thompson, and S. J. Plimpton,, A; Thompson, AP; Plimpton, SJ (2008). "A constant-time kinetic Monte Carlo algorithm for simulation of large biochemical reaction networks". Journal of Chemical Physics 128 (20): 205101. doi:10.1063/1.2919546. PMID 18513044. 
• (Bratsun 2005): D. Bratsun, D. Volfson, J. Hasty, L. Tsimring, (2005). "Delay-induced stochastic oscillations in gene regulation". PNAS 102 (41): 14593–8. doi:10.1073/pnas.0503858102. PMC 1253555. PMID 16199522. 
• (Cai 2007): X. Cai, (2007). "Exact stochastic simulation of coupled chemical reactions with delays". J. Chem. Phys. 126 (12): 124108. doi:10.1063/1.2710253. PMID 17411109. 
• Hartmann, A.K. (2009). Practical Guide to Computer Simulations. World Scientific. ISBN 978-981-283-415-7. 
• Press, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007). "Section 17.7. Stochastic Simulation of Chemical Reaction Networks". Numerical Recipes: The Art of Scientific Computing (3rd ed.). New York: Cambridge University Press. ISBN 978-0-521-88068-8. 
• (Ramaswamy 2009): R. Ramaswamy, N. Gonzalez-Segredo, I. F. Sbalzarini, (2009). "A new class of highly efficient exact stochastic simulation algorithms for chemical reaction networks". J. Chem. Phys. 130 (24): 244104. doi:10.1063/1.3154624. PMID 19566139. 
• (Ramaswamy 2010): R. Ramaswamy, I. F. Sbalzarini, (2010). "A partial-propensity variant of the composition-rejection stochastic simulation algorithm for chemical reaction networks". J. Chem. Phys. 132 (4): 044102. doi:10.1063/1.3297948. PMID 20113014. 
• (Ramaswamy 2011): R. Ramaswamy, I. F. Sbalzarini, (2011). "A partial-propensity formulation of the stochastic simulation algorithm for chemical reaction networks with delays". J. Chem. Phys. 134 (1): 014106. doi:10.1063/1.3521496. PMID 21218996. 

External links

Software

• Cain - Stochastic simulation of chemical kinetics. Direct, next reaction, tau-leaping, hybrid, etc.
• PDM - C++ implementations of all partial-propensity methods.
• StochPy - Stochastic modelling in Python
• STEPS - STochastic Engine for Pathway Simulation using swig to create Python interface to C/C++ code