Gillespie algorithm

Mathematical formulation

Consider a system whose state at time ${\displaystyle t}$ is described by a vector ${\displaystyle {\boldsymbol {x}}}$. Each possible event or reaction ${\displaystyle j}$ has a propensity function ${\displaystyle 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 ${\displaystyle {\boldsymbol {x}}}$.

Define the total propensity

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

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

${\displaystyle 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 ${\displaystyle r_{1}}$ and ${\displaystyle r_{2}}$, uniformly distributed on ${\displaystyle (0,1)}$, the waiting time can be generated as

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

The next reaction ${\displaystyle j}$ is the smallest integer satisfying

${\displaystyle \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:

1. Initialize the time ${\displaystyle t=t_{0}}$ and state ${\displaystyle {\boldsymbol {x}}={\boldsymbol {x}}_{0}}$.
2. Evaluate each propensity ${\displaystyle a_{j}({\boldsymbol {x}})}$ and their total ${\displaystyle a_{0}({\boldsymbol {x}})}$.
3. Generate the waiting time ${\displaystyle \tau }$ and select the next reaction ${\displaystyle j}$.
4. Advance the time according to ${\displaystyle t\leftarrow t+\tau }$.
5. Update the state according to ${\displaystyle {\boldsymbol {x}}\leftarrow {\boldsymbol {x}}+{\boldsymbol {\nu }}_{j}}$, where ${\displaystyle {\boldsymbol {\nu }}_{j}}$ is the state-change vector associated with reaction ${\displaystyle 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.

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.

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.

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:

${\displaystyle {\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 ${\displaystyle k_{\mathrm {D} }}$ and the rate constant for dissociation be ${\displaystyle k_{\mathrm {B} }}$. If the system contains ${\displaystyle n_{\mathrm {A} }}$ molecules of A, ${\displaystyle n_{\mathrm {B} }}$ molecules of B, and ${\displaystyle n_{\mathrm {AB} }}$ dimers, the two propensities are

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

and

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

The total reaction propensity is therefore

${\displaystyle 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 ${\displaystyle 1/a_{0}}$. Once the waiting time has been generated, the probability that the next event is dimer formation is

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

while the probability of dissociation is

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

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