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 A recent review (Gillespie, 2007) outlines three different, but equivalent formulations; the direct, first-reaction, and first-family methods, whereby the former two are special cases of the latter. The formulation of the direct and first-reaction methods is centered on performing the usual Monte-Carlo inversion steps on the so-called "fundamental premise of stochastic chemical kinetics", which mathematically is the function , where each of the terms are propensity functions of an elementary reaction, whose argument is , the vector of species counts. The parameter is the time to the next reaction (or sojourn time), and is the current time. To paraphrase Gillespie, this expression is read as "the probability, given , that the system's next reaction will occur in the infinitesimal time interval , and will be of stoichiometry corresponding to the th reaction". This formulation provides a window to the direct and first-reaction methods by implying is an exponentially-distributed random variable, and is "a statistically independent integer random variable with point probabilities ". Thus, the Monte-Carlo generating method is simply to draw two pseudorandom numbers, and on , and compute , and the smallest integer satisfying . Utilizing this generating method for the sojourn time and next reaction, the direct method algorithm is stated by Gillespie as 1. Initialize the time and the system's state 2. With the system in state at time , evaluate all the and their sum 3. Effect the next reaction by replacing and 4. Record as desired. Return to step 1, or else end the simulation. This family of algorithms is computationally expensive and thus many modifications and adaptations exist, including the next reaction method (Gibson & Bruck), 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