The Gillespie Direct Method for Continuous-Time Markov Chains

From Event Rates to a Simulation Step

A continuous-time Markov chain can change state at irregular times, and its possible moves may have different rates at different states. The Gillespie direct method simulates such a path one event at a time. At the current state, it combines the event rates into one total rate, draws a waiting time, chooses which event occurs, and updates the state. The key practical point is that rates must be recalculated after every state change.

Let \((\Omega,\mathcal{F},\mathbb{P})\) be a probability space carrying the random draws used by the simulation. Let \((X_t)_{t\ge0}\) be the simulated process with state space \(S\), where \(S\) is finite or countable, and let the time index \(t\) be continuous. Suppose there are \(R\) event channels. At state \(x\in S\), channel \(r\) has propensity \(a_r(x)\ge0\) and state-change vector \(\nu_r\); if that channel is selected, the new state is \(x+\nu_r\). A channel that is impossible at \(x\) has propensity zero. The total propensity is \(a_0(x)=\sum_{r=1}^R a_r(x)\), assumed finite at every state visited. The propensities are rates, not probabilities: they need not sum to one.

After the update, the method returns to the rate calculation using the new state. Thus the waiting time and event choice are based on the propensities at the state currently occupied, not on the propensities from the start of the simulation. For a simulation up to a fixed time \(T\), stop if the next event time would exceed \(T\): record the state as unchanged from the current time through \(T\), and do not apply that event.

Why the Two Random-Number Conversions Work

The waiting-time calculation is an inverse-transform step. The exponential distribution with rate \(a>0\) has survival probability \(e^{-as}\) at \(s\ge0\). The event-choice calculation is a categorical draw with weights proportional to the propensities. These facts justify the two random-number transformations themselves; connecting them to the waiting time and event identity of a continuous-time Markov chain is a separate question.

Proof. Since \(U\) is uniform on \((0,1)\), for any \(s\ge0\), \(\mathbb{P}(W>s)=\mathbb{P}(-\ln(U)/a>s)=\mathbb{P}(U<e^{-as})=e^{-as}\). The last equality follows because \(e^{-as}\in(0,1]\) and a uniform random variable on \((0,1)\) lies below a threshold \(c\in[0,1]\) with probability \(c\). This is the exponential survival function, including \(\mathbb{P}(W>0)=1\). \(\square\)

Proof. Define \(c_0=0\) and \(c_r=\sum_{j=1}^{r}a_j/a_0\) for \(r=1,\ldots,R\). Then \(0=c_0\le c_1\le\cdots\le c_R=1\). Event \(r\) is selected exactly when \(c_{r-1}\le U<c_r\), apart from endpoints, which have probability zero under a continuous uniform distribution. The interval has length \(c_r-c_{r-1}=a_r/a_0\), so that is the probability of selecting event \(r\). If \(a_r=0\), the interval has length zero. These intervals cover \((0,1)\) except for their shared endpoints, so one event is selected with probability one. \(\square\)

In implementation, the strict inequality in the selection rule gives a definite convention at a cumulative boundary: a value exactly on the boundary goes to the next positive-width interval. Since a continuous uniform draw equals any one of the finitely many boundaries with probability zero, this convention does not alter the probabilities. A random-number generator should supply values strictly between zero and one for the waiting-time formula; if an implementation can return zero, it should use a documented endpoint-safe conversion rather than evaluate \(\ln(0)\).

Worked Simulations

Implementing the Method Reliably

A direct implementation follows the same order at every step. Let \(t\) and \(x\) denote the current time and state. If the requested simulation horizon is \(T\), the implementation should first ensure \(t\le T\), then evaluate every propensity at \(x\). The following checklist makes the order explicit:

The state space and event updates must be compatible: a selected event should not produce an invalid state. Often this is ensured by assigning zero propensity to a reaction whenever its reactants or required population are unavailable. The algorithm also assumes that the total propensity can be computed at each visited state. In models where infinitely many events might occur before a finite time, additional care is needed: the event-by-event procedure could require infinitely many updates before reaching the horizon. A finite or suitably controlled simulation requires that such accumulation of event times does not occur.

Common mistake: treating rates as probabilities. In the infection example, a mistaken calculation might use \(\tfrac45\) and \(1\) directly as event probabilities. They sum to \(\tfrac95\), not \(1\), so they cannot be probabilities. The corrected probabilities are \((\tfrac45)/(\tfrac95)=\tfrac49\) for infection and \(1/(\tfrac95)=\tfrac59\) for recovery; their sum is \(\tfrac49+\tfrac59=1\). Another error is to divide the waiting-time numerator by the selected event's propensity. The waiting time uses the total rate \(a_0\), because one first uses the combined rate to draw the time and then uses the relative propensities to choose an event.

Finally, a simulated event time is not itself a fixed time step. The method advances the clock by the random waiting time generated at the current state. If an event changes the state, the rate landscape changes with it, so the next step begins with a fresh total propensity and fresh random draws. Keeping those two operations—advance by a waiting time, then recompute rates after the state update—separate is central to a correct implementation.