Gillespie Stochastic Simulation Algorithm

The kinetics of chemical reaction networks are often modeled using differential equations, treating the amounts of chemical species as continuous and deterministic quantities. If the number of reacting molecules is small, random molecular fluctuations become significant and a discrete, stochastic modeling approach is often more appropriate. Under assumptions of a "well-mixed" environment, such a system can be modeled as a continuous-time Markov process, whose probability distribution obeys what is called a chemical "master equation". The Gillespie algorithm (or SSA) is a discrete-event simulation algorithm that produces single realizations of the stochastic process that are in exact statistical agreement with the master equation.
This submission includes simple implementations of the two original versions of the SSA (Direct and First-Reaction Method). The user provides a model consisting of a matrix of reaction stoichiometries, a handle to a function that computes reaction propensities (probability/unit time), an initial state vector, and initial/final times. The ssa functions returns the sequence of event times and species amounts.