Gillespie method

Gillespie direct method :

For a system in a given state, Gillespie's direct algorithm asks two questions:
Which reaction occurs next?
When does it occur?
Both of these questions must be answered probabilistically.
The probability that a certain reaction μ will take place in the next instant of time dt is given by:

                    (for example)

So we can specify the probability density P( μ ,τ) that the next reaction is μ and it occurs at time τ.
It can be shown that

This equation leads directly to the answers of the two aforementioned questions.
First, what is the probability distribution for reactions? Integrating P( μ ,τ) over all τ from 0 to ∞ results in:

    *

Second, what is the probability distribution for times? Summing P( μ ,τ) over all μ results in

  **

These two distributions lead to Gillespie's direct algorithm:

 Exact Stochastic Simulation - Gillespie's  Direct Method:

1. Initialize (i.e., set initial numbers of molecules, set t = 0, set stopping criteria e.g. maximum time or minimum reactants numbers )  .

2. Calculate the propensity function, a_i, for all i.

3. Choose μ according to the distribution in eq *.

4. Choose t according to an exponential with parameter   (as in eq **).

5. Change the number of molecules to reflect execution of reaction μ . Set t = t + τ .

6. If stopping criteria are not met, go to Step 2.