Gillespie algorithm

The Gillespie algorithm generates a statistically correct trajectory (possible solution) of a stochastic equation. It was developed and published by Dan Gillespie in 1977 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 for now):
# Initialization: Initialize the number of molecules in the system, reactions constants, and random number generators.
# Monte Carlo Step: Generate random numbers to determine the next reaction to occur as well as the time interval.
# Update: Increase the time step by the randomly generated time in Step 1. Update the molecule count based on the reaction that occurred.
# Iterate: Go back to Step 1 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). See the articles cited below for details.

Further reading

*cite journal| author=Daniel T. Gillespie | title=Exact Stochastic Simulation of Coupled Chemical Reactions | journal=The Journal of Physical Chemistry| volume=81 |issue=25| pages= 2340-2361 | year=1977 | doi = 10.1021/j100540a008
*cite journal| author=Daniel T. Gillespie | title=A General Method for Numerically Simulating the Stochastic Time Evolution of Coupled Chemical Reactions| journal= Journal of Computational Physics | volume=22| issue=4 | pages=403-434| year=1976 | doi = 10.1016/0021-9991(76)90041-3
*cite journal| author=M. Rathinam, L. R. Petzold, Y. Cao, and Daniel T. Gillespie, | title=Stiffness in stochastic chemically reacting systems: The implicit tau-leaping method | journal= Journal of Chemical Physics| volume= 119| issue=24 | pages= 12784-12794 | year=2003| doi = 10.1063/1.1627296
*cite journal| author=M. A. Gibson and J. Bruck, | url=http://www.cs.caltech.edu/courses/cs191/paperscs191/JPhysChemA(2000-104)1876.pdf | title=Efficient Exact Stochastic Simulation of Chemical Systems with Many Species and Many Channels | journal= J. Phys. Chem. A | volume= 104 | pages= 1876-1889 | year=2000 | doi = 10.1021/jp993732q
*cite journal| author=H. Salis, and Y. N. Kaznessis, | title=Accurate hybrid stochastic simulation of a system of coupled chemical or biochemical reactions | journal= Journal of Chemical Physics| volume= 122 | pages= 054103 | year=2005| doi = 10.1063/1.1835951

* (Slepoy 2008): cite journal| author=A. Slepoy, A. P. Thompson, and S. J. Plimpton, | title=A constant-time kinetic Monte Carlo algorithm for simulation of large biochemical reaction networks | journal= Journal of Chemical Physics| volume= 128| issue=20 | pages= 205101 | year=2008| doi = 10.1063/1.2919546