Gillespie algorithm

The Gillespie algorithm generates a statistically correct trajectory (possible solution) of a stochastic equation. It was developed and published by Dan Gillespie in 1977 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 typically number in the tens of molecules (or less). Mathematically, it is a variety of a dynamic Monte Carlo method and similar to the kinetic Monte Carlo methods. It is used heavily in computational systems biology.

Additional recommended knowledge

How to ensure accurate weighing results every day?

How to quickly check pipettes?

What is the Sensitivity of my Balance?

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. When simulated, a Gillespie realization represents a random walk that exactly represents the distribution 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

Below is a summary of the steps to run the Gillespie algorithm (math omitted for now):

• 0) Initialization: Initialize the number of molecules in the system, reactions constants, and random number generators.
• 1) Monte Carlo Step: Generate random numbers to determine the next reaction to occur as well as the time interval.
• 2) Update: Increase the time step by the randomly generated time in Step 1. Update the molecule count based on the reaction that occurred.
• 3) Iterate: Go back to Step 1 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 Gibson & Bruck 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. See the articles under further reading for details.

Further reading

• Summary of Gillespie Algorithm with Matlab examples
• Interactive Evaluation of the Gillespie Algorithm for Protein Synthesis
• Daniel T. Gillespie, "Exact Stochastic Simulation of Coupled Chemical Reactions". The Journal of Physical Chemistry, Vol. 81, No. 25, pp. 2340-2361 (1977).
• Daniel T. Gillespie, "A General Method for Numerically Simulating the Stochastic Time Evolution of Coupled Chemical Reactions". Journal of Computational Physics 22(4), 403-434 (1976).
• M. Rathinam, L. R. Petzold, Y. Cao, and Daniel T. Gillespie, "Stiffness in stochastic chemically reacting systems: The implicit tau-leaping method". Journal of Chemical Physics, Volume 119, Issue 24, pp. 12784-12794 (2003).
• M. A. Gibson and J. Bruck, "Efficient Exact Stochastic Simulation of Chemical Systems with Many Species and Many Channels". J. Phys. Chem. A 104: 1876-1889 (2000).