SSA - the original stochastic simulation algorithm
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)
Next Reaction variant of the SSA (optimizations)
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)
HyJCMSS - The Hybrid Jump/Continuous Markov Stochastic Simulator
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)
Numerical Solution of Stochastic Differential Equations
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)
Multiscale Hy3S: Hybrid Stochastic Simulation for Supercomputers
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)
An Equation-Free Probabilistic Steady State Approximation
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)
Numerical Simulation of Stochastic Gene Circuits
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
Model-Driven Designs of an Oscillating Gene Network
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)