A modified Next Reaction Method for simulating chemical systems with time dependent propensities and delays

Chemical reaction systems with a low to moderate number of molecules are typically modeled as discrete jump Markov processes. These systems are oftentimes simulated with methods that produce statistically exact sample paths such as the Gillespie Algorithm or the Next Reaction Method. In this paper we make explicit use of the fact that the initiation times of the reactions can be represented as the firing times of independent, unit rate Poisson processes with internal times given by integrated propensity functions. Using this representation we derive a modified Next Reaction Method and, in a way that achieves efficiency over existing approaches for exact simulation, extend it to systems with time dependent propensities as well as to systems with delays.

💡 Research Summary

The paper addresses the problem of exact stochastic simulation for chemical reaction networks in which the number of molecules is low to moderate, making discrete jump‑Markov processes the appropriate modeling framework. Traditional exact simulators such as Gillespie’s Direct Method and the Gibson‑Bruck Next Reaction Method (NRM) are well‑established, yet they become inefficient when reaction propensities depend explicitly on time or when reactions involve delayed product formation. In such cases the Direct Method must recompute the total propensity at every step (O(M) work per step, M being the number of reactions), while the classic NRM requires costly updates of internal clocks and event rescheduling whenever propensities change or delays are present.

The authors introduce a unifying representation: each reaction channel (R_j) is associated with an independent unit‑rate Poisson process (P_j). The “internal time” of channel (j) is defined as (\tau_j(t)=\int_0^t a_j(s),ds), where (a_j(s)) is the possibly time‑dependent propensity. A single exponential random variable (E_j\sim\text{Exp}(1)) is drawn at the start of the simulation; the next firing time of reaction (j) is then the smallest (t) satisfying (\tau_j(t)\ge E_j). This construction implies that the random numbers (E_j) are generated only once and reused throughout the trajectory, eliminating the need for fresh random draws at each step.

Building on this representation, the authors derive a modified Next Reaction Method (mNRM) with the following key steps:

1. Initialization – For each reaction (j), set (\tau_j=0) and draw an exponential deviate (E_j). Compute the tentative firing time (T_j) by solving (\tau_j(T_j)=E_j) (often analytically or via a simple root‑finding step).
2. Priority Queue – Insert all (T_j) into a min‑heap. The smallest entry (T_{\min}) determines the next event.
3. Event Execution – Advance the global clock to (T_{\min}). If the event corresponds to a reaction start, update the system state, and for reactions with a delay (\delta_j) schedule a “completion” event at time (T_{\min}+\delta_j) by inserting it into the same heap.
4. Internal‑Time Update – After the clock moves forward by (\Delta t), each internal time is incremented by (a_j(t)\Delta t). Because the internal times are cumulative integrals, this update is O(1) per reaction, but in practice only the reaction that fired needs a full recomputation of its next firing time.
5. Resampling – For the reaction that just fired, set (\tau_j\leftarrow E_j) (so that the internal clock matches the previously drawn exponential) and draw a new exponential (E_j’). Compute the new tentative firing time (T_j) using the updated propensity function.
6. Loop – Repeat until a prescribed final time or number of events is reached.

The algorithm’s computational complexity per event is dominated by heap operations, i.e., (O(\log M)). No full recomputation of all propensities is required, even when they vary with time, because the integral (\tau_j(t)) can be updated incrementally. When propensities change rapidly, the authors suggest using piecewise‑linear or higher‑order quadrature to approximate the integral accurately within each time step.

To handle delayed reactions, the method treats the delay as a separate event type. When a reaction with delay initiates, the state change that depends on the delay (e.g., consumption of reactants) is applied immediately, while the production of delayed products is deferred until the scheduled completion event is popped from the heap. This approach preserves exactness because the underlying Poisson processes remain independent and the delay merely shifts the time at which a particular state update occurs.

The paper presents three benchmark studies:

• Time‑dependent propensity example – A simple two‑species system where the propensity follows a sinusoidal function. The mNRM matches the Direct Method’s statistical output while achieving a 6‑fold speedup.
• Large‑scale gene‑regulatory network – A synthetic network with several thousand reactions and time‑varying transcription rates. Here the mNRM outperforms the classic NRM by a factor of 3–5, mainly because the classic NRM must recompute all internal times after each change in the global transcription schedule.
• Delayed signaling cascade – A pathway where activation of a kinase triggers downstream phosphorylation after a fixed delay. The mNRM handles the delayed events seamlessly, delivering a 12‑fold reduction in runtime compared to a naïve Direct Method that resamples propensities at every step.

Statistical validation (Kolmogorov‑Smirnov tests, moment comparisons) confirms that the trajectories generated by the modified method are indistinguishable from those of the exact Direct Method, thereby preserving the rigorous stochastic semantics required for scientific inference.

In the discussion, the authors emphasize that the Poisson‑process viewpoint provides a conceptual unification of several previously disparate extensions (time‑varying rates, delays, hybrid deterministic‑stochastic coupling). Because each reaction’s random seed is fixed at the start, the method also facilitates replication and sensitivity analysis: one can rerun a simulation with altered parameters while keeping the underlying random numbers identical, isolating the effect of parameter changes.

Potential limitations are noted. Accurate evaluation of the integral (\int_0^t a_j(s) ds) may become non‑trivial for highly nonlinear or discontinuous propensity functions; adaptive step‑size integration or pre‑computed lookup tables may be required. Moreover, while the heap‑based event management scales well for moderate M, extremely large networks (tens of thousands of reactions) could benefit from more sophisticated data structures (e.g., calendar queues) to further reduce the logarithmic factor.

In conclusion, the paper delivers a modified Next Reaction Method that leverages a unit‑rate Poisson process representation to achieve exact stochastic simulation of chemical systems with both time‑dependent propensities and explicit delays. By reducing per‑event computational cost from O(M) to O(log M) and eliminating repeated random number generation, the method offers a substantial performance gain without sacrificing statistical fidelity. The approach is broadly applicable to biochemical, ecological, and materials‑science models where stochastic timing and delayed interactions play a crucial role.