Gillespie algorithm
The Gillespie stochastic simulation algorithm generates individual random trajectories of a continuous-time Markov chain (CTMC) directly from its event rates. It is widely used for epidemic, population, biochemical and other event-based biological models.
Why do we need it?
The Kolmogorov equations can describe the complete probability distribution of a CTMC. But a biological model with many compartments can have an enormous number of possible states, making the full probability system difficult to solve.
Gillespie simulation takes a different approach. It follows one possible realisation of the process, event by event.
What does “exact” mean here?
The direct Gillespie algorithm is called an exact stochastic simulation algorithm because, under the CTMC assumptions, it samples the next event time and event type from the correct distributions. It does not approximate time by fixed steps such as \(\Delta t=0.01\).
Start with a simple SIS epidemic
Let \(i\) be the current number of infectious individuals in a population of size \(N\). Two events are possible:
| Event | State change | Rate |
|---|---|---|
| infection | \(i\to i+1\) | \(b(i)=\beta(N-i)i/N\) |
| recovery | \(i\to i-1\) | \(d(i)=\gamma i\) |
The algorithm does not decide in advance that an infection or recovery happens every fixed number of minutes. Both the event time and event type are random.
Step 1 — calculate all current event rates
Suppose the possible events have rates
\[a_1,a_2,\ldots,a_m.\]Calculate the total rate
\[\boxed{a_0=a_1+a_2+\cdots+a_m}.\]For SIS,
\[\boxed{a_0=b(i)+d(i)}.\]The total rate controls how quickly some event occurs.
Step 2 — generate the waiting time
While the current state is unchanged, the total rate \(a_0\) is constant. Therefore the time \(\tau\) until the next event is exponentially distributed:
\[\tau\sim\operatorname{Exp}(a_0).\]Generate
\[U_1\sim\operatorname{Uniform}(0,1)\]and calculate
\[\boxed{\tau=-\frac{\ln U_1}{a_0}}.\]The simulation clock then advances from \(t\) to
\[t+\tau.\]Step 3 — choose which event occurs
The total rate tells us when something happens, but not what happens. Event \(k\) is selected with probability
\[\boxed{P(\text{event }k\text{ next})=\frac{a_k}{a_0}}.\]For SIS,
\[P(\text{infection next})=\frac{b(i)}{b(i)+d(i)},\]\[P(\text{recovery next})=\frac{d(i)}{b(i)+d(i)}.\]A second independent random number
\[U_2\sim\operatorname{Uniform}(0,1)\]selects between these possibilities.
Why are two random numbers used?
The two random numbers have different jobs:
| Random number | Question answered | Used for |
|---|---|---|
| \(U_1\) | When does the next event occur? | exponential waiting time |
| \(U_2\) | Which event occurs? | event selection |
Step 4 — update the biological state
If infection is selected,
\[i\leftarrow i+1.\]If recovery is selected,
\[i\leftarrow i-1.\]The simulation time becomes
\[t\leftarrow t+\tau.\]Step 5 — recalculate the rates
This step is essential. In SIS, the rates depend on the current infectious count:
\[b(i)=\beta\frac{(N-i)i}{N},\qquad d(i)=\gamma i.\]After an event changes \(i\), both rates may change. Therefore the algorithm must recalculate them before generating the next waiting time.
The complete algorithm
A complete numerical SIS step
Take
\[N=100,\qquad i=10,\qquad\beta=0.30,\qquad\gamma=0.10.\]The infection rate is
\[b(10)=0.30\frac{(100-10)(10)}{100}=2.7\text{ per day}.\]The recovery rate is
\[d(10)=0.10(10)=1.0\text{ per day}.\]Therefore
\[a_0=2.7+1.0=3.7\text{ per day}.\]Generate the event time
Suppose the first random number is
\[U_1=0.40.\]Then
\[\tau=-\frac{\ln(0.40)}{3.7}\approx0.248\text{ days}.\]So the next event occurs about 0.248 days after the current time.
Choose the event
The event probabilities are
\[P(\text{infection})=\frac{2.7}{3.7}\approx0.730,\]\[P(\text{recovery})=\frac{1.0}{3.7}\approx0.270.\]So divide the unit interval as
\[0\le U_2<0.730:\quad\text{infection},\]\[0.730\le U_2\le1:\quad\text{recovery}.\]Suppose
\[U_2=0.62.\]Then infection is selected, because \(0.62<0.730\). The new state is
\[i=11.\]Now the rates are recalculated using \(i=11\), and the procedure begins again.
How the event-selection interval works
What does a Gillespie trajectory look like?
The state remains constant while the process waits. At each random event time, the state jumps by the amount specified by that event. An SIS infectious count therefore produces a step-like trajectory.
This shape is important. The CTMC does not change smoothly between events. Nothing happens during the waiting time, then a discrete biological event changes the state.
Why is the trajectory irregular?
There are two sources of visible randomness:
| Random feature | Cause |
|---|---|
| unequal horizontal lengths | random exponential waiting times |
| upward or downward jumps | random selection between infection and recovery |
Even if two simulations start from exactly the same state with exactly the same parameters, they usually follow different trajectories.
One trajectory is not the model prediction
A single Gillespie path is only one possible realisation of the stochastic model. Another run uses different random numbers and can produce a different epidemic history.
Extinction appears naturally
In an SIS epidemic without imported infection, if
\[i=0,\]then
\[b(0)=0,\qquad d(0)=0.\]No further event can occur. The process has reached an absorbing state and the epidemic is extinct.
A deterministic SIS model may approach a positive endemic level, while a finite stochastic SIS model can still reach \(i=0\) through random fluctuations. Gillespie simulation makes this possibility visible trajectory by trajectory.
Why not just use a very small fixed time step?
A fixed-step simulation divides time into intervals \(\Delta t\) and approximates event probabilities within each interval. This can work when \(\Delta t\) is sufficiently small, but it introduces a time-discretisation approximation.
The direct Gillespie method instead jumps directly from one event time to the next.
| Small fixed-step method | Direct Gillespie method |
|---|---|
| advance by chosen \(\Delta t\) | advance by random \(\tau\) |
| approximates short-time event probabilities | samples CTMC event time exactly under model assumptions |
| may contain many steps with no event | jumps directly to the next event |
| accuracy depends on time-step choice | no arbitrary simulation time step |
Connection to the generator matrix
If the CTMC is in state \(i\), the generator contains rates \(q_{ij}\) for all possible jumps \(i\to j\). The total leaving rate is
\[a_0=-q_{ii}=\sum_{j\ne i}q_{ij}.\]Gillespie then uses
\[\tau\sim\operatorname{Exp}(a_0)\]and chooses destination \(j\) with probability
\[\frac{q_{ij}}{a_0}.\]So the Gillespie algorithm is a direct simulation procedure for the CTMC encoded by \(Q\).
Connection to Poisson and exponential ideas
The earlier sections now fit together:
The Poisson-process idea explains random event occurrence at a constant rate while the state is fixed. The exponential distribution gives the next-event waiting time. Gillespie repeatedly applies these ideas while recalculating the rates after each state change.
General reaction or compartment models
The method is not limited to SIS. Suppose the state is a vector
\[\mathbf X=(X_1,X_2,\ldots,X_r)\]and event \(k\) changes the state by a vector \(\boldsymbol\nu_k\). If its current rate is \(a_k(\mathbf X)\), then after selecting event \(k\),
\[\boxed{\mathbf X\leftarrow\mathbf X+\boldsymbol\nu_k}.\]This notation allows the same algorithm to handle SIR, SEIR, predator–prey systems, biochemical reaction networks and many other compartment models.
SIR example of event vectors
For state
\[\mathbf X=(S,I,R),\]an infection event
\[S+I\to2I\]changes the compartment counts by
\[\boldsymbol\nu_1=(-1,+1,0),\]while a recovery event
\[I\to R\]has
\[\boldsymbol\nu_2=(0,-1,+1).\]Once the corresponding event rates are defined, the same Gillespie steps apply without changing the underlying algorithm.
When the direct method can become expensive
If event rates are extremely large, the waiting times can become very short and the algorithm may need to simulate a huge number of individual events. Larger systems can therefore be computationally expensive.
Approximate methods such as tau-leaping or diffusion/SDE approximations may then be useful, but they trade some exact CTMC event-level detail for computational efficiency.
What Gillespie simulation gives us
One run gives one possible biological history. Repeated runs give a sample from the model's distribution of possible histories.
| From repeated simulations | Example biological question |
|---|---|
| extinction frequency | What is the probability the infection dies out? |
| outbreak frequency | What is the probability it reaches a specified outbreak level? |
| distribution of peak size | How uncertain is the epidemic peak? |
| distribution of peak time | When might the peak occur? |
| capacity exceedance frequency | What is the probability demand exceeds a hospital threshold? |
The complete picture
What comes next?
To estimate probabilities and summary statistics reliably, we need many independent stochastic runs rather than one. The next section introduces Monte Carlo simulation and shows how repeated Gillespie trajectories can be turned into quantitative estimates of biological uncertainty.