Exact Simulation of Delayed Reaction Events

Why Delays Require an Event Queue

In an ordinary stochastic simulation, each reaction changes the state as soon as it fires. A fixed-delay reaction has two relevant times: the time it is initiated and the later time when its completion effect occurs. Between those times, other reactions may fire and change which reactions are possible. An exact simulation must therefore represent both the current state and the reactions already scheduled to complete.

The key computational distinction is between a reaction initiation and a reaction completion. An initiation is selected stochastically using the current propensities. A completion occurs at a previously determined time; it is not a new random reaction firing. The simulator must compare the next possible stochastic initiation with the earliest scheduled completion and process whichever comes first.

This two-update description makes the reaction semantics explicit. For example, a reaction that consumes reactants when it begins and produces products after a delay has \(\nu_j^{I}\) equal to the reactant loss and \(\nu_j^{C}\) equal to the product gain. A different biological assumption may reserve reactants until completion; that requires different state updates and feasibility rules. The scheduling algorithm does not decide which interpretation is biologically appropriate.

Assume that, between events, the state is constant and each initiation propensity is therefore constant. Given state \(x\), the waiting time to the next initiation is exponentially distributed with rate \(a_0(x)=\sum_j a_j(x)\), and channel \(j\) is selected with probability \(a_j(x)/a_0(x)\). Propensities must be nonnegative and consistent with the model’s available molecules. The updates at initiation and completion must also leave the state valid. These assumptions are the basis for exact event scheduling.

The Event-Scheduling Algorithm

Maintain the current time \(t\), the state \(x\), and a priority queue of pending completions ordered by their scheduled times. At each step, calculate the next stochastic initiation time from the current propensities, then compare it with the earliest pending completion. A queued completion is processed before any later proposed initiation. Once a state-changing event occurs, the old proposed initiation is no longer valid: calculate propensities again and draw a new candidate.

A reaction with no delay can be represented as an immediate state change in its initiation step. The description above treats delayed channels as having positive fixed delays, avoiding the special case of a zero-delay event being inserted at the current time.

Why the Comparison Gives an Exact Simulation

Suppose the current propensities are constant up to the next scheduled completion, which occurs after \(d\) time units. Before that completion, the probability that no initiation occurs for \(u\) time units is \(e^{-a_0(x)u}\). Thus the candidate waiting time drawn from an exponential distribution is the correct waiting time for the next initiation if the state remains unchanged.

If the candidate time is later than the completion, the state changes at the completion and may have new propensities. The candidate drawn under the old state must not be used to select a later reaction. Instead, recompute the rates in the new state. The exponential waiting-time distribution is memoryless: conditional on no initiation before a fixed time \(d\), the remaining wait has the same exponential distribution with rate \(a_0(x)\). Therefore drawing afresh after processing the completion gives the correct post-completion waiting time.

At an initiation, the channel probabilities \(a_j(x)/a_0(x)\) are the standard competing-hazards probabilities. The channel selection is made using the pre-initiation state; after applying its initiation update, the propensity calculation is repeated. Together, this race among initiations and the exact ordering of scheduled completions reproduce the model’s event process, provided the stated assumptions hold.

New Initiations Do Not Cancel Earlier Completions

A scheduled completion remains in the queue when another reaction initiates, unless the model explicitly specifies a cancellation mechanism. In particular, if reactants are consumed at initiation, later changes in their availability do not undo a completion that was already scheduled. This is why the initiation update must encode what is consumed immediately and the completion update must encode what is produced later.

Multiple Completions, Diagnostics, and Alternatives

Two or more events can have the same completion time, especially when they share a fixed delay and were initiated simultaneously or at the same recorded time. If their completion updates are additive and the model defines them as simultaneous, apply their combined update as one batch, then recalculate propensities. For example, two pending completions each adding one \(Y\) change \(Y\) by two at their shared time. If updates are not additive, or their order could affect validity, the model must specify how coincident completions are resolved; an arbitrary software ordering should not silently become a biological assumption.

Under a continuous exponential waiting-time model, an independently drawn initiation has probability zero of landing at exactly a fixed completion time. Nonetheless, ties can occur from finite-precision arithmetic or when multiple completions share a scheduled time. A practical implementation should use a documented tie rule: batch compatible simultaneous completions, apply their updates consistently, and then recalculate propensities. It should not execute an initiation based on rates calculated before those completions.

Useful diagnostics follow directly from the event semantics. Check that completion times are no earlier than their initiations, that the queue is processed in time order, and that every initiation uses propensities from the current state. Track the queue in a small hand-worked case and verify that each scheduled completion occurs once. If species counts become negative, investigate the reaction updates and feasibility rules; do not repair the path by truncating counts to zero. Also check that a newly initiated event does not accidentally overwrite or erase an unrelated pending event.

This method differs from tau-leaping, which approximates groups of reaction firings over a time step, and from deterministic delay equations, which describe continuous-valued population dynamics rather than individual stochastic events. Transit-chain models can represent delay dispersion through extra stochastic or deterministic stages, rather than one fixed completion time per initiation. Those approaches may be useful for different scientific questions, but they do not replace exact scheduling when the intended model has discrete events with fixed completion delays. Exactness here means exact simulation of the specified stochastic model, not proof that the model is biologically correct.