Gillespie’s direct algorithm step by step

Gillespie’s direct algorithm repeatedly combines exact exponential waiting times with rate-weighted event selection to generate one complete continuous-time Markov-chain trajectory.

What the algorithm does

Gillespie’s direct algorithm, also called the stochastic simulation algorithm, generates the next CTMC jump directly from the current event rates. It advances from event to event rather than through small fixed time steps.

\[ (t_0,X_0)\longrightarrow (t_1,X_1)\longrightarrow (t_2,X_2)\longrightarrow\cdots, \]

The jump times satisfy \(t_0<t_1<t_2<\cdots\), but the gaps are random and generally unequal.

The complete algorithm

1. Current stateRead time and compartment counts.

2. Event ratesCalculate every currently possible biological rate.

3. Total rateAdd the event rates.

4. Waiting timeGenerate \(\tau=-\log U_1/a_0\).

5. Event typeSelect using \(U_2\) and rate proportions.

6. Jump and repeatAdvance time, update state, store, and return to step 2.

The rates are recalculated after every event because they depend on the newly updated state.

Why recalculation is essential

Suppose an infection occurs. Then \(S\) decreases and \(E\) increases. Consequently, both the future infection rate and progression rate may change.

Do not calculate the rates once outside the loop and reuse them. That would simulate a different process whose event intensities ignore the changing epidemic state.

No fixed time step

DTMC or Euler calculation

Time advances by a chosen fixed increment such as dt = 0.1.

Gillespie CTMC

Time advances by a newly generated positive waiting time \(\tau\) after every state.

time = time + waiting_time

The code therefore has no numerical time-step parameter. The observation horizon limits how long to simulate, but it does not set jump times.

Why use a `while` loop?

The number of events is not known before simulation. A while loop is appropriate because repetition continues while biological and time conditions remain true:

while time < max_time:
    # calculate and execute the next event

The loop may stop after very few events because of early extinction or after many events in a large outbreak.

Stopping conditions

Absorbing stateAll event rates are zero, so no future jump is possible.

Biological extinctionFor basic SEIR, \(E=0\) and \(I=0\).

Observation horizonThe next generated jump would occur after max_time.

The total-rate test is the most general absorbing-state check. The \(E=I=0\) test gives its biological meaning for this SEIR model.

Do not execute an event beyond the horizon

After generating \(\tau\), first check:

if time + waiting_time > max_time:
    stopping_reason = "time horizon reached"
    break

If the event lies outside the observation window, do not update the state and do not store it as an observed event.

The state remains unchanged from the last jump until the horizon. A plotting copy can be extended to max_time without pretending that another biological event occurred.

Store two related datasets

State path

Store initial time and state, then every post-event time and state. This produces the step plot.

Event log

Store one row per actual jump: event number, waiting time, event type, new time and post-event state.

The initial state is part of the path but is not an event. Label it separately rather than counting it as infection, progression or recovery.

Scenario

Use a closed SEIR population:

\[ (S_0,E_0,I_0,R_0)=(195,3,2,0), \qquad N=200, \] \[ \beta=0.3,\qquad \sigma=0.2,\qquad \gamma=0.1\text{ per day}. \]

The observation horizon is 200 days. The simulation stops earlier if both exposed and infectious populations reach zero.

Interactive Python laboratory

Click Run code to generate one complete event-driven trajectory. The output includes a summary, first and last events, event counts, a compartment step plot and the sequence of random waiting times.

Interactive PythonComplete Gillespie SEIR trajectory

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

def seir_rates(S, E, I, N, beta, sigma, gamma):
    return np.array([
        beta * S * I / N,
        sigma * E,
        gamma * I
    ])

# Model inputs
S, E, I, R = 195, 3, 2, 0
N = S + E + I + R

beta = 0.30
sigma = 0.20
gamma = 0.10

time = 0.0
max_time = 200.0
seed = 12
rng = np.random.default_rng(seed)

# State-path storage includes the initial state
times = [time]
S_values = [S]
E_values = [E]
I_values = [I]
R_values = [R]

# One row will be added per actual event
event_rows = []
event_number = 0
stopping_reason = "time horizon reached"

while time < max_time:
    # 1. Recalculate rates from the current state
    rates = seir_rates(
        S, E, I, N, beta, sigma, gamma
    )

if (rates < 0).any():
        raise ValueError("An event rate became negative.")

total_rate = rates.sum()

# 2. Stop if the current state is absorbing
    if total_rate == 0:
        stopping_reason = "absorbing state reached"
        break

# 3. Generate the next waiting time
    u_time = 1.0 - rng.random()
    waiting_time = -np.log(u_time) / total_rate

# 4. Do not execute an event beyond the horizon
    if time + waiting_time > max_time:
        stopping_reason = "time horizon reached"
        break

# 5. Select the event using an independent draw
    probabilities = rates / total_rate
    cumulative = np.cumsum(probabilities)
    u_event = rng.random()

if u_event < cumulative[0]:
        event = "infection"
        S -= 1
        E += 1
    elif u_event < cumulative[1]:
        event = "progression"
        E -= 1
        I += 1
    else:
        event = "recovery"
        I -= 1
        R += 1

# 6. Advance time only when the event occurs
    time += waiting_time
    event_number += 1

# 7. Validate the updated state
    if min(S, E, I, R) < 0:
        raise ValueError("A compartment became negative.")

if S + E + I + R != N:
        raise ValueError("Population conservation failed.")

# 8. Store the post-event state and event record
    times.append(time)
    S_values.append(S)
    E_values.append(E)
    I_values.append(I)
    R_values.append(R)

event_rows.append({
        "event number": event_number,
        "time": time,
        "waiting time": waiting_time,
        "event": event,
        "S": S,
        "E": E,
        "I": I,
        "R": R
    })

# 9. Biological interpretation of absorption
    if E == 0 and I == 0:
        stopping_reason = "infection became extinct"
        break

trajectory = pd.DataFrame({
    "time": times,
    "S": S_values,
    "E": E_values,
    "I": I_values,
    "R": R_values
})

event_log = pd.DataFrame(event_rows)

# Trajectory-level summaries
peak_I = trajectory["I"].max()
peak_time = trajectory.loc[
    trajectory["I"].idxmax(), "time"
]
final_state = (
    trajectory[["S", "E", "I", "R"]]
    .iloc[-1]
    .astype(int)
    .tolist()
)

print("Gillespie trajectory summary")
print(f"Seed: {seed}")
print(f"Stopping reason: {stopping_reason}")
print(f"Final simulated event time: {time:.6f} days")
print(f"Number of events: {event_number}")
print(f"Peak infectious count: {peak_I}")
print(f"First time at peak: {peak_time:.6f} days")
print(f"Final state [S, E, I, R]: {final_state}")

if len(event_log) > 0:
    print()
    print("Event counts")
    print(event_log["event"].value_counts().to_string())
    print()
    print("First 8 events")
    print(event_log.head(8).round(6).to_string(index=False))
    print()
    print("Last 8 events")
    print(event_log.tail(8).round(6).to_string(index=False))

# Extend only the plotted state to the horizon when needed
plot_times = times.copy()
plot_S = S_values.copy()
plot_E = E_values.copy()
plot_I = I_values.copy()
plot_R = R_values.copy()

if stopping_reason == "time horizon reached" and time < max_time:
    plot_times.append(max_time)
    plot_S.append(S)
    plot_E.append(E)
    plot_I.append(I)
    plot_R.append(R)

# Figure 1: irregular-time compartment path
plt.figure(figsize=(9.5, 5.2))
plt.step(plot_times, plot_S, where="post",
         label="Susceptible", color="#2563eb")
plt.step(plot_times, plot_E, where="post",
         label="Exposed", color="#f59e0b")
plt.step(plot_times, plot_I, where="post",
         label="Infectious", color="#dc2626")
plt.step(plot_times, plot_R, where="post",
         label="Removed", color="#39815e")
plt.xlabel("Continuous time (days)")
plt.ylabel("Number of people")
plt.title("One Gillespie SEIR trajectory")
plt.legend()
plt.grid(alpha=0.25)
plt.tight_layout()

# Figure 2: waiting times in event order
if len(event_log) > 0:
    event_colours = {
        "infection": "#2563eb",
        "progression": "#f59e0b",
        "recovery": "#39815e"
    }

colours = [
        event_colours[event]
        for event in event_log["event"]
    ]

plt.figure(figsize=(9.5, 4.8))
    plt.scatter(
        event_log["event number"],
        event_log["waiting time"],
        c=colours,
        s=18,
        alpha=0.75
    )
    plt.xlabel("Event number")
    plt.ylabel("Waiting time before event (days)")
    plt.title("Random waiting times between Gillespie events")
    plt.grid(alpha=0.25)
    plt.tight_layout()

Output

Run the code to see the result.

Understand the new loop code

Code	Meaning
`while time < max_time`	Repeats events while the process remains inside the observation window.
`rates = seir_rates(...)`	Recalculates all event rates from the current post-event state.
`time += waiting_time`	Advances time by a random rather than fixed increment.
`event_rows.append({...})`	Adds one labelled record for each actual biological event.
`times.copy()`	Creates separate plotting lists so extending the plotted state does not falsify the event log.
`.value_counts()`	Counts how many infection, progression and recovery events occurred.

Why the step plot is still correct

Time is continuous, but compartment counts remain constant between jumps. A step plot shows this piecewise-constant path accurately.

The difference from the DTMC step plot is the horizontal spacing: DTMC jumps are considered at fixed grid times, whereas Gillespie jumps occur at irregular random times.

Try these experiments

Change the seed and compare the number of events, peak and extinction time.
Use \((S,E,I,R)=(197,3,0,0)\). Confirm that progression can create infectious people.
Use \((S,E,I,R)=(190,0,0,10)\). Confirm immediate absorption with zero events.
Shorten max_time and observe safe stopping before an event beyond the horizon.
Temporarily move the rate calculation outside the loop. Explain why the resulting process is biologically incorrect, then restore it.

Biological interpretation

The trajectory is one possible epidemic history. Event times cluster when the current total rate is large and become more widely separated when the total rate is small.

Infection, progression and recovery counts reflect the particular random path. Another seed can produce early extinction or a substantially different peak even with identical biological parameters.

What this lesson has—and has not—done

You can now:

implement every step of Gillespie’s direct algorithm;
recalculate rates after every event;
advance through unequal continuous waiting times;
store both a state path and event log;
stop at absorption, biological extinction or a horizon; and
plot a piecewise-constant continuous-time trajectory.

Next: Lessons 7–9 reuse this generic algorithm while focusing on the distinctive state structures and biological outcomes of SIS, SIR and SEIR CTMC models.