Gillespie algorithm

In probability theory, the Gillespie algorithm (or occasionally the Doob-Gillespie algorithm) generates a statistically correct trajectory (possible solution) of a stochastic equation. It was created by Joseph L. Doob in 1945 and developed and popularized by Dan Gillespie in 1977 in a paper where he uses it to simulate chemical or biochemical systems of reactions efficiently and accurately using limited computational power. As computers have become faster, the algorithm has been used to simulate increasingly complex systems. The algorithm is particularly useful for simulating reactions within cells where the number of reagents typically number in the tens of molecules (or less). Mathematically, it is a variety of a dynamic Monte Carlo method and similar to the kinetic Monte Carlo methods. It is used heavily in computational systems biology.

Idea behind the algorithm

Traditional continuous and deterministic biochemical rate equations do not accurately predict cellular reactions since they rely on bulk reactions that require the interactions of millions of molecules. They are typically modeled as a set of coupled ordinary differential equations. In contrast, the Gillespie algorithm allows a discrete and stochastic simulation of a system with few reactants because every reaction is explicitly simulated. When simulated, a Gillespie realization represents a random walk that exactly represents the distribution of the master equation.

The physical basis of the algorithm is the collision of molecules within a reaction vessel. It is assumed that collisions are frequent, but collisions with the proper orientation and energy are infrequent. Therefore, all reactions within the Gillespie framework must involve at most two molecules. Reactions involving three molecules are assumed to be extremely rare and are modeled as a sequence of binary reactions. It is also assumed that the reaction environment is well mixed.

Algorithm

Gillespie developed two different, but equivalent formulations; the direct method and the first reaction method. Below is a summary of the steps to run the algorithm (math omitted):

1. Initialization: Initialize the number of molecules in the system, reactions constants, and random number generators.
2. Monte Carlo step: Generate random numbers to determine the next reaction to occur as well as the time interval. The probability of a given reaction to be chosen is proportional to the number of substrate molecules.
3. Update: Increase the time step by the randomly generated time in Step 2. Update the molecule count based on the reaction that occurred.
4. Iterate: Go back to Step 2 unless the number of reactants is zero or the simulation time has been exceeded.

The algorithm is computationally expensive and thus many modifications and adaptations exist, including the next reaction method (Gibson & Bruck), tau-leaping, as well as hybrid techniques where abundant reactants are modeled with deterministic behavior. Adapted techniques generally compromise the exactitude of the theory behind the algorithm as it connects to the Master equation, but offer reasonable realizations for greatly improved timescales. Recently, an exact version of the algorithm with constant-time scaling has been developed enabling efficient simulation of systems with very large numbers of reaction channels (Slepoy 2008). The generalized Gillespie algorithm that accounts for the non-Markovian properties of random biochemical events with delay has been developed by Bratsun et al. 2005. See the articles cited below for details.

Further reading

• A Kolmogorov (1931). "Über die analytisehen Methoden in der Wahrseheinliehkeitsreehnung" (PDF). Mathematische Annalen. 104: 415.
• W Feller (1940). "On the Integro-Differential Equations of Purely Discontinous Markoff Processes". Transactions of the American Mathematical Society. 48: 4885–15.
• Jacob L. Doob (1942). "Topics in the Theory of Markoff Chains". Transactions of the American Mathematical Society. 52 (1): 37–64.
• Jacob L. Doob (1945). "Markoff chains – Denumerable case". Transactions of the American Mathematical Society. 58 (3): 455–473. doi:10.2307/1990339.
• David G. Kendall (1950). "An Artificial Realization of a Simple "Birth-and-Death" Process". Journal of the Royal Statistical Society. Series B (Methodological),. 12 (1): 116–119.{{cite journal}}: CS1 maint: extra punctuation (link)
• Maurice S. Bartlett (1953). "Stochastic Processes or the Statistics of Change". : Journal of the Royal Statistical Society. Series C (Applied Statistics),. 2 (1): 44–64.{{cite journal}}: CS1 maint: extra punctuation (link)
• Daniel T. Gillespie (1977). "Exact Stochastic Simulation of Coupled Chemical Reactions". The Journal of Physical Chemistry. 81 (25): 2340–2361. doi:10.1021/j100540a008.
• Daniel T. Gillespie (1976). "A General Method for Numerically Simulating the Stochastic Time Evolution of Coupled Chemical Reactions". Journal of Computational Physics. 22 (4): 403–434. doi:10.1016/0021-9991(76)90041-3.
• M. Rathinam, L. R. Petzold, Y. Cao, and Daniel T. Gillespie, (2003). "Stiffness in stochastic chemically reacting systems: The implicit tau-leaping method". Journal of Chemical Physics. 119 (24): 12784–12794. doi:10.1063/1.1627296.{{cite journal}}: CS1 maint: extra punctuation (link) CS1 maint: multiple names: authors list (link)
• Sinitsyn, NA; Hengartner, N; Nemenman, I (2009). "Adiabatic coarse-graining and simulations of stochastic biochemical networks" (PDF). Proceed. Nat. Acad. Science U. S. A. 106 (20): 10546–10551. doi:10.1073/pnas.0809340106. PMC 2705573. PMID 19525397. {{cite journal}}: More than one of |author= and |last1= specified (help)CS1 maint: extra punctuation (link)
• M. A. Gibson and J. Bruck, (2000). "Efficient Exact Stochastic Simulation of Chemical Systems with Many Species and Many Channels" (PDF). J. Phys. Chem. A. 104: 1876–1889. doi:10.1021/jp993732q.{{cite journal}}: CS1 maint: extra punctuation (link)
• Salis, H; Kaznessis, Y (2005). "Accurate hybrid stochastic simulation of a system of coupled chemical or biochemical reactions". Journal of Chemical Physics. 122 (5): 054103. doi:10.1063/1.1835951. PMID 15740306. {{cite journal}}: More than one of |author= and |last1= specified (help)CS1 maint: extra punctuation (link)

• (Slepoy 2008): Slepoy, 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. {{cite journal}}: More than one of |author= and |last1= specified (help)CS1 maint: extra punctuation (link)

• (Bratsun 2005): D. Bratsun, D. Volfson, J. Hasty, L. Tsimring, (2005). "Delay-induced stochastic oscillations in gene regulation". PNAS. 102 (41): 14593. doi:10.1073/pnas.0503858102. PMC 1253555. PMID 16199522.{{cite journal}}: CS1 maint: extra punctuation (link) CS1 maint: multiple names: authors list (link)

External links

Software

• Cain - Stochastic simulation of chemical kinetics. Direct, next reaction, tau-leaping, hybrid, etc.
• SynBioSS - Stochastic simulation of chemical kinetics using the exact SSA as well as an SSA/Langevin hybrid. Both MPI-parallel (supercomputer) and GUI (desktop) versions are provided.
• GillespieSSA - R package for Gillespie algorithm
• [1] - Mathematica code and applet for stochastic simulation of chemical kinetics.