Validating Gillespie Simulations with Rate and Time Checks

What a Gillespie Validation Check Can Tell You

A Gillespie simulation produces both a waiting time and an event type at each step. These outputs provide two distinct checks: at a fixed state, holding times should follow an exponential distribution with rate equal to the sum of the propensities, while the event type should follow a categorical distribution proportional to those propensities. A third check is to compare summaries of a small simulated model with calculations from its generator.

These checks are diagnostic, not demands for exact agreement. A finite sample fluctuates, so a discrepancy is evidence to investigate rather than automatic proof of a coding error. In particular, compare each simulated event with the predictions at the state before that event, and separate sampling variability from errors in rates, state updates, or bookkeeping.

Consider a continuous-time process \((X_t)_{t\ge0}\) on a state space \(S\), defined on a probability space \((\Omega,\mathcal{F},\mathbb{P})\). Let \((\mathcal{F}_t)_{t\ge0}\) denote the history of the process up to time \(t\). At state \(i\), suppose there are \(r\) possible reaction channels with propensities \(a_1(i),\ldots,a_r(i)\), all finite and nonnegative. Write \(a_0(i)=\sum_{\ell=1}^r a_\ell(i)\). The direct Gillespie step is defined when \(a_0(i)>0\): it draws a holding time \(\tau\) and an event label \(J\in\{1,\ldots,r\}\). If \(a_0(i)=0\), no event can occur from \(i\), so the process remains there.

Proof. Fix the state \(i\), and abbreviate \(a_\ell=a_\ell(i)\) and \(a_0=a_0(i)>0\). One construction of the direct step assigns independent exponential clocks with rates \(a_1,\ldots,a_r\) to the channels with positive rates; a zero-rate channel never rings. For \(t\ge0\), the event that the minimum clock exceeds \(t\) is the event that every positive-rate clock exceeds \(t\). Independence and the exponential survival formula give \(\mathbb{P}(\tau>t)=\prod_{\ell:a_\ell>0}e^{-a_\ell t}=e^{-a_0t}\). For a channel \(\ell\) with \(a_\ell>0\), the joint density that its clock rings at time \(u\) and all other clocks exceed \(u\) is \(a_\ell e^{-a_0u}\) for \(u>0\). Integrating gives \(\mathbb{P}(J=\ell)=\int_0^\infty a_\ell e^{-a_0u}\,du=a_\ell/a_0\); for \(a_\ell=0\), the probability is zero, also equal to \(a_\ell/a_0\). The joint density of \(\tau=u,J=\ell\) is \(a_\ell e^{-a_0u}=(a_0e^{-a_0u})(a_\ell/a_0)\). This factors into the density of \(\tau\) and the probability of \(J=\ell\), proving conditional independence. \(\square\)

The theorem turns implementation details into concrete tests. At each recorded pre-event state, save the propensities, their sum, the holding time, and the event label. Check times against that state’s total rate, and labels against that state’s propensity proportions. Pooling observations with different rates can conceal an error or create an apparent discrepancy.

Checking Holding Times

For a fixed rate \(a_0>0\), an exponential holding time has mean \(1/a_0\), variance \(1/a_0^2\), and survival probability \(e^{-a_0t}\). Thus, across \(n\) holding times generated at the same state (or at states with the same total rate), a natural rate estimate is \(n/\sum_{k=1}^n\tau_k\). The exponential model also implies that \(1-e^{-a_0\tau_k}\) has a uniform distribution on \([0,1]\): for \(0\le u\le1\), its cumulative probability is \(\mathbb{P}(1-e^{-a_0\tau_k}\le u)=u\). This transform provides a diagnostic beyond the sample mean.

These checks require care about which intervals are included. Include the holding time ending in an event, but do not treat a simulation’s final, unfinished interval as a completed exponential observation without accounting for its censoring. A state with total rate zero contributes no event-terminated holding time. When rates change from state to state, use \(a_0(i)\) for each interval’s pre-event state rather than one pooled rate.

Checking Event-Type Frequencies

At state \(i\), the expected relative frequency of channel \(\ell\), conditional on an event occurring there, is \(a_\ell(i)/a_0(i)\). For \(n_i\) recorded departures from that state, the expected count for that channel is \(n_i a_\ell(i)/a_0(i)\). The events need not be independent across an entire changing-state trajectory; the useful comparison is conditional on the pre-event state. Repeated visits to the same state provide repeated categorical draws with the same probabilities, provided the propensities are unchanged there.

Event labels can also reveal bookkeeping bugs that waiting times will not detect. For example, the holding-time law depends on the sum of propensities, so two different errors in channel-specific rates could leave that sum unchanged. The event-type check can expose the incorrect allocation. Conversely, correct event proportions do not validate the total rate used to draw times.

Generator-Based Checks for a Small Model

For a finite-state continuous-time Markov chain, the generator \(Q=(q_{ij})\) records the rates of transitions between distinct states; its diagonal entries make every row sum to zero. In a reaction model, rates for channels that move from \(i\) to \(j\) are added to obtain \(q_{ij}\). A generator can give benchmarks beyond individual Gillespie steps, including state probabilities and expected occupation times.

Proof. For a finite-state chain, the semigroup property and the definition of the generator give \(P(h)=I+hQ+o(h)\) as \(h\downarrow0\), entry by entry. Thus the row distribution \(\mu(t)=\mu_0P(t)\) satisfies \(\mu'(t)=\mu(t)Q\) and \(\mu(0)=\mu_0\). The matrix exponential \(e^{tQ}\), defined by its convergent power series, satisfies this differential equation and initial condition. Uniqueness for a finite-dimensional linear differential equation gives \(\mu(t)=\mu_0e^{tQ}\).

For occupation time, the indicator is nonnegative, so integration and expectation can be interchanged; its expectation at time \(t\) is \(\mathbb{P}(X_t=j)=\mu_j(t)\). This gives the stated integral. For transition counts, conditional on the history at time \(t\) and \(X_t=i\), the probability of a transition \(i\to j\) in an interval of length \(h\) is \(q_{ij}h+o(h)\); the probability of two or more jumps is \(o(h)\) for a finite-state chain with bounded rates. Summing these conditional expected increments over a partition of \([0,T]\) and taking the mesh to zero gives the integral of \(q_{ij}\mathbb{P}(X_t=i)\). Hence the expected count is \(q_{ij}\int_0^T\mu_i(t)\,dt\). \(\square\)

Common Mistakes and a Practical Checklist

Common mistake: pooling state-specific event probabilities. Suppose there are two channels. At state \(A\), their rates are \(2\) and \(1\), so the first-channel probability is \(2/3\). At state \(B\), their rates are \(1\) and \(3\), so that probability is \(1/4\). If a validation set contains 60 departures from each state, the conditional expected first-channel count is \(60(2/3)+60(1/4)=40+15=55\), or \(55/120=11/24\) of the events. An incorrect calculation pools the channel rates to get \(60(2+1)=180\) and \(60(1+3)=240\), then reports \(180/(180+240)=3/7\). That calculation weights states by their total rates instead of comparing the observed event labels conditional on the recorded states. Keep the states separate and use the propensity proportions at each one.

Before trusting a Gillespie output, check the following:

• At each state, recompute all propensities and verify that the total rate is their sum.
• Compare holding times with the exponential rate for their own pre-event state.
• Compare event counts with the categorical probabilities for their own pre-event state.
• For a small model, check that the generator rows sum to zero and compare simulated summaries with generator-based benchmarks.
• Use a suitable amount of simulation and interpret finite-sample differences as variability unless there is a persistent or substantial mismatch.

A good validation plan uses independent checks: waiting times, event types, and model-level summaries test different parts of the implementation. Agreement across them is stronger evidence than agreement with any single statistic; disagreement helps locate whether the issue lies in rates, event selection, state changes, or simulation-time accounting.