Gillespie's Direct Method from Exponential Waiting Times

From Propensities to a Stochastic Event

In a cellular continuous-time Markov chain (CTMC), a state records the system’s current configuration, and each possible transition has a propensity that gives its instantaneous event rate in that state. The earlier material on formulating cellular CTMCs establishes how to specify those states, transitions, and propensities. Gillespie’s direct method turns that specification into a simulation: it determines both when the next event occurs and which event occurs.

The method is exact for a CTMC whose propensities remain constant between state changes. Its central idea is a race: at the current state, each enabled event channel has an exponentially distributed waiting time. The first clock to ring determines the next transition. Since the state changes at that event, all propensities must then be evaluated again at the new state.

Definition: Let \(X(t)=x\) be the current state, with \(M\) possible event channels. For channel \(j\), let \(a_j(x)\geq 0\) be its propensity and let \(\nu_j\) be its state-change vector. The total propensity is \(a_0(x)=\sum_{j=1}^{M}a_j(x)\). A propensity has units of inverse time and is the instantaneous rate of that channel at state \(x\).

A zero propensity means that the corresponding channel cannot fire in the current state. If every channel has zero propensity, then \(a_0(x)=0\): no event can occur, so the state is absorbing under the specified model. Otherwise, the direct method draws a waiting time and selects one channel from the current propensities.

Deriving the Waiting Time and Event Type

Fix the current state \(x\), and suppose each channel \(j\) has an independent exponential clock \(T_j\) with rate \(a_j(x)\). A channel with zero propensity has \(T_j=\infty\) for this race. For a positive rate, its survival probability is \(P(T_j>t)=e^{-a_j(x)t}\). The next event time is the minimum of the clocks, \(T=\min_j T_j\).

The event that no channel fires by time \(t\) is the event that every clock exceeds \(t\). Independence therefore gives the survival probability of the minimum:

P(T>t)=\prod_{j=1}^{M} e^{-a_j(x)t} =e^{-a_0(x)t}, \qquad a_0(x)=\sum_{j=1}^{M}a_j(x).

Thus the time to the next event is exponential with rate \(a_0(x)\). Its density is \(a_0(x)e^{-a_0(x)t}\), and its mean is \(1/a_0(x)\). To find the event type as well, the joint density that channel \(j\) fires at time \(t\) while all other clocks have not yet fired is \(a_j(x)e^{-a_0(x)t}\). Integrating over all possible event times gives the probability that channel \(j\) is first:

P(\text{next channel}=j\mid X(t)=x) =\int_0^\infty a_j(x)e^{-a_0(x)t}\,dt =\frac{a_j(x)}{a_0(x)}.

In fact, the joint density factors into the waiting-time density and the event-type probability:

a_j(x)e^{-a_0(x)t} = \bigl(a_0(x)e^{-a_0(x)t}\bigr) \frac{a_j(x)}{a_0(x)}.

This factorization shows that, conditional on the current state, the waiting time and identity of the next channel are independent. The wait depends on the sum of the propensities; the relative propensities determine which channel wins.

The Direct-Method Calculation

Two independent draws from a uniform distribution on \((0,1)\) implement these results. If \(U_1\) determines the waiting time and \(U_2\) selects the channel, then:

\tau=-\frac{\ln U_1}{a_0(x)}.

For the channel choice, form cumulative sums \(A_j(x)=\sum_{k=1}^{j}a_k(x)\), with \(A_0(x)=0\). Select the smallest index \(j\) for which \(A_j(x)>U_2a_0(x)\). This selects channel \(j\) with probability \(a_j(x)/a_0(x)\). The strict inequality is a convenient convention; exact equality has probability zero for a continuous uniform draw.

Evaluate propensities.
At the current state \(x\), calculate every \(a_j(x)\) and their sum \(a_0(x)\). Check that all propensities are nonnegative and finite.

Check for an absorbing state.
If \(a_0(x)=0\), no event can occur under the model. Otherwise, continue.

Draw the waiting time.
Generate \(U_1\) uniformly on \((0,1)\) and set \(\tau=-\ln(U_1)/a_0(x)\).

Select a channel.
Generate an independent \(U_2\) uniformly on \((0,1)\). Choose the first \(j\) whose cumulative propensity exceeds \(U_2a_0(x)\).

Update and repeat.
Advance time by \(\tau\), update the state to \(x+\nu_j\), and recalculate all propensities at that new state before drawing again.

Worked Calculations

Worked Example: Competing Cellular Events

Hypothetical data. At a particular state, a cell has 20 free receptors, 12 free ligand molecules, and 3 receptor–ligand complexes. Consider three channels: binding, unbinding, and complex internalization. Let the binding rate constant be \(0.002\) cell\(^{-1}\) min\(^{-1}\), and let the unbinding and internalization rates be \(0.15\) min\(^{-1}\) and \(0.05\) min\(^{-1}\), respectively. Suppose the propensities are \(a_{\text{bind}}=k_{\text{on}}RL\), \(a_{\text{unbind}}=k_{\text{off}}C\), and \(a_{\text{int}}=k_{\text{int}}C\), where \(R,L,C\) are the current counts.

Calculate the channel propensities.

a_{\text{bind}}=(0.002)(20)(12)=0.48\ \text{min}^{-1}, \quad a_{\text{unbind}}=(0.15)(3)=0.45\ \text{min}^{-1}, \quad a_{\text{int}}=(0.05)(3)=0.15\ \text{min}^{-1}.

The total propensity is \(a_0=0.48+0.45+0.15=1.08\) min\(^{-1}\). Therefore the time to the next event is exponential with rate \(1.08\) min\(^{-1}\), and its mean is \(1/1.08=0.9259\) min. The probabilities for the next event to be binding, unbinding, and internalization are \(0.48/1.08=0.4444\), \(0.45/1.08=0.4167\), and \(0.15/1.08=0.1389\), respectively. They sum to 1, as required.

Worked Example: Turning Uniform Draws into an Event

Use the hypothetical state and propensities from the preceding example. Suppose the independent random draws are \(U_1=0.25\) and \(U_2=0.70\).

Waiting time. With \(a_0=1.08\) min\(^{-1}\),

\tau=-\frac{\ln(0.25)}{1.08} =\frac{1.3863}{1.08} =1.284\ \text{min}.

Channel. The channel order is binding, unbinding, then internalization. The cumulative propensities are \(0.48\), \(0.48+0.45=0.93\), and \(1.08\) min\(^{-1}\). The selection threshold is \(U_2a_0=(0.70)(1.08)=0.756\) min\(^{-1}\). Since \(0.48<0.756\leq0.93\), unbinding is selected. If the unbinding transition converts one complex into one free receptor and one free ligand, the updated counts are \(R=21\), \(L=13\), and \(C=2\). The simulation must calculate new propensities from those updated counts.

Worked Example: Probability of an Event Before a Fixed Time

Return to the same hypothetical state, where \(a_0=1.08\) min\(^{-1}\). Find the probability that any event occurs within \(0.5\) min, and the probability that the first event is unbinding and occurs within that time.

The probability that no event occurs by \(t=0.5\) min is the exponential survival probability \(e^{-a_0t}\). Thus:

P(T\leq0.5)=1-e^{-(1.08)(0.5)} =1-e^{-0.54} =0.4173.

For the joint event that unbinding is first and occurs by \(0.5\) min, integrate its joint density from zero to \(0.5\). Equivalently, multiply the probability that unbinding is the first channel by the probability that some event occurs by that time:

P(\text{unbinding first},\,T\leq0.5) =\frac{0.45}{1.08}\left(1-e^{-(1.08)(0.5)}\right) =(0.4167)(0.4173) =0.1739.

This is a joint probability, not the probability of unbinding conditional on an event having occurred by \(0.5\) min. Under these constant propensities, that conditional probability is \(0.45/1.08=0.4167\).

Assumptions, Diagnostics, and Alternatives

The direct method’s exactness depends on the process being represented correctly as a CTMC. In particular, the current state must contain the information needed to determine future event rates; while the state remains unchanged, each propensity must be constant; and the specified propensities must be valid nonnegative rates. The competing-exponential derivation then applies at each state. After a jump, the new state determines a new race. The exponential memoryless property is why this stepwise procedure is consistent with the CTMC.

Several implementation checks help detect errors before interpreting a simulation:

Check rates and units. Every propensity contributing to \(a_0\) must have units of inverse time. A negative or non-finite propensity signals a problem in the model or its implementation.
Check state feasibility. A channel that requires a molecule absent from the state should have zero propensity. A selected transition should not create an impossible state if the transition definition is valid.
Check the cumulative selection. The cumulative sum must end at \(a_0\), and the selected channel must satisfy \(A_{j-1}\leq U_2a_0<A_j\). Incorrect channel ordering is harmless only if the same ordering is used consistently in the sum and selection.
Recompute after every event. Using propensities from the previous state after a transition gives the wrong event-time and event-type distributions whenever that transition changes a rate.
Handle the zero-total case. Dividing by \(a_0=0\) is undefined; it is not a very long waiting time. The specified process has no further events from that state.

For a small number of channels, the first-reaction method is a direct alternative: draw a separate exponential waiting time for each channel and choose the smallest. Under the same constant-propensity assumptions, this gives the same race distribution, but it requires one draw per channel rather than a draw for the total waiting time and another for channel selection. A deterministic differential-equation model is a different kind of alternative: it describes continuous population-level change rather than generating individual random event times. It does not reproduce the event-by-event stochastic trajectory of the CTMC.

The direct method in this form is not appropriate if propensities vary continuously with time while the state is fixed, or if the model’s waiting times are not exponential and the current state does not capture the relevant history. In such cases, the constant-rate race derived here does not describe the event-time distribution. The model or simulation method must account for the actual rate or waiting-time structure rather than treating a changing propensity as constant.

Common Pitfall

Common mistake: “The binding propensity is the largest, so binding must occur next.” A larger propensity makes binding more likely, but does not guarantee it. Here binding is selected with probability \(0.4444\), while the other channels together have probability \(0.5556\). Also, the waiting time is governed by the total propensity \(1.08\) min\(^{-1}\), not by the propensity of whichever channel is selected.

A related error is to advance time by \(1/a_0\) on every step. That quantity is the mean waiting time, not a fixed interval. The method draws a new exponential waiting time at each state, so simulated intervals vary even when the state-dependent total propensity happens to be the same.

Check Your Understanding

Try each question, then check the short solutions below.

At a state, two channels have propensities \(0.3\) min\(^{-1}\) and \(0.7\) min\(^{-1}\). What are the total propensity, mean waiting time, and probability that channel 2 fires next?
For the propensities in question 1, use \(U_1=0.5\). Calculate the waiting time to three significant figures.
Explain why a channel with zero propensity has zero probability of being selected when \(a_0>0\).
What should the simulator do if every propensity at the current state is zero?

Solutions. (1) \(a_0=0.3+0.7=1.0\) min\(^{-1}\); the mean wait is \(1/a_0=1.00\) min; channel 2 has probability \(0.7/1.0=0.700\). (2) \(\tau=-\ln(0.5)/1.0=0.693\) min. (3) Its selection probability is \(a_j/a_0=0/a_0=0\). (4) It should identify an absorbing state under the specified model and stop generating events from that state; it should not divide by zero or invent a waiting time.