These two papers describe a simple and very popular algorithm to exactly simulate a jump Markov process along with some interesting examples. This type of simulation has been referred to as "stochastic simulation", "kinetic Monte Carlo", or "continuous time Monte Carlo". While there were some earlier algorithms that did similar things (N-fold, etc), the SSA did not approximate the waiting times between events using some fixed time step and is thus considered the "exact" solution. Of course, Monte Carlo itself has a long history going back to the 50s, but the SSA is the first rigorous algorithm to simulate non-equilibrium systems described using a Master equation.

Daniel T. Gillespie
Journal of Computational Physics, v22, (1976)
Journal of Physical Chemistry, v81 n25 (1977)

Only 23 (!) years later, interest in continuous time Monte Carlo simulations of chemical networks renewed. This paper describes two useful data structures to improve the efficiency of the original SSA and the usage of the gamma-distribution in representing non-Markovian, multistep processes. There are other variants to the SSA, but the Next Reaction variant has good scaling properties and does not require pre-computation tricks.

Michael Gibson, Joshua Bruck
Journal of Physical Chemistry A, v104, (2000)

If a system contains one or more frequently occurring reactions then the SSA or any of its variants will require significant computational time. To speed up the simulation, one may approximate any reaction that is both frequently occurring and acting on a large number of reactant or product molecules as a continuous Markov process. The dynamics of these "fast/continuous" reactions may then be described using the chemical Langevin equation, a type of stochastic differential equation (SDE). However, one must still combine the actions of the remaining "slow/discrete" reactions with the "fast/continuous" ones. This system is classified as a coupled hybrid jump/continuous Markov process and the coupling between the discrete and continuous processes must be explicitly accounted for.

This method generates the dynamics of a hybrid jump/continuous Markov process by converting the entire problem into the simultaneous solution of two different types of SDEs. The chemical Langevin equation is a system of Ito-type SDEs describing the rates of change of the number of molecules of each species affected by the "fast/continuous" reactions. To describe the times when "slow/discrete" reaction events occur, we derive a system of differential Jump equations, which are Jump-type SDEs, and are coupled to the chemical Langevin equation. By simultaneously solving these two different types of SDEs, we can explicitly include the coupling between the jump and continuous Markov processes without resorting to hand-waiving. We may also take advantage of existing SDE numerical integrators and use current theory to quantify the numerical error. Of course, most importantly, if the system contains one or more fast reactions the method speeds up the stochastic simulation by orders of magnitude.

Howard Salis, Yiannis Kaznessis
Journal of Chemical Physics, v122 n5 (2005)

This is the current Bible for the numerical solution of SDEs. It includes basic probability, the theory behind each numerical method (including a good treatment of stochastic integrals), algorithmic descriptions of each method with pictures, and some applications/examples. It does not cover adaptive time step methods for SDEs, however, which are relatively new.

Peter Kloeden, Eckhard Platen
Springer, 2nd edition (2000)

A detailed description of the GUI, the simulation programs, and the numerical methods behind each program. There are some interesting examples, including a biochemical network with multiple, disparate timescales that spontaneously transitions between two stable states. We also include a very large system benchmark with up to 20000 reactions and 50000 species!

Howard Salis, Vassilios Sotiropoulos, Yiannis Kaznessis
BMC Bioinformatics (in press for 2006)

This method is an additional, complementary way to speed up the stochastic simulation of a system of bio/chemical reactions. If a reaction has a high rate, but acts on a small number of product or reactant molecules (called "fast/discrete") then it may not be validly approximated as a continuous Markov process. Instead, one may dynamically determine if a subset of reactions in the network has reached a "probabilistic steady state". A probabilistic steady state occurs when the joint distribuion of a subset of species (or, equivalently, a subset of reaction extents) is constant or slowly changing for a long period of time. By dynamically determining when a quasi steady-state occurs and sampling from this distribution, the method very accurately predicts the future dynamics of both the "fast/discrete" and "slow/discrete" reactions without executing the numerous extraneous fast/discrete reaction events. The method works on arbitrary bio/chemical systems, including ones with highly non-linear kinetics. Depending on the separation in time scales, the method may speed up the original stochastic simulation algorithm by many orders of magnitude and removes the "dynamical stiffness" from the system. This method is not yet integrated with the hybrid jump/continuous Markov process simulator, but will be sometime in the future. Its stand alone code is available on request. Stay tuned for further developments!

Howard Salis, Yiannis Kaznessis
Journal of Chemical Physics, v123 n21 (2005)

The quality of a bistable switch gene network (similar to the prototype built by Gardner, Cantor, and Collins in 2000) is studied using a systematic computational approach. The model is a highly detailed one using components from the lac operon. Design parameters, including the number and affinity of repressor-binding operators, the extent of DNA looping, and the strength of auto negative feedback, are varied and the "certainty" and response time of the bistable switch are measured. Design principles are proposed that summarize the results and enable greater efficiency in building gene networks that behave as bistable switches. The simulations use an early form of Hy3S.

Howard Salis and Yiannis Kaznessis
Computers & Chemical Engineering, Vol 29(3), 2005

A systematic computational study of an oscillating gene network composed from components of the lac, tet, and ara operons in E. coli (similar in connectivity to the repressilator built by Elowitz and Leibler in 2000). The design parameters of this gene network, including the number and affinity of repressor-binding operators and the half-lives of mRNAs and proteins, are varied and the average and standard deviation of the period of oscillations are measured by assuming that the oscillations are an almost cyclostationary signal. Design principles are proposed to summarize the best ways to construct a robust oscillating gene network.

Lisa Tuttle, Howard Salis, Jonathan Tomshine, and Yiannis Kaznessis
Biophysical Journal, v89 (2005)