Another approach that resolves this accuracy problem was proposed by Daniel T. Gillespie in the mid 1970s. His algorithm, called the First Reaction Method, treats the systems as a stochastic process with discrete variables that represent the populations, not concentrations, of chemical species.

To understand this algorithm fully, we must first define the term propensity used in step 5. Consider the following chemical reaction:

A + B --> C (1)

The probability that the reaction given in equation 1 occurs, or the probability that a given molecule A reacts with a given molecule B , in a small time dt is

P1 = a1dt + o(dt) (2)

As dt approaches zero, the propensity term 1 a dominates equation 2. The propensity may be a function of volume, temperature, concentration, etc. In step 4 of the First Reaction Method we calculate the propensity of the reaction based on the stochastic rate constant associated with the reaction 1 k and the current populations of the reactants A X and B X , using the following equation:

a1 = XA . AB . k1 (3)

Multiplying XA and XB together in equation 3 reflects the number of combinations by which the reaction could occur, thus making the propensity depend on the concentrations of the chemical reactants. The input rate constant 1 k is used to account for all other factors (volume, temperature, etc.) that may determine the propensity of the reaction.

Step 5 of the algorithm uses the propensity generated in step 4 to generate a putative time or the amount of time it will take for this reaction to occur based on the current state of the system. This is accomplished by generating a uniformly distributed random number between 0 and 1 (URN), scaling it to fit the exponential distribution, and multiplying that number by the propensity, as shown in the equation below.

t = (-1/a)log(URN) (4)

Step 6 selects the reaction with the smallest putative time using a linear search. Step 7 updates the number of molecules based on the stoichiometry of the reaction. For example, if we were executing the reaction given in equation 1, we would decrement the values XA and XB and increment the value of C X . Step 8 updates the current time based on the putative time selected in Step 6.