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.

The method is widely used in computational systems biology, particularly when the numbers of reacting molecules are small enough that stochastic fluctuations are important. 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.

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.

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

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, and Maurice S. Bartlett applied stochastic-process methods to epidemic modelling.

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.

Method

Mathematical formulation

Consider a system whose state at time $t$ is described by a vector $𝒙$. Each possible event or reaction $j$ has a propensity function $a_j(𝒙)$. For chemical systems, the propensity represents the instantaneous rate at which a particular reaction can occur when the system is in state $𝒙$.

Define the total propensity

$a_0(𝒙)=\sum _j{a}_j(𝒙).$

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

$p(\tau ,j\mid 𝒙,t)=a_j(𝒙)\exp [-a_0(𝒙)\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(𝒙)}\ln \left(\frac{1}{r_1}\right).$

The next reaction $j$ is the smallest integer satisfying

$\sum _{j^{\prime }=1}^j{a}_{j^{\prime }}(𝒙)>r_2{a}_0(𝒙).$

Direct method

The direct form of the Gillespie algorithm can be written as follows:

1. Initialize the time $t=t_0$ and state $𝒙={𝒙}_0$.
2. Evaluate each propensity $a_j(𝒙)$ and their total $a_0(𝒙)$.
3. Generate the waiting time $\tau$ and select the next reaction $j$.
4. Advance the time according to $t\leftarrow t+\tau$.
5. Update the state according to $𝒙\leftarrow 𝒙+{𝝂}_j$, where ${𝝂}_j$ is the state-change vector associated with reaction $j$.
6. 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.

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.

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. 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.

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

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.

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

Partial-propensity techniques have additionally been extended to reaction networks containing delays.

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. 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.

Example

Reversible binding of A and B to form AB dimers

A simple chemical system illustrates the direct method. Consider molecules of two types, Template:Math and Template:Math, which reversibly bind to form Template:Math dimers:

$\mathrm{A}+\mathrm{B}\stackrel{-\rightharpoonup }{\leftharpoondown -}\mathrm{A}\mathrm{B}$

There are two possible reactions. An Template:Math molecule and a Template:Math molecule can combine to produce an Template:Math 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 Template:Math, $n_{\mathrm{B}}$ molecules of Template:Math, and $n_{\mathrm{A}\mathrm{B}}$ dimers, the two propensities are

$a_1=k_{\mathrm{D}}{n}_{\mathrm{A}}{n}_{\mathrm{B}}$

and

$a_2=k_{\mathrm{B}}{n}_{\mathrm{A}\mathrm{B}}.$

The total reaction propensity is therefore

$a_0=k_{\mathrm{D}}{n}_{\mathrm{A}}{n}_{\mathrm{B}}+k_{\mathrm{B}}{n}_{\mathrm{A}\mathrm{B}}.$

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}\mathrm{+}\mathrm{B}\mathrm{\to }\mathrm{A}\mathrm{B})=\frac{k_{\mathrm{D}}{n}_{\mathrm{A}}{n}_{\mathrm{B}}}{a_0},$

while the probability of dissociation is

$P(\mathrm{A}\mathrm{B}\mathrm{\to }\mathrm{A}\mathrm{+}\mathrm{B})=\frac{k_{\mathrm{B}}{n}_{\mathrm{A}\mathrm{B}}}{a_0}.$

If dimer formation occurs, $n_{\mathrm{A}}$ and $n_{\mathrm{B}}$ each decrease by one and $n_{\mathrm{A}\mathrm{B}}$ 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.

See also

• Continuous-time Markov chain
• Kinetic Monte Carlo
• Master equation
• Monte Carlo method
• Stochastic process
• Tau-leaping

Further reading

• Template:Cite book

• Template:Cite book

• Template:Cite journal