Gillespie algorithm
It was created by Joseph L. Doob and others (circa 1945), presented by Dan Gillespie in 1976, and popularized in 1977 in a paper where he uses it to simulate chemical or biochemical systems of reactions efficiently and accurately using limited computational power (see stochastic simulation).
[1] As computers have become faster, the algorithm has been used to simulate increasingly complex systems.
[citation needed] The process that led to the algorithm recognizes several important steps.
It was William Feller, in 1940, who found the conditions under which the Kolmogorov equations admitted (proper) probabilities as solutions.
In his Theorem I (1940 work) he establishes that the time-to-the-next-jump was exponentially distributed and the probability of the next event is proportional to the rate.
As such, he established the relation of Kolmogorov's equations with stochastic processes.
Later, Doob (1942, 1945) extended Feller's solutions beyond the case of pure-jump processes.
Gillespie (1977) obtains the algorithm in a different manner by making use of a physical argument.
The rate is determined by the number of molecules in each chemical species.
The key assumptions are that Given the two assumptions, the random waiting time for some reaction is exponentially distributed, with exponential rate being the sum of the individual reaction's rates.
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.
A trajectory corresponding to a single Gillespie simulation represents an exact sample from the probability mass function that is the solution of the master equation.
The physical basis of the algorithm is the collision of molecules within a reaction vessel.
A review (Gillespie, 2007) outlines three different, but equivalent formulations; the direct, first-reaction, and first-family methods, whereby the former two are special cases of the latter.
The formulation of the direct and first-reaction methods is centered on performing the usual Monte Carlo inversion steps on the so-called "fundamental premise of stochastic chemical kinetics", which mathematically is the function where each of the
is "a statistically independent integer random variable with point probabilities
Thus, the Monte Carlo generating method is simply to draw two pseudorandom numbers,
This family of algorithms 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.
The computational cost of exact versions of the algorithm is determined by the coupling class of the reaction network.
An exact version of the algorithm with constant-time scaling for weakly coupled networks has been developed, enabling efficient simulation of systems with very large numbers of reaction channels (Slepoy Thompson Plimpton 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 and independently Barrio et al. 2006, as well as (Cai 2007).
Partial-propensity formulations, as developed independently by both Ramaswamy et al. (2009, 2010) and Indurkhya and Beal (2010), are available to construct a family of exact versions of the algorithm whose computational cost is proportional to the number of chemical species in the network, rather than the (larger) number of reactions.
These formulations can reduce the computational cost to constant-time scaling for weakly coupled networks and to scale at most linearly with the number of species for strongly coupled networks.
A partial-propensity variant of the generalized Gillespie algorithm for reactions with delays has also been proposed (Ramaswamy Sbalzarini 2011).
, then the time, δt, until the next reaction occurs is a random number drawn from exponential distribution function with mean
The Gillespie algorithm just repeats these two steps as many times as needed to simulate the system for however long we want (i.e., for as many reactions).
dimers and 2 of A and B but due to the small numbers of molecules fluctuations around these values are large.
The Gillespie algorithm is often used to study systems where these fluctuations are important.
