Gillespie algorithm


In probability theory, the Gillespie algorithm generates a statistically correct trajectory of a stochastic equation. It was created by Joseph L. Doob and others, 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. 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 is low and keeping track of the position and behaviour of individual molecules is computationally feasible. Mathematically, it is a variant of a dynamic Monte Carlo method and similar to the kinetic Monte Carlo methods. It is used heavily in computational systems biology.

History

The process that led to the algorithm recognizes several important steps. In 1931, Andrei Kolmogorov introduced the differential equations corresponding to the time-evolution of stochastic processes that proceed by jumps, today known as Kolmogorov equations . It was William Feller, in 1940, who found the conditions under which the Kolmogorov equations admitted probabilities as solutions. In his Theorem I 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 extended Feller's solutions beyond the case of pure-jump processes. The method was implemented in computers by David George Kendall using the Manchester Mark 1 computer and later used by Maurice S. Bartlett in his studies of epidemics outbreaks. Gillespie obtains the algorithm in a different manner by making use of a physical argument.

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. 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. 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 :
  1. Initialization: Initialize the number of molecules in the system, reaction constants, and random number generators.
  2. Monte Carlo step: Generate random numbers to determine the next reaction to occur as well as the time interval. The probability of a given reaction to be chosen is proportional to the number of substrate molecules, the time interval is exponentially distributed with mean.
  3. Update: Increase the time by the randomly generated time in Step 2. Update the molecule count based on the reaction that occurred.
  4. Iterate: Go back to Step 2 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, 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. In weakly coupled networks, the number of reactions that is influenced by any other reaction is bounded by a small constant. In strongly coupled networks, a single reaction firing can in principle affect all other reactions. 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. 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. See the articles cited below for details.
Partial-propensity formulations, as developed independently by both Ramaswamy et al. and Indurkhya and Beal, 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 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. The use of partial-propensity methods is limited to elementary chemical reactions, i.e., reactions with at most two different reactants. Every non-elementary chemical reaction can be equivalently decomposed into a set of elementary ones, at the expense of a linear increase in network size.

Simple example: Reversible binding of A and B to form AB dimers

A simple example may help to explain how the Gillespie algorithm works. Consider a system of molecules of two types, and. In this system, and reversibly bind together to form dimers such that two reactions are possible: either A and B react reversibly to form an dimer, or an dimer dissociates into and. The reaction rate constant for a given single A molecule reacting with a given single molecule is, and the reaction rate for an dimer breaking up is.
If at time t there is one molecule of each type then the rate of dimer formation is, while if there are molecules of type and molecules of type, the rate of dimer formation is. If there are dimers then the rate of dimer dissociation is.
The total reaction rate,, at time t is then given by
So, we have now described a simple model with two reactions. This definition is independent of the Gillespie algorithm. We will now describe how to apply the Gillespie algorithm to this system.
In the algorithm, we advance forward in time in two steps: calculating the time to the next reaction, and determining which of the possible reactions the next reaction is. Reactions are assumed to be completely random, so if the reaction rate at a time t is, then the time, δt, until the next reaction occurs is a random number drawn from exponential distribution function with mean. Thus, we advance time from t to t + δt.
The probability that this reaction is an molecule binding to a molecule is simply the fraction of total rate due to this type of reaction, i.e.,
the probability that reaction is
The probability that the next reaction is an dimer dissociating is just 1 minus that. So with these two probabilities we either form a dimer by reducing and by one, and increase by one, or we dissociate a dimer and increase and by one and decrease by one.
Now we have both advanced time to t + δt, and performed a single reaction. The Gillespie algorithm just repeats these two steps as many times as needed to simulate the system for however long we want. The result of a Gillespie simulation that starts with and at t=0, and where and, is shown at the right. For these parameter values, on average there are 8 dimers and 2 of and 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.
That was just a simple example, with two reactions. More complex systems with more reactions are handled in the same way. All reaction rates must be calculated at each time step, and one chosen with probability equal to its fractional contribution to the rate. Time is then advanced as in this example.

Another example: The SIR epidemic without vital dynamics

The SIR model is a classic biological description of how certain diseases permeate through a fixed-size population. In its simplest form there are members of the population, whereby each member may be in one of three states -- susceptible, infected, or recovered -- at any instant in time, and each such member transitions irreversibly through these states according to the directed graph below. We can denote the number of susceptible members as, the number of infected members as, and the number of recovered members as. Therefore we may also conclude that for any point in time.
Further, a given susceptible member will transition to the infected state by coming into contact with any of the infected members, and so infection occurs with rate . A given member of the infected state recovers without dependence on any of the three states, which is specified by rate . Given this basic scheme, it possible to construct the following non-linear system.
This system has no analytical solution. However, with the Gillespie algorithm, it can be simulated many times, and a regression technique such as least-squares may be applied to fit a polynomial over all of the trajectories. As the number of trajectories increases, such polynomial regression will asymptotically behave like an analytic solution. In addition to estimating the solution to an intractable problem like the SIR epidemic, the stochastic nature of each trajectory allows one to compute statistics other than.
The trajectory presented in the above figure was simulated with the following Python implementation of the Gillespie algorithm.

import math
import random
  1. Input parameters ####################
  2. int; total population
N = 350
  1. float; maximum elapsed time
T = 100.0
  1. float; start time
t = 0.0
  1. float; spatial parameter
V = 100.0
  1. float; rate of infection after contact
_alpha = 10.0
  1. float; rate of cure
_beta = 0.5
  1. int; initial infected population
n_I = 1
  1. ########################################
  2. Compute susceptible population, set recovered to zero
n_S = N - n_I
n_R = 0
  1. Initialize results list
SIR_data =
SIR_data.append)
  1. Main loop
while t < T:
if n_I 0:
break
w1 = _alpha * n_S * n_I / V
w2 = _beta * n_I
W = w1 + w2
dt = -math.log / W
t = t + dt
if random.uniform < w1 / W:
n_S = n_S - 1
n_I = n_I + 1
else:
n_I = n_I - 1
n_R = n_R + 1
SIR_data.append)
with open as fp:
fp.write