Why Test Delayed and Stochastic Code?
A simulation can produce a smooth trajectory or a plausible collection of sample paths and still be wrong. A delayed solver may query the history at the wrong time; a stochastic simulator may consume random numbers inconsistently; an event queue may process a completion late or twice. These errors can survive visual inspection, particularly when the code is used inside a larger QSP model.
A useful test starts from a specific model property and turns it into a checkable prediction. Some predictions are exact, such as the state change caused by one scheduled event. Others are approximate, such as agreement with a known solution up to a stated numerical tolerance. Tests should be small enough that a failure points to a particular part of the implementation.
For a delayed or stochastic model, a test specification should identify the model assumptions, the input state and history, the relevant event or random inputs, the expected output, and the tolerance if the output is numerical. A tolerance is not a substitute for an expected result: it should reflect the solver's accuracy and the scale of the quantity being checked.
Test Histories and Delayed Queries
A delay equation can require a state value at a time before the current integration time. Testing this operation separately from a full simulation helps isolate incorrect history initialization, boundary handling, and interpolation. Begin with an equation and a history for which the solution can be derived directly.
Worked Example: Checking a Delayed-History Query
Consider the scalar delay equation \(dx/dt=-k x(t-\tau)\), with \(x(t)=10\) nM for \(t\leq 0\), \(k=0.2\) day\(^{-1}\), and \(\tau=2\) days. During the first delay interval, \(0\leq t\leq 2\), the queried value is always the prescribed history value \(x(t-\tau)=10\) nM. Therefore,
At \(t=1\) day, the expected state is \(x(1)=10-2(1)=8\) nM; at \(t=2\) days it is \(x(2)=10-2(2)=6\) nM. A test can check both the history query, which should return 10 nM at the required delayed times, and the integrated state at these checkpoints. For a numerical solver, compare with the analytic values using an absolute tolerance appropriate to the solver configuration, for example \(10^{-4}\) nM in a deliberately tight-tolerance test.
The history check and the trajectory check diagnose different failures. If the query at a negative time is wrong, the history initialization or time indexing is suspect. If the queried history is correct but \(x(1)\) is not near 8 nM, inspect the integration and units. A test at the boundary \(t=\tau\) is particularly useful: the required delayed value is \(x(0)\), so an off-by-one or strict-versus-inclusive boundary error may become visible there.
For histories represented by stored points, also test a query between points, where interpolation is used, and a query exactly at a stored point. Test times just to either side of the boundary between prescribed history and simulated trajectory. The expected interpolation rule and boundary convention must be explicit; different valid implementations can use different internal representations, but they must return the values required by the stated model.
Check Scheduled Events and Random Inputs
In a delayed stochastic model, a reaction may create a future event that changes the state when its scheduled time arrives. A focused test should separate initiation from completion and verify the state and event queue at both stages. The test should reflect the model's specified delay semantics: for example, whether reactants are consumed at initiation and products appear at completion.
Worked Example: Verifying One Scheduled Completion
Suppose a delayed reaction consumes one \(A\) molecule at initiation and produces one \(B\) molecule after a fixed delay \(\tau=2\) min. Start with \(A=5\), \(B=0\), and no pending events. A test fixture forces an initiation at \(t=1.2\) min. Immediately after initiation, the expected counts are \(A=4\), \(B=0\), with one pending completion scheduled for \(1.2+2=3.2\) min.
Advance the event-processing component to \(t=3.199\) min. The completion must not yet have occurred: \(A=4\), \(B=0\), and one event remains pending. Advance to \(t=3.2\) min under the implementation's documented boundary convention. The event should complete exactly once, giving \(A=4\), \(B=1\), with no pending completion. Checking the time, state change, and pending-event count together detects errors that a final-state-only test could miss, such as an event being processed twice.
If a scheduled completion and a newly sampled reaction occur at the same time, the implementation needs a specified tie-handling convention. There is no need to impose one universal ordering in a test; instead, test that the code consistently applies its documented convention and does not lose or duplicate either event.
Stochastic unit tests should control random inputs whenever possible. A fixed random seed can support a regression test for a particular implementation, but a seed alone may not provide a stable test across software versions or changes in the order of random-number draws. Injecting the uniforms or using a test random-number generator with a known sequence makes the expected path explicit.
Worked Example: Testing a Stochastic Reaction Choice
Consider two reactions: \(A\to B\), with propensity \(a_1=0.5A\) min\(^{-1}\), and \(B\to C\), with propensity \(a_2=1.5B\) min\(^{-1}\). Start from \((A,B,C)=(4,2,0)\), so \(a_1=2\) min\(^{-1}\), \(a_2=3\) min\(^{-1}\), and \(a_0=a_1+a_2=5\) min\(^{-1}\). For a direct-method step, supply \(U_{\text{time}}=e^{-0.5}\) and \(U_{\text{select}}=0.7\), both in \((0,1)\).
The waiting time calculated from the supplied time uniform is
The reaction-selection threshold is \(U_{\text{select}}a_0=0.7(5)=3.5\) min\(^{-1}\). Because \(3.5>a_1=2\) and \(3.5\leq a_1+a_2=5\), the selected reaction is \(B\to C\). The expected next state is \((4,1,1)\) at \(0.1\) min after the starting time. This single test checks the waiting-time transform, the propensity interval used for selection, and the state-update rule. It should also verify that the implementation rejects or flags invalid uniform inputs such as 0, 1, or a value outside \([0,1)\), according to its documented random-number interface.
Test Limiting Cases and Conservation
Limiting cases offer checks that are independent of many implementation details. Setting a delay to zero should recover the corresponding no-delay model when the model's semantics make that limit well-defined. Setting a reaction propensity to zero should prevent that reaction from firing. Such tests can expose incorrect branches, stale events, or units errors without requiring a large simulation.
Worked Example: Testing the Zero-Delay Limit
Use the equation and history from the earlier example, but set \(\tau=0\). The model reduces to \(dx/dt=-kx\), with \(x(0)=10\) nM and \(k=0.2\) day\(^{-1}\). Its solution at \(t=1\) day is
A limiting-case test compares the delay solver with \(\tau=0\) against this value, or against a separately tested no-delay solver using the same numerical tolerances. Agreement supports the zero-delay branch, but does not prove that positive-delay history handling is correct; that requires tests such as the first-interval calculation above. For stochastic code, an analogous check can compare a zero-delay reaction implementation with its specified immediate-transition version, while ensuring that no zero-time event loop occurs.
Conservation tests turn stoichiometry into a property that can be checked after each transition. Suppose the delayed reaction in the scheduled-event example consumes \(A\) at initiation and creates \(B\) at completion. Let \(P\) be the number of pending completions. Then \(A+B+P\) is conserved: initiation changes the vector by \((-1,0,+1)\), and completion changes it by \((0,+1,-1)\). Both changes sum to zero.
Worked Example: Testing an Extended Conservation Invariant
Start with \((A,B,P)=(3,0,0)\), so \(A+B+P=3\). After one delayed initiation, the state is \((2,0,1)\), and \(2+0+1=3\). When that event completes, the state becomes \((2,1,0)\), and \(2+1+0=3\). A property-based test can check \(A+B+P=3\) after every initiation and completion across many generated event sequences. Counts should also remain nonnegative.
The pending-event count belongs in the invariant for this particular delayed-reaction convention. Checking only \(A+B\) would incorrectly report a loss at initiation, even though the molecule is represented by a pending completion. Conversely, a conservation check alone cannot catch every defect: exchanging the identities of two products might preserve the total. Pair invariant checks with event-specific assertions.
Build a Test Strategy and Diagnose Failures
A practical test suite combines different kinds of evidence rather than relying on one large simulation. The following sequence moves from isolated mechanics toward integrated behavior.
Check history values, propensity calculations, event-time arithmetic, and state updates with inputs whose outputs are known.
Include times at and around delay boundaries, zero propensity, zero delay where defined, and empty or full event queues.
Check nonnegative counts, conservation laws, and consistency between queued events and state bookkeeping after each update.
Run small end-to-end cases, then use repeated simulations and statistical checks when a property concerns distributions rather than one path.
Exact path tests and statistical tests serve different purposes. With controlled random inputs, an exact path test checks whether a specified sequence of draws produces the expected sequence of transitions. It does not establish that the overall output distribution is correct. A statistical test, by contrast, compares summaries over many replicates with a known or independently calculated distribution; its result has sampling variability and should not be treated as a deterministic pass/fail oracle without a justified tolerance.
| Approach | Best use | Main limitation |
|---|---|---|
| Analytic or hand-calculated oracle | Small cases with exact or known solutions | May not exist for a complex model |
| Controlled random inputs | Deterministic tests of a stochastic code path | Does not by itself test distributional correctness |
| Invariants and properties | Checks that should hold over many states or events | Some incorrect implementations can still preserve them |
| Statistical validation | Comparisons of simulated distributions or moments | Needs enough replicates and a tolerance that accounts for Monte Carlo variation |
When a test fails, first identify whether the discrepancy is in a deterministic calculation, a time boundary, a random draw, or a numerical approximation. Log the queried delayed time, current state, propensity values, selected event, event queue, and random inputs for a minimal failing case. For numerical comparisons, report the absolute error and the tolerance used. A tolerance that is too loose can conceal an indexing defect; one that is unrealistically tight can make a correct solver fail because of ordinary floating-point error.
Tests also have limits. A finite set of cases cannot prove that all event schedules are correct, and a conserved total cannot certify every state update. Statistical tests can miss small biases or flag correct code by chance. Keep tests small and diagnostic, document their assumptions, and use more than one independent check for important model behavior.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions.
- For a delayed history test, why should a query at the delay boundary be checked separately from an interior query?
- A reaction consumes one molecule at initiation and creates one product at completion. What quantity should be included in a conservation check while the event is pending?
- In the stochastic example, what reaction is selected if \(U_{\text{select}}=0.3\), with the same propensities?
- What does a test with injected random inputs establish, and what does it not establish?
Solutions. (1) At the boundary the queried time may switch between prescribed history and simulated values, so an indexing or boundary-convention error can appear there even when interior interpolation is correct. (2) Include the number of pending completions, \(P\), so the invariant is \(A+B+P\) under the stated reaction semantics. (3) The threshold is \(0.3(5)=1.5\) min\(^{-1}\), which lies in the first channel's interval \((0,2]\), so \(A\to B\) is selected. (4) It checks the code's response to a specified sequence of random inputs; by itself, it does not establish that the implementation produces the correct probability distribution over all paths.