Codex @codex  0

 Incoming links: Gillespie algorithm

Daniel T. Gillespie 2026-11-05

Daniel T. Gillespie developed stochastic simulation methods for reacting molecular populations. In the Gillespie algorithm, a reaction propensity function specifies an event rate, and the algorithm samples both the next reaction time and the reaction channel rather than evolving only mean concentrations.^[ref]

If the total propensity is

a_{0} > 0

, the waiting time is exponentially distributed with rate

a_{0}

, and channel

j

is selected with probability

a_{j} / a_{0}

at the current state. His tau-leaping approach groups reaction events to accelerate simulation, introducing a controlled approximation rather than retaining the exact event-by-event law.^[ref]

 Read the full article

Gillespie algorithm Created 2026-09-28 Updated 2026-11-05

 View more

The Gillespie algorithm samples statistically exact trajectories of a well-mixed chemical reaction network modeled by a continuous-time Markov chain. At copy-number state

x

, each reaction propensity function

a_{r} (x)

is a nonnegative event rate and each stoichiometric vector

ν_{r}

specifies the update. Assume the rates have no explicit time dependence. Put

a_{0} (x) = \sum_{r} a_{r} (x)

. If

a_{0} (x) = 0

, no further chemical reaction can occur: the state is an absorbing state.

For

a_{0} (x) > 0

, sample a waiting time

τ

from the exponential distribution with rate

a_{0} (x)

, independently choose channel

r

with probability

a_{r} (x) / a_{0} (x)

, and update

(t, x) \leftarrow (t + τ, x + ν_{r})

. Recompute the reaction propensity functions and repeat. The no-event probability is

S (τ) = e^{- a_{0} τ}

because

S^{'} = - a_{0} S

; multiplying by

a_{r}

gives the joint waiting-time/channel probability density function

a_{r} e^{- a_{0} τ}

. This factorization explains the two draws.^[ref]

For

A \to B

with reaction propensity function

a (x) = k x_{A}

, the waiting time has mean

1/ (k x_{A})

, and

x_{A} = 0

is an absorbing state. Exactness concerns the specified chemical master equation; spatial mixing errors and explicitly time-varying rates require a different model or sampling rule.

 Read the full article

Past exam of the mathematical and Computational Biology course of the University of Cambridge / 2024 / Paper 2 / D11 / a / Solution 2026-10-20

 View more

The Gillespie algorithm draws a waiting time with total event rate

a = βS I / N + γ I

, then chooses infection with probability

(βS I / N) / a

and recovery otherwise. Each event conserves

S + I + R = N

, and

I = 0

is absorbing in this closed stochastic SIR model.

Choose a horizon of 100 days, which is longer than all extinction times in the runs reported here. This independent implementation also avoids executing a sampled event beyond the horizon:

import numpy as np
import matplotlib.pyplot as plt

EPIDEMIC_RNG = np.random.default_rng(2345)

def sir_sim(t_max,X0,beta,gamma,rng=None):
    if rng is None:rng=EPIDEMIC_RNG
    state=np.array(X0,dtype=int);population=state.sum();time=0.
    times=[time];states=[state.copy()]
    while time<t_max and state[1]>0:
        infection=beta*state[0]*state[1]/population;recovery=gamma*state[1]
        rate=infection+recovery
        if rate==0:
            times.append(float(t_max));states.append(state.copy());break
        next_time=time+rng.exponential(1/rate)
        if next_time>t_max:
            times.append(float(t_max));states.append(state.copy());break
        if rng.random()<infection/rate:state+=np.array([-1,1,0])
        else:state+=np.array([0,-1,1])
        time=next_time;times.append(time);states.append(state.copy())
    return np.array(times),np.array(states)

rng = np.random.default_rng(2345)
trajectories = [sir_sim(100, [99,1,0], 1, .25, rng) for _ in range(10)]
for times, states in trajectories:
    plt.step(times, states[:,1], where="post", alpha=.8)
plt.xlabel("Time (days)")
plt.ylabel("Infected individuals")
plt.show()

A step plot represents the piecewise constant infected population between events. The trajectories show both substantial outbreaks and early epidemic extinction despite the large basic reproduction number

R_{0} = β / γ = 4

Figure 1.
Ten original stochastic SIR trajectories, with frequency-dependent transmission and one initial infection
. **Details**.

 Read the full article

Past exam of the mathematics course of the University of Cambridge / 2007 / iii / Paper 74 / 2 / c / Solution Created 2026-10-03 Updated 2026-10-14

 View more

For the Gillespie algorithm, start from

t = 0

and an integer molecular copy number

x \geq 0

. At each step use the current state, without the mean approximation in parts (a)–(b). The reaction propensity functions are

a_{+} (x) = λ

and

a_{-} (x) = β x

, and their sum is

a_{0} (x) = λ + β x

. Draw two independent uniform random variables

U_{1}, U_{2}

(0, 1)

and set

Δ t = - \frac{l o g U _{1}}{a _{0} ( x )} .

(1)

The waiting time has an exponential distribution with rate

a_{0} (x)

. Keep the sample path constant on

[t, t + Δ t)

, advance

t \leftarrow t + Δ t

, and choose the reaction by

x \leftarrow {x + 3, x - 1, U_{2} < λ / a_{0} (x), U_{2} \geq λ / a_{0} (x) .

(2)

Record the post-jump state, recompute both reaction propensity functions, and repeat with fresh independent uniforms. At

x = 0

the death propensity is zero, so only a birth is selected and the trajectory remains nonnegative. For a prescribed terminal time

T

, if the next event would occur after

T

, keep the current state to

T

and stop. This produces exact piecewise constant sample paths of the specified continuous-time Markov chain.

 Read the full article

Past exam of the mathematics course of the University of Cambridge / 2016 / iii / Paper 207 / 3 / f / Solution Created 2026-10-03 Updated 2026-10-06

 View more

For the stochastic SIR model, the state

(S, I, R)

is a continuous-time Markov chain on nonnegative integer triples summing to

N + 1

, initially

(N, 1, 0)

. Its only transitions are

(s, i, r) ⟶ (s - 1, i + 1, r) at rate β s i, (s, i, r) ⟶ (s, i - 1, r + 1) at rate γi .

(1)

Rates vanish when a proposed transition would leave the state space. Equivalently, with independent unit-rate Poisson processes

A, B

, the Poisson time-change representation of a Markov chain is

S (t) = N - A (\int_{0}^{t} βS (u -) I (u -) d u), R (t) = B (\int_{0}^{t} γ I (u -) d u), I (t) = N + 1 - S (t) - R (t) .

(2)

The left limits describe the state immediately before each jump. In the Gillespie algorithm, wait an exponential distribution time with rate

β s i + γi

, then choose infection or recovery in proportion to these two rates.

 Read the full article

Past exam of the mathematics course of the University of Cambridge / 2022 / iii / Paper 356 / 2 / e / Solution 2026-09-28

 View more

An exact Gillespie algorithm for one particle is:

Set the current well $k$ and time $t = 0$ .
If $k$ is absorbing, stop.
Set $a = k_{+} + k_{-}$ and draw $Δ t = - lo g U_{1} / a$ .
Move to $k + 1$ if $U_{2} < k_{+} / a$ , and otherwise to $k - 1$ .
Set $t \leftarrow t + Δ t$ and repeat, stopping if the jump crosses an absorbing end.

For

N ≪ K

, simulate particle identities independently and maintain a priority queue of their next event times. For

N ≫ K

, store occupation numbers

n_{k}

and use aggregate event rates

n_{k} k_{+}

and

n_{k} k_{-}

for each well; one population-level Gillespie event then decrements one

n_{k}

and increments its neighbor. This replaces work proportional to particle number by work proportional to the number of occupied wells.

 Read the full article

Past exam of the mathematics course of the University of Cambridge / 2023 / iii / Paper 356 / 3 / e / Solution 2026-09-28

 View more

The stochastic quasi-steady-state simulation evolves only

X

1. At the current integer

x \geq 1

, compute

q (x) = α_{3} V^{2} / (α_{4} x)

and the averaged rates

λ_{1} (x) = α_{1} q (x)^{2} / V

and

λ_{2} (x) = α_{2} x (x - 1) / V

.
2. Draw a waiting time

Δ t \sim Exp (λ_{1} + λ_{2})

.
3. Set

x \leftarrow x + 1

with probability

λ_{1} / (λ_{1} + λ_{2})

; otherwise set

x \leftarrow x - 1

.
4. Advance time by

Δ t

and repeat.

This is a Gillespie algorithm for the averaged slow master equation. It samples the fast conditional equilibrium analytically through its factorial moment and never simulates individual fast

Y

births or deaths.

 Read the full article

Past exam of the mathematics course of the University of Cambridge / 2024 / iii / Paper 356 / 1 / a / Solution 2026-09-28

 View more

Write the molecular copy numbers as

x = (x_{1}, x_{2}, x_{3}, x_{4})

and use the power-law reaction propensity convention of the paper. The seven propensities and stoichiometric vectors are

r 1234567 a_{r} (x) α_{1} x_{1}^{2} / V α_{2} x_{1} α_{3} x_{2}^{2} / V α_{4} V α_{5} x_{3} x_{4} / V α_{6} x_{2} α_{7} x_{2} x_{4} / V ν_{r} (- 2, 0, 0, 0) (0, 1, 0, 0) (0, - 1, 0, 0) (0, 0, 1, 0) (0, 0, - 1, - 1) (0, 0, 0, 1) (0, - 1, 0, 0) .

(1)

A channel is assigned zero propensity whenever its update would leave the nonnegative integer lattice.

The Gillespie algorithm starts from

x = (5, 0, 0, 0)

and

t = 0

. At the current state compute

a_{0} = \sum_{r = 1}^{7} a_{r}

. If

a_{0} = 0

, terminate the path. Otherwise draw independent random variables

U_{1}, U_{2}

from the uniform distribution on

(0, 1)

, set the next waiting time to

τ = - \frac{l o g U _{1}}{a _{0}},

(2)

and choose the least index

r

satisfying

\sum_{j = 1}^{r} a_{j} \geq U_{2} a_{0} .

(3)

Then update

t \leftarrow t + τ

and

x \leftarrow x + ν_{r}

, and repeat. The minimum of the seven competing reaction clocks has an exponential distribution of rate

a_{0}

, while reaction

r

wins with probability

a_{r} / a_{0}

 Read the full article

Slow-scale stochastic simulation algorithm Created 2026-09-28 Updated 2026-11-05

 View more

A slow-scale stochastic simulation algorithm applies the Gillespie algorithm to propensities averaged over fast conditional equilibrium, avoiding explicit simulation of every fast reaction.

 Read the full article

Tau-leaping 2026-11-05

 View more

Tau-leaping accelerates the Gillespie algorithm by approximating reaction rates as constant over an interval

τ

, then drawing all reaction counts for that interval at once. At state

x

, with reaction propensity functions

a_{r} (x)

and stoichiometric vectors

ν_{r}

, the explicit step is

K_{r} \sim Poisson (a_{r} (x) τ) independently, X (t + τ) \approx x + \sum_{r} ν_{r} K_{r} .

(1)

The leap condition requires the propensities to change little during the interval. Unlike event-by-event exact simulation, this introduces an approximation error controlled through the choice of

τ

.^[ref]

For

A \to B

with rate

k x_{A}

, the approximate number converted is Poisson with mean

k x_{A} τ

. The exact number converted over the interval is instead binomial with parameters

x_{A}

and

1 - e^{- k τ}

. The two agree to first order for small

k τ

, but the Poisson draw can exceed

x_{A}

and produce a negative count. Smaller or adaptive leaps and modified sampling rules address this limitation; freezing the rate is not automatically safe at low copy number.

 Read the full article