Why Gillespie Waiting Times and Event Choices Are Exact

Why a Single Waiting Time and One Event Choice Work

In a continuous-time reaction system, several events may be possible at the same state. Each has its own rate, or propensity, and the next change is whichever event happens first. Gillespie’s direct method replaces this collection of competing clocks with one waiting-time draw and one event-choice draw. The method is exact because those two draws reproduce the joint distribution of the first clock to ring and its event label.

Let \((\Omega,\mathcal{F},\mathbb{P})\) be a probability space and let \((X_t)_{t\ge 0}\) be a continuous-time process with state space \(S\), taken here to be finite or countable. At state \(x\), suppose there are finitely many event channels \(r\in\{1,\dots,m\}\). Channel \(r\) has propensity \(a_r(x)\ge 0\), and when it occurs it changes the state according to a specified rule \(\phi_r(x)\). Assume the propensities are finite. A channel with propensity zero cannot occur at that state. The total propensity is

For the moment, fix a state \(x\) and write \(a_r=a_r(x)\). Model each possible event by an independent exponential clock \(T_r\) with rate \(a_r\): for \(a_r>0\), \(\mathbb{P}(T_r>t)=e^{-a_rt}\) for \(t\ge 0\). A zero-rate clock is interpreted as \(T_r=\infty\). The next event time is \(T=\min_r T_r\), and its label is the index \(J\) of the first clock to ring.

Proof. Independence of the clocks gives, for \(t\ge 0\),

This is the survival function of an exponential random variable with rate \(a_0\). For any channel \(r\) with \(a_r>0\), the probability that its clock rings in a short interval from \(t\) to \(t+dt\), while all other clocks have not rung by time \(t\), has density \(a_re^{-a_0t}\). Thus the joint density of \(T\) and the event \(J=r\) is \(a_re^{-a_0t}\). Integrating over \(t\ge 0\) gives \(\mathbb{P}(J=r)=\int_0^\infty a_re^{-a_0t}\,dt=a_r/a_0\). Dividing the joint density by the density \(a_0e^{-a_0t}\) of \(T\) gives \(\mathbb{P}(J=r\mid T=t)=a_r/a_0\), which does not depend on \(t\); hence \(T\) and \(J\) are independent. A positive-rate exponential variable has a continuous distribution, so two positive-rate clocks tie with probability zero. A zero-rate clock never rings. Finally, if \(a_0=0\), every nonnegative \(a_r\) is zero, so all clocks equal infinity and no event occurs. \(\square\)

The theorem gives the direct-method step at state \(x\). If \(a_0(x)>0\), draw a uniform random variable \(U\) on \((0,1)\) and set \(\Delta t=-\log(U)/a_0(x)\). Separately choose channel \(r\) with probability \(a_r(x)/a_0(x)\), for example by partitioning the unit interval into consecutive intervals of those lengths. The waiting time and the event label must be generated independently. Then advance time by \(\Delta t\), apply \(\phi_r\), recompute propensities at the new state, and repeat. If \(a_0(x)=0\), the state has no possible event and the process remains there.

Why State-Dependent Propensities Preserve the Law

After an event, the state usually changes, so the propensities must be recalculated. The clock calculation applies conditionally at every state: given the current state \(x\), the next waiting time and event label have the law just derived using \(a_r(x)\). If fresh independent random draws are used after each event, this gives the same conditional evolution as restarting the competing clocks at the new state. In a reaction model with finite state space, this constructs the continuous-time Markov chain with the specified event rates. For countably many states, the same stepwise construction applies; to describe a process for all finite times, one also requires that infinitely many events do not accumulate in finite time.

Proof. Fix a state \(x\). There are finitely many states and channels, so the total rates are bounded above by some finite \(R\). In an interval of length \(h>0\), the probability of at least one event is \(1-e^{-a_0(x)h}=a_0(x)h+O(h^2)\), by the exponential expansion. The probability of two or more events is \(O(h^2)\): the probability of a first event is at most \(Rh\), and, conditional on its occurrence, the probability of another event within the remaining interval is at most \(Rh\); their product is at most \(R^2h^2\). Therefore, for \(y\ne x\), the probability of being at \(y\) at time \(h\), given state \(x\) at time zero, is the probability that a single event leading to \(y\) occurs, up to an error \(O(h^2)\). Summing the rates of those channels gives \(\mathbb{P}(X_h=y\mid X_0=x)=q_{xy}h+O(h^2)\). The probability of remaining at \(x\) is consequently \(1-\sum_{y\ne x}q_{xy}h+O(h^2)=1+q_{xx}h+O(h^2)\). These are the infinitesimal transition rates, and the definition of \(q_{xx}\) gives \(\sum_y q_{xy}=0\) for every row. \(\square\)

If a channel changes no state, it can still be simulated as a reaction event, but it is not a state change. In that case the total rate for the next reaction event includes that channel, while the off-diagonal generator rates include only channels that actually change the state. The proposition’s assumption excludes this distinction so that its total event rate is also the total state-change rate.

Exactness and a Common Mistake

Here, “exact” means that, conditional on the current state, the algorithm produces the same waiting-time and event-label distribution as the competing-clock model. It does not mean that a computer represents a continuous random number or a logarithm without numerical approximation. In a finite-state model, repeating the exact conditional step gives the intended continuous-time Markov chain. A simulation remains random: exactness concerns its probability law, not whether an individual simulated path matches a predetermined path.

Common mistake: choose an event first, then use that event’s rate for the waiting time. With rates \(1\) and \(3\), the event probabilities are \(1/4\) and \(3/4\), but the waiting time is exponential with total rate \(4\), regardless of which event wins. If one instead selects a channel first and then uses an exponential waiting time with that channel’s rate, the mean waiting time becomes \((1/4)(1)+(3/4)(1/3)=1/4+1/4=1/2\). The correct mean is \(1/4\), since the correct total rate is \(1+3=4\). Equivalently, the incorrect survival probability is the mixture \((1/4)e^{-t}+(3/4)e^{-3t}\), not the correct \(e^{-4t}\). The waiting time and event label are independent under the competing-clock law, but the waiting time is governed by the sum of the rates, not the selected channel’s rate.

The method also has a degenerate case that should not be forced through the formula. If all propensities at \(x\) are zero, then \(a_0(x)=0\), and dividing by \(a_0(x)\) is undefined. The correct interpretation is that no event can occur from \(x\); the state is absorbing under these event rules. If some but not all propensities are zero, the total is positive and the zero-rate channels receive event probability zero.

The practical check at each simulation step is therefore simple: add all current propensities; use that sum as the rate of one exponential waiting time; choose a channel with probability equal to its propensity divided by that same sum; then update the state and recompute. The clock proof explains why each part is required and why the resulting step is exact.