MM:NUCM
Gillespie Algorithm
Dr Rich Bingham
Gillespie Algorithm
In previous practicals, you have used Poisson processes to model birth/death processes. While these
can be useful, when multiple events overlap, the approximations used to justify them break down.
A powerful & popular alternative to the Poisson process simulation is the Gillespie algorithm, or the
stochastic simulation algorithm.
Idea
If we can calculate the probabilities of every event that can happen in a system, we can know what is
most likely to happen next. By randomly picking the next event weighted according to the
probabilities, we can reconstruct likely trajectories of the system while still allowing for variation.
System states
Sampling frequencies
Gillespie Algorithm
The algorithm is deceptively simple.
1 - Set up the system (Initialise)
Setting up the number of each constituents & setting rates & time.
2 - Calculate all probabilities
Using the description of the model we can evaluate the ‘probability’ of every process.
3 – Pick Reaction
Using the probabilities to weight the choice, randomly pick a process to run.
4 – Run picked process & iterate time (Fire)
We use the total probability to scale the time step.
5 – Repeat until you are done!
Rate reaction networks
The Gillespie algorithm is well suited to simulate ODE systems that can be discretized, such as
population models e.g. Lotka-Volterra. We do this by changing the ODEs into a set of reaction rate
equations.
[1]
[2]
[3]
Probability of [1] = Probability of [2] = Probability of [3] =
Running the algorithm
The timestep scaling ensures that your simulation happens at the correct speed. The random
reaction selection ensures proper sampling.
Step 2 (after initializing - n_t = 100, beta=0.1, etc…)
P(1) = P(2)= P(3)=
Total Prob = T_P = P(1)+P(2)+P(3)
=
T_P P(1) P(2) P(3)
Step 3
Pick a random number (r) & scale by the probability
r * T_P =
This means we are going to fire reaction 3
P(1) P(2) P(3)
Running the algorithm
The timestep scaling ensures that your simulation happens at the correct speed. The random
reaction selection ensures proper sampling.
Step 4
We are firing reaction 3
P(3)
Iterate the populations
n_t = n_t -1 n_v = n_v -1 n_i = n_i + 1
Iterate the timestep
we use a random number (r) & the probability to scale the time step tau
tau = ln (1/r)
tau = tau/T_P
time = time + tau
This repeats (2-4) until an endpoint (size or population) is reached.
Stochastic Output
Any single run is stochastic, but the average of many runs (usually) gives the same output as a
deterministic calculation.
We often are interested in the
state of the system after an
initial equilibration.
Repeated simulations show the
range of outcomes that can
happen
Give it a go yourself!
