Modeling Delayed Reactions in Stochastic Cellular Systems

Why the Timing Convention Matters

A delayed reaction has at least two relevant times: when it is initiated and when its effects are completed. The interval between them is not merely a timing detail. If reactants are removed at initiation, other reactions cannot use them during the delay. If molecular counts change only at completion, the model must still say whether those reactants are available to other reactions while the event is pending. Different choices define different stochastic systems and can produce different reaction histories.

The earlier material on cellular continuous-time Markov chains establishes the role of a propensity: conditional on the current state, a reaction channel with propensity \(a_j\) has initiation probability \(a_j\,dt+o(dt)\) over a sufficiently short interval \(dt\). For a delayed channel, this is an initiation propensity. It determines when a reaction begins, not when its products appear.

Propensities for Eligible Reactants

Let \(x_i\) be the number of molecules of species \(i\), and let \(\nu^-_{ij}\) be the number of molecules of species \(i\) required by reaction channel \(j\). In a stochastic mass-action model, a common propensity is proportional to the number of eligible reactant combinations. For a channel whose rate per eligible combination is \(c_j\), define its total channel initiation hazard directly by

where \({x_i\choose \nu^-_{ij}}=0\) if \(x_i<\nu^-_{ij}\). The factor counts unordered selections of the required molecules of each species. For example, \(A\to\cdots\) has propensity \(c_jx_A\), while \(2A\to\cdots\) has propensity \(c_j{x_A\choose 2}\). The units of \(c_j\) are inverse time per eligible combination, so \(a_j(x)\) has units of inverse time.

This formula specifies the channel's total hazard. It does not follow merely from assigning each possible molecular combination a marginal small-interval initiation probability: those events could be dependent. If a model instead constructs the channel from separate clocks for combinations, it must specify a construction such as independent clocks, for which simultaneous initiations have probability \(o(dt)\). Otherwise, define the total channel hazard directly, as above.

The standard hazard interpretation assumes that, conditional on the current state and pending-event information, the chance of one channel initiation in \(dt\) is \(a_j\,dt+o(dt)\), and the chance of two or more initiations is \(o(dt)\). Propensities must be nonnegative and must be zero when the required reactants are unavailable. The delay may be fixed, \(\tau_j\), or sampled from a specified distribution. Unless the model says otherwise, the delay is determined at initiation and does not change in response to later state changes.

Two State-Change Conventions

Write \(\nu^-_j\) for the vector of reactants consumed and \(\nu^+_j\) for the vector of products formed. The conventional net stoichiometric change is \(\nu^+_j-\nu^-_j\), but a delayed channel divides that net effect across time differently under the two interpretations below.

For completion-only reactions with reservations, let \(r_i\) be the number of molecules of species \(i\) reserved by pending events. The available, unreserved count is \(x_i-r_i\), so the propensity uses those counts:

At initiation, a new pending event reserves its required reactants, increasing \(r\) by \(\nu^-_j\); total molecular counts \(x\) do not change. At completion, total counts change by \(-\nu^-_j+\nu^+_j\), and the reservation is released. Thus, the molecules committed to that event are removed exactly once. Channels not requiring a reserved molecule may still use other available molecules.

Deriving and Applying the Channel Hazard

For a mass-action channel, the binomial factors count eligible combinations. With \(x_A=4\) and a reaction requiring two \(A\) molecules, there are \({4\choose 2}=6\) eligible pairs. If the model defines the channel hazard as \(c_j\) per eligible pair, its total propensity is \(6c_j\). More generally, the propensity formula is a direct definition of the total rate at which that channel initiates from its current eligible pool. This avoids relying on an unstated assumption about dependence between possible molecular combinations.

Once a delayed event has initiated, its completion is part of the evolving system. For fixed delay \(\tau_j\), an event initiated at time \(t\) completes at \(t+\tau_j\). For a random delay \(L_j\), the model needs a stated distribution for \(L_j\). Pending events and their scheduled completion times, or an equivalent representation, are needed to determine future state changes. The molecular counts alone generally do not determine the future when pending events exist.

Competing Reactions Reveal the Difference

The convention matters most when other channels can use the reactants during the delay. Consider two channels: delayed \(A\to B\) with initiation rate \(0.30\) per minute per available \(A\), and immediate \(A\to C\) with rate \(0.20\) per minute per available \(A\). There are initially three \(A\) molecules. Both channels are assumed to draw from the same available pool; the delayed channel has a fixed 4-minute delay.

Assumptions, Diagnostics, and Alternatives

A well-defined model states which counts each propensity uses, how reactants are committed, and what happens at completion. If a completion-only model intentionally allows reactants to be used again before completion, then it needs another rule to prevent negative counts or to handle a completion that cannot be fulfilled. Silently allowing an impossible completion is not a harmless approximation; it changes or invalidates the process.

Useful checks include the following:

• Availability: Verify that every propensity is zero when too few uncommitted reactants are available.
• Stoichiometric accounting: Confirm that an initiation-time event removes each reactant once and a completion-only event removes each reserved reactant once, at completion.
• Pending-event consistency: At any time, reservations should equal the reactants committed by all pending completion-only events.
• Units: Each total propensity should have units of inverse time; delay times should have units of time.
• Boundary cases: If the delay is zero, initiation and completion effects occur at the same time, so the convention should reduce to the intended instantaneous reaction. If no reactants are available, no new event may initiate.

The main alternative is not a different propensity formula but a different biological interpretation. Initiation-time consumption is appropriate when reactants are committed or depleted as soon as the process begins. Completion-only count changes are appropriate when the molecules remain present but are committed to an ongoing process; the reservation bookkeeping represents that commitment. An Erlang transit chain or another distributed-delay representation instead introduces intermediate states and reaction steps. Such models can represent delay dispersion through explicit stages, but they are not automatically equivalent to a fixed-delay channel with either convention.

Neither convention is universally superior. The choice should reflect when molecules become unavailable and when products become observable in the biological process being represented. A model may also divide effects differently—for example, modifying one species at initiation and producing another at completion—but it must state those changes explicitly and calculate eligibility from the resulting available state.