Diagnosing Stochastic Simulation Bias and Convergence

Why Separate Diagnostics Matter

A stochastic simulation can produce plausible-looking trajectories and still be wrong in more than one way. Event times may be generated from the wrong rate, a reaction may violate a conserved quantity, a numerical approximation may distort the distribution, or a small number of replicates may leave the estimated mean imprecise. These problems require different checks: increasing the replicate count does not repair a biased algorithm, and reducing a step size does not correct an incorrectly implemented event rate.

The diagnostics here apply to stochastic cellular models described by reaction channels, state changes, and propensities, including the continuous-time Markov chain (CTMC) formulation used in earlier material on cellular CTMCs and Gillespie’s direct method. They are checks against a specified model, not evidence that the model represents biology adequately.

For a CTMC in state \(x\), let reaction channel \(j\) have propensity \(a_j(x)\), and let \(a_0(x)=\sum_j a_j(x)\) be the total propensity. In a state where these propensities remain constant until the next event, the exact waiting time \(\Delta t\) has an exponential distribution with rate \(a_0(x)\). Conditional on an event occurring, channel \(j\) is selected with probability \(a_j(x)/a_0(x)\). These statements provide direct checks of event timing and channel selection.

Check Event Timing and Event Selection

A useful event-time diagnostic transforms each observed waiting time into a dimensionless residual:

For an exact direct-method simulation with fixed propensities between events, the residuals should follow an exponential distribution with rate 1, written \(R_i\sim \mathrm{Exp}(1)\). Its mean and variance are both 1. If propensities vary continuously with time, the corresponding residual is the integrated total hazard over the waiting interval, \(\int_{t_i}^{t_{i+1}} a_0(X(t),t)\,dt\). Residuals that are consistently too large or too small can indicate that waiting times are being generated at the wrong rate. Also inspect their spread and distribution: matching the mean alone is not enough to establish a good fit.

Event selection can be checked separately. At a fixed state \(x\), repeated events should select channel \(j\) with probability \(a_j(x)/a_0(x)\). Compare observed channel proportions with these probabilities in controlled tests or in groups of events where the state and rate structure make the comparison meaningful. In an unconstrained trajectory, state changes alter the probabilities, so pooling all events and comparing with a single probability can be misleading.

Check Conservation Laws Against the State Changes

Some reaction networks conserve a total amount, such as the combined number of free and bound molecules. Write the state change for channel \(j\) as a column of the stoichiometric matrix \(S\). A vector \(c\) defines a conserved quantity when \(c^\mathsf{T}S=0\): every permitted reaction leaves \(c^\mathsf{T}x\) unchanged. Thus, if no external input, removal, or other explicitly modeled process changes that quantity, every simulated state should satisfy \(c^\mathsf{T}x(t)=c^\mathsf{T}x(0)\).

For a delayed reaction, the quantity being checked may need to include material in a pending event or an intermediate state. A molecule that has left one compartment but has not yet completed its delayed transition should not be counted as lost if the model treats it as still part of the conserved total. Define the accounting boundary from the model’s reaction scheme before interpreting a discrepancy.

Diagnose Step-Size Bias in Tau-Leaping

The exact SSA advances from one reaction event to the next; it does not approximate the event process by taking a fixed time step. Tau-leaping instead advances by intervals of length \(\tau\) and approximates reaction counts over each interval. Its results can therefore depend systematically on \(\tau\). A practical check is to simulate the same model and outputs with successively smaller step sizes, using enough replicates that Monte Carlo uncertainty does not obscure the differences. Where possible, compare with an exact SSA result or an analytical benchmark.

Consider a pure-birth process in which each cell produces one additional cell at rate \(k\) per cell, so the propensity is \(a(x)=kx\). Over a step of length \(\tau\), a Poisson tau-leap draws a number of births with mean \(kx\tau\). Consequently, its conditional mean at the end of the step is \(x+kx\tau=x(1+k\tau)\). The exact process instead has mean \(xe^{k\tau}\) after that interval. The local difference follows from the exponential expansion:

For positive \(k\) and small \(\tau\), the exact conditional mean exceeds the one-step tau-leap mean. Across multiple steps the discrepancy accumulates. This comparison is a justification for checking step-size sensitivity, not a universal direction-of-bias rule for every reaction network or every output.

Separate Replicate Convergence From Algorithmic Accuracy

For \(n\) independent replicate outputs \(Y_1,\ldots,Y_n\), the estimated mean is \(\bar{x}\), the sample standard deviation is \(s\), and the Monte Carlo standard error is \(SE=s/\sqrt{n}\). Under conditions that make the sampling distribution of the mean approximately normal, an approximate 95% Monte Carlo interval for the expected output is \(\bar{x}\pm1.96SE\). It describes uncertainty in the estimated expectation due to finite replication, conditional on the simulator and model being used. It is not an interval for an individual trajectory, and it does not account for algorithmic bias or model uncertainty.

A useful convergence check is to track the estimate and its standard error as independent replicates accumulate, and to inspect important outputs separately. A stable estimate with a small standard error indicates good Monte Carlo precision for that summary, not that the simulation algorithm is unbiased. For probabilities of rare events, the number of observed events and the uncertainty in the estimated probability matter; a narrow-looking interval for a common outcome does not establish precision for a rare one.

A Practical Diagnostic Sequence and Its Limits

These checks complement rather than replace one another. Comparing an approximate method with SSA can reveal disagreement but cannot, by itself, identify which result is wrong if the model or implementation is flawed. Agreement across step sizes can also be misleading if all runs share the same event-handling error. The exact SSA is exact for the specified CTMC under its assumptions; it is not automatically exact for a different biological system. Likewise, conservation checks apply only to quantities the model truly conserves. A quantity can change legitimately when the model includes sources, sinks, or material outside the chosen accounting boundary.

Common Pitfall

Wrong reasoning: “The mean changed very little when I doubled the number of replicates, so the algorithm has converged and is unbiased.” Replicate stability addresses Monte Carlo variability in an estimate; it does not test whether event rates, stoichiometric updates, or tau-leap approximations are correct.

Correction: Check event timing and state updates against the model, assess conserved quantities where applicable, and test step-size sensitivity for approximate algorithms. Then use replicate standard errors to quantify the remaining Monte Carlo uncertainty. Report these as distinct diagnostics.