Gillespie algorithm

In probability theory and computational science, the Gillespie algorithm, also known as the stochastic simulation algorithm (SSA), is a method for generating statistically exact sample trajectories of certain continuous-time Markov jump processes. It is especially associated with the simulation of coupled chemical reaction systems. Daniel Gillespie presented the method in 1976 and developed it further in 1977.^[1]^[2]

The method is widely used in computational systems biology, particularly when the numbers of reacting molecules are small enough that stochastic fluctuations are important.^[3] Gillespie-type algorithms are also used in network epidemiology to simulate stochastic epidemic spreading on complex networks, including continuous-time SI, SIS, and SIR processes.^[4]

Mathematically, the algorithm is a form of dynamic Monte Carlo method and is closely related to kinetic Monte Carlo methods.

History

The mathematical foundations of the Gillespie algorithm lie in the theory of stochastic jump processes. In 1931, Andrey Kolmogorov developed differential equations describing the time evolution of Markov processes, now known as the Kolmogorov equations.^[5]

William Feller subsequently studied purely discontinuous Markov processes and established conditions under which the Kolmogorov equations yield proper transition probabilities.^[6] Joseph L. Doob extended the mathematical theory of Markov processes during the 1940s.^[7]^[8]

Early computer simulations of related stochastic processes predated Gillespie's work. David George Kendall simulated a simple birth-and-death process using the Manchester Mark 1 computer in 1950,^[9] and Maurice S. Bartlett applied stochastic-process methods to epidemic modelling.^[10]

Gillespie's 1976 and 1977 papers derived an exact simulation procedure for systems of coupled chemical reactions using probabilistic arguments based on the physical interpretation of molecular reaction events.^[1]^[2]

Method

Mathematical formulation

Consider a system whose state at time $t$ is described by a vector ${\boldsymbol {x}}$ . Each possible event or reaction $j$ has a propensity function $a_{j}({\boldsymbol {x}})$ . For chemical systems, the propensity represents the instantaneous rate at which a particular reaction can occur when the system is in state ${\boldsymbol {x}}$ .^[3]

Define the total propensity

a_{0}({\boldsymbol {x}})=\sum _{j}a_{j}({\boldsymbol {x}}).

Provided that the process is Markovian, the waiting time $\tau$ until the next event is exponentially distributed with rate $a_{0}({\boldsymbol {x}})$ . The joint probability density that the next event occurs after a waiting time $\tau$ and is reaction $j$ is

p(\tau ,j\mid {\boldsymbol {x}},t)=a_{j}({\boldsymbol {x}})\exp[-a_{0}({\boldsymbol {x}})\tau ].

Thus, the simulation requires two random choices: the time of the next event and which event occurs.

Using two independent pseudorandom numbers $r_{1}$ and $r_{2}$ , uniformly distributed on $(0,1)$ , the waiting time can be generated as

\tau ={\frac {1}{a_{0}({\boldsymbol {x}})}}\ln \left({\frac {1}{r_{1}}}\right).

The next reaction $j$ is the smallest integer satisfying

\sum _{j'=1}^{j}a_{j'}({\boldsymbol {x}})>r_{2}a_{0}({\boldsymbol {x}}).

Direct method

The direct form of the Gillespie algorithm can be written as follows:^[2]^[3]

Initialize the time $t=t_{0}$ and state ${\boldsymbol {x}}={\boldsymbol {x}}_{0}$ .
Evaluate each propensity $a_{j}({\boldsymbol {x}})$ and their total $a_{0}({\boldsymbol {x}})$ .
Generate the waiting time $\tau$ and select the next reaction $j$ .
Advance the time according to $t\leftarrow t+\tau$ .
Update the state according to ${\boldsymbol {x}}\leftarrow {\boldsymbol {x}}+{\boldsymbol {\nu }}_{j}$ , where ${\boldsymbol {\nu }}_{j}$ is the state-change vector associated with reaction $j$ .
Record the state if required and repeat until the chosen stopping condition is reached.

Each realization produced in this way is a sample trajectory of the corresponding stochastic process.

Comparison with deterministic models

Traditional deterministic descriptions of chemical kinetics generally use systems of coupled ordinary differential equations for continuous concentrations. Such models are often effective when populations are large and stochastic fluctuations are relatively small.

When some molecular species occur in low copy numbers, however, fluctuations associated with individual reaction events may become significant. The Gillespie algorithm represents molecules and reaction events discretely and therefore provides a way of simulating these intrinsic stochastic fluctuations.^[3]

For the standard chemical interpretation of the algorithm, the reaction system is normally assumed to be sufficiently well mixed that spatial correlations can be neglected and that the propensity of each reaction depends only on the current state.

Computational variants

The original direct method may become computationally expensive when a model contains a large number of reaction channels. Consequently, numerous exact and approximate variants have been developed.

The next reaction method of Gibson and Bruck uses data structures such as dependency graphs and priority queues to reduce the computational effort needed to determine subsequent reaction events.^[11]

Approximate methods such as tau-leaping advance the system through multiple reaction events over a finite time interval rather than simulating every individual event, potentially providing substantial computational savings.^[12] Hybrid methods can likewise combine stochastic descriptions of low-copy-number species with deterministic approximations for abundant species.

For weakly coupled reaction networks, exact algorithms with computational cost independent of the total number of reaction channels under appropriate conditions have been developed.^[13]

Extensions have also been proposed for systems in which events are delayed or otherwise depart from the assumptions of a simple Markov reaction process.^[14]^[15]^[16]

Partial-propensity methods

Partial-propensity formulations reduce the computational work needed to evaluate and select reaction channels by exploiting the structure of their propensity functions.^[17]^[18]

Related approaches based on reaction factoring and dependency structures have also been proposed for large biochemical systems.^[19]

Partial-propensity techniques have additionally been extended to reaction networks containing delays.^[20]

Applications

Chemical and biochemical reaction networks

The Gillespie algorithm was originally developed for coupled chemical reaction systems and remains widely used to investigate stochastic chemical and biochemical dynamics.^[2]^[3] Applications include systems in which fluctuations caused by small numbers of molecules may influence the system's behaviour, such as gene-expression and regulatory networks.

Epidemic processes on networks

Epidemic models can also be represented as continuous-time stochastic event systems. For a network-based epidemic model, events may include transmission across an edge between susceptible and infected nodes and recovery of an infected node.

Gillespie-based simulation can therefore be used to generate continuous-time realizations of epidemic models on complex contact networks. Kuryliak, Emmerich, and Dosyn applied efficient Gillespie implementations to SI, SIS, and SIR dynamics on heterogeneous networks and studied how network structure and interventions influence epidemic peak size and timing.^[4]

Example

Reversible binding of A and B to form AB dimers

A simple chemical system illustrates the direct method. Consider molecules of two types, $A$ and $B$ , which reversibly bind to form $AB$ dimers:^[3]

{\ce {A + B <=> AB}}

There are two possible reactions. An $A$ molecule and a $B$ molecule can combine to produce an $AB$ dimer, or an existing dimer can dissociate.

Let the rate constant for dimer formation be $k_{\mathrm {D} }$ and the rate constant for dissociation be $k_{\mathrm {B} }$ . If the system contains $n_{\mathrm {A} }$ molecules of $A$ , $n_{\mathrm {B} }$ molecules of $B$ , and $n_{\mathrm {AB} }$ dimers, the two propensities are

a_{1}=k_{\mathrm {D} }n_{\mathrm {A} }n_{\mathrm {B} }

and

a_{2}=k_{\mathrm {B} }n_{\mathrm {AB} }.

The total reaction propensity is therefore

a_{0}=k_{\mathrm {D} }n_{\mathrm {A} }n_{\mathrm {B} }+k_{\mathrm {B} }n_{\mathrm {AB} }.

The waiting time until the next reaction is sampled from an exponential distribution with mean $1/a_{0}$ . Once the waiting time has been generated, the probability that the next event is dimer formation is

P(\mathrm {A+B\rightarrow AB} )={\frac {k_{\mathrm {D} }n_{\mathrm {A} }n_{\mathrm {B} }}{a_{0}}},

while the probability of dissociation is

P(\mathrm {AB\rightarrow A+B} )={\frac {k_{\mathrm {B} }n_{\mathrm {AB} }}{a_{0}}}.

If dimer formation occurs, $n_{\mathrm {A} }$ and $n_{\mathrm {B} }$ each decrease by one and $n_{\mathrm {AB} }$ increases by one. If dissociation occurs, the opposite update is made. Repeatedly sampling the waiting time and reaction type generates one stochastic trajectory of the system.

References

1 2 Gillespie, Daniel T. (1976). "A General Method for Numerically Simulating the Stochastic Time Evolution of Coupled Chemical Reactions". Journal of Computational Physics. 22 (4): 403–434. Bibcode:1976JCoPh..22..403G. doi:10.1016/0021-9991(76)90041-3.
1 2 3 4 Gillespie, Daniel T. (1977). "Exact Stochastic Simulation of Coupled Chemical Reactions". The Journal of Physical Chemistry. 81 (25): 2340–2361. doi:10.1021/j100540a008. S2CID 2606191.
1 2 3 4 5 6 Gillespie, Daniel T. (2007-05-01). "Stochastic Simulation of Chemical Kinetics". Annual Review of Physical Chemistry. 58 (1): 35–55. Bibcode:2007ARPC...58...35G. doi:10.1146/annurev.physchem.58.032806.104637. PMID 17037977.
1 2 Kuryliak, Yulian; Emmerich, Michael T. M.; Dosyn, Dmytro (2025). "Simulating epidemic peak dynamics on complex networks using efficient Gillespie algorithms". Infection, Genetics and Evolution. 132 105768. doi:10.1016/j.meegid.2025.105768. PMID 40441473.
↑ Kolmogorov, Andrey N. (1931). "Über die analytischen Methoden in der Wahrscheinlichkeitsrechnung" [On Analytical Methods in the Theory of Probability]. Mathematische Annalen. 104: 415–458. doi:10.1007/BF01457949. S2CID 119439925.
↑ Feller, William (1940). "On the Integro-Differential Equations of Purely Discontinuous Markoff Processes". Transactions of the American Mathematical Society. 48 (3): 488–515. doi:10.2307/1990095. JSTOR 1990095.
↑ Doob, Joseph L. (1942). "Topics in the Theory of Markoff Chains". Transactions of the American Mathematical Society. 52 (1): 37–64. doi:10.1090/S0002-9947-1942-0006633-7. JSTOR 1990152.
↑ Doob, Joseph L. (1945). "Markoff Chains—Denumerable Case". Transactions of the American Mathematical Society. 58 (3): 455–473. doi:10.2307/1990339. JSTOR 1990339.
↑ Kendall, David G. (1950). "An Artificial Realization of a Simple "Birth-and-Death" Process". Journal of the Royal Statistical Society, Series B. 12 (1): 116–119. doi:10.1111/j.2517-6161.1950.tb00048.x. JSTOR 2983837.
↑ Bartlett, Maurice S. (1953). "Stochastic Processes or the Statistics of Change". Journal of the Royal Statistical Society, Series C. 2 (1): 44–64. doi:10.2307/2985327. JSTOR 2985327.
↑ Gibson, Michael A.; Bruck, Jehoshua (2000). "Efficient Exact Stochastic Simulation of Chemical Systems with Many Species and Many Channels". Journal of Physical Chemistry A. 104 (9): 1876–1889. Bibcode:2000JPCA..104.1876G. doi:10.1021/jp993732q.
↑ Rathinam, Muruhan; Petzold, Linda R.; Cao, Yang; Gillespie, Daniel T. (2003). "Stiffness in stochastic chemically reacting systems: The implicit tau-leaping method". Journal of Chemical Physics. 119 (24): 12784–12794. Bibcode:2003JChPh.11912784R. doi:10.1063/1.1627296.
↑ Slepoy, Alexander; Thompson, Aidan P.; Plimpton, Steven J. (2008). "A constant-time kinetic Monte Carlo algorithm for simulation of large biochemical reaction networks". Journal of Chemical Physics. 128 (20): 205101. Bibcode:2008JChPh.128t5101S. doi:10.1063/1.2919546. PMID 18513044.
↑ Bratsun, Dmitri; Volfson, Dmitri; Hasty, Jeff; Tsimring, Lev S. (2005). "Delay-induced stochastic oscillations in gene regulation". Proceedings of the National Academy of Sciences. 102 (41): 14593–14598. doi:10.1073/pnas.0503858102. PMC 1253555. PMID 16199522.
↑ Barrio, Manuel; Burrage, Kevin; Leier, André; Tian, Tianhai (2006). "Oscillatory Regulation of hes1: Discrete Stochastic Delay Modelling and Simulation". PLOS Computational Biology. 2 (9) e117. doi:10.1371/journal.pcbi.0020117. PMC 1560403. PMID 16965175.
↑ Cai, Xiaodong (2007). "Exact stochastic simulation of coupled chemical reactions with delays". Journal of Chemical Physics. 126 (12): 124108. Bibcode:2007JChPh.126l4108C. doi:10.1063/1.2710253. PMID 17411109.
↑ Ramaswamy, Rajesh; González-Segredo, Nélido; Sbalzarini, Ivo F. (2009). "A new class of highly efficient exact stochastic simulation algorithms for chemical reaction networks". Journal of Chemical Physics. 130 (24): 244104. arXiv:0906.1992. doi:10.1063/1.3154624. PMID 19566139.
↑ Ramaswamy, Rajesh; Sbalzarini, Ivo F. (2010). "A partial-propensity variant of the composition-rejection stochastic simulation algorithm for chemical reaction networks". Journal of Chemical Physics. 132 (4): 044102. doi:10.1063/1.3297948. PMID 20113014.
↑ Indurkhya, Sagar; Beal, Jacob S. (2010). "Reaction Factoring and Bipartite Update Graphs Accelerate the Gillespie Algorithm for Large-Scale Biochemical Systems". PLOS ONE. 5 (1) e8125. doi:10.1371/journal.pone.0008125. PMC 2798956. PMID 20066048.
↑ Ramaswamy, Rajesh; Sbalzarini, Ivo F. (2011). "A partial-propensity formulation of the stochastic simulation algorithm for chemical reaction networks with delays". Journal of Chemical Physics. 134 (1): 014106. doi:10.1063/1.3521496. PMID 21218996.

A page that knows where it met you.