PaperlayerView exact revision on arXiv

Source-linked AI summary

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

David F. Anderson

arXiv:0708.0370v2q-bio.MNq-bio.QM

TL;DR

Exact simulation of low- to moderate-molecule chemical systems needs methods for statistically exact sample paths, including systems with time-dependent propensities and delays. The paper uses independent unit-rate Poisson processes with integrated-propensity internal times to derive a modified Next Reaction Method and extend it to these settings. The resulting delayed-system method is reported as computationally more efficient than the Rejection Method, while avoiding wasted random numbers and complicated delayed-reaction machinery.

  • Problem

    Existing exact-simulation methods must be extended to chemical systems with propensities that depend explicitly on time and with delays between reaction initiation and completion.

  • Method

    The paper represents reaction initiations as firings of independent unit-rate Poisson processes with integrated propensity functions as internal times, then derives a modified Next Reaction Method.

  • Results

    The method is extended to systems with delays and reported as computationally efficient over the Rejection Method of Bratsun et al.

  • Takeaways & Limitations

    The modified Next Reaction Method provides a single simulation framework for time-dependent propensities and delayed reactions.

Abstract

from arXiv · show

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.

I. INTRODUCTION

The paper develops an exact-simulation framework for discrete stochastic chemical systems by deriving a modified Next Reaction Method from Poisson-process representations. It extends this approach to time-dependent propensities and delayed reactions, with claimed efficiency benefits over existing methods.

  • Motivation: The paper targets discrete stochastic chemical systems whose state distributions are generally impractical to solve analytically.Simulation methods are therefore used to generate statistically exact sample paths that approximate those distributions.
  • Approach: The authors represent reaction times as firings of independent unit-rate Poisson processes whose internal times are integrated propensity functions.They identify this as a random time change representation and use it to derive new simulation methods.
  • Contributions: The representation naturally yields a modified Next Reaction Method for chemical systems.The paper presents this as the first major consequence of the representation.
  • Contributions: The modified method is designed for propensities that depend explicitly on time, such as changes caused by variable temperature or cell volume.The paper concludes that it is the preferable algorithm for such cases.
  • Contributions: The method extends to reactions with delays between initiation and completion while avoiding wasted random numbers and complicated machinery for stored delayed reactions.The delayed-reaction algorithm is described as an efficient extension over existing approaches.
  • Connections: The paper also relates its ideas to generalized semi-Markov processes and stochastic Petri nets and notes possible extensions to accurate and efficient approximate tau-leaping methods.These connections situate the proposed framework within broader stochastic-simulation methods.

II. THE GILLESPIE ALGORITHM

The Gillespie Algorithm simulates statistically exact reaction paths by repeatedly determining the next reaction time and reaction identity from propensity functions. Each step uses two independent uniform random numbers to perform these selections.

  • Setup: The algorithm models a system of chemical species undergoing reactions, with each reaction assigned a propensity function or intensity.The system state records the number of molecules of each species.
  • Reaction sampling: Under the fundamental premise of chemical kinetics, the next reaction time is exponentially distributed with parameter a0(X(t)).The probability that reaction k is next is ak(X(t))/a0(X(t)).
  • Algorithm: After advancing time by the sampled interval, the algorithm updates molecular counts according to the selected reaction and repeats.The propensity functions are recalculated before the next iteration.
  • Algorithm: The Gillespie Algorithm generates two independent uniform(0,1) random numbers at every step.One determines when the next reaction occurs, and the other determines which reaction occurs.
  • Comparison: The paper uses the Gillespie Algorithm as the baseline for developing a Next Reaction Method that generates exact sample paths with fewer random numbers per step.The Next Reaction Method is introduced as requiring one random number per subsequent step.

III. REPRESENTATION USING POISSON PROCESSES

The paper separates system randomness from chemical state by representing reactions through independent unit-rate Poisson processes run at integrated propensity times. This internal-time view reframes simulation as calculating and comparing Poisson firing times.

  • Representation: Reaction times are represented as firing times of independent unit-rate Poisson processes with internal times given by integrated propensity functions.The reaction-count processes and the chemical state can thereby be expressed through this stochastic representation.
  • Representation: All randomness is contained in the independent Poisson processes and separated from the system state.The state changes only when one of these processes changes.
  • Simulation perspective: The relevant simulation task is to calculate each Poisson process’s next firing time and translate it into a chemical reaction time.This provides a common perspective on different simulation algorithms.
  • Internal time: Each reaction has an absolute time frame and its own internal time defined by the integral of its propensity function.The paper calls this quantity the internal time for reaction k.
  • Gillespie reinterpretation: The Gillespie Algorithm finds the next reaction by determining which Poisson process fires first, using the process intensities given by the propensities.Its subsequent steps use the loss-of-memory property of Poisson processes.
  • Next Reaction perspective: The Next Reaction Methods instead calculate each process’s next firing time and select the minimum, retaining internal times rather than invoking loss of memory.This nearly halves the number of random variables required.

IV. A MODIFIED NEXT REACTION METHOD

The modified Next Reaction Method explicitly tracks internal firing times and absolute firing times, producing an algorithm equivalent to the original Next Reaction Method. Its formulation preserves speed while extending smoothly to time-dependent propensities and non-exponential waiting times.

  • Algorithm: At each step, the method advances to the minimum predicted firing time and updates the state, absolute time, internal times, and propensities.The selected reaction is the channel attaining the minimum.
  • Internal-time update: For reactions that do not fire, the next internal firing time remains unchanged even though their absolute firing times change when propensities change.The updated absolute time is obtained from the corresponding internal-time difference.
  • Efficiency: After initialization, the Next Reaction Method requires one random number per timestep instead of the two required by Gillespie’s method.This reduction follows from retaining scheduled firing information for subsequent steps.
  • Modified formulation: The modified method avoids the time conversions required by the original Next Reaction Method by tracking P_k and T_k explicitly.P_k denotes the next firing time in reaction k’s internal time frame.
  • Time-dependent propensities: Algorithms 2 and 3 have the same simulation speeds because they are equivalent, while Algorithm 3 extends smoothly to time-dependent rate constants.The paper contrasts this with the original formulation’s lack of smooth extension.
  • Non-exponential waiting times: The framework can also handle non-exponential waiting times by changing the waiting-time draws in the relevant algorithm steps.This extension follows from replacing exponential waiting times with draws from the correct distribution.

V. TIME DEPENDENT PROPENSITY FUNCTIONS

The modified Next Reaction Method represents reaction initiations through independent Poisson processes with integrated internal times, enabling exact simulation with time-dependent propensities. It reuses stored internal firing times and extends naturally to non-exponential waiting times.

  • Motivation: Time-dependent rate constants arise when chemical-system conditions such as temperature or volume change, making propensity functions explicitly time dependent.The resulting propensities are no longer constant between reactions.
  • Algorithm: For each reaction channel, the next absolute firing time is obtained by solving an integral equation for the required internal-time increment.The method accounts for the fact that the state remains constant during intervals without reactions.
  • Algorithm: Stored internal firing times let the algorithm update the channel that fired while retaining next-firing information for channels that did not fire.For nonfiring channels, the internal time of the next firing remains available without generating another random number.
  • Efficiency: The modified method calculates each nonfiring channel’s next firing without another random draw, while solving its integral is no more difficult than the original method’s calculation.The efficiency gain comes from avoiding the random-variable generation required by the original Next Reaction Method.
  • Conclusions: The authors conclude that the modified method is preferable to the Gillespie Algorithm and Next Reaction Method for propensities depending explicitly on time.It also extends to non-exponential waiting times between reactions by changing only the generation of internal firing times.

VI. SYSTEMS WITH DELAYS

The paper extends reaction-system simulation to delayed reactions by separating initiation from completion and classifying reactions according to when they change the system state. It explicitly excludes more complex delay structures involving random delays or multiple completion times per initiation.

  • Reaction classes: Reactions are partitioned into ND, CD, and ICD according to whether state changes occur at initiation, completion, or both.ND reactions update at initiation; CD reactions update only at completion; ICD reactions update at both times.
  • Delay structure: Delayed systems require separate handling of reaction initiations and completions because the governing assumption concerns initiation times.The paper therefore distinguishes three reaction cases.
  • Reaction classes: ND reactions consume reactants and produce products at initiation, whereas CD reactions make both changes at completion time t + τk.The timing of state updates differs between the two classes.
  • Reaction classes: ICD reactions remove reactants at initiation and add products at completion time t + τk.This class combines initiation-time and completion-time state updates.
  • Scope: The paper does not consider random delay times or systems with multiple completion times for a single initiation.The authors state that such extensions are straightforward but omit them for clarity.

A. Current Algorithms

Existing delayed-reaction algorithms trade discarded random draws against additional bookkeeping. The Rejection Method is simpler but may discard time steps, whereas Cai’s Direct Method avoids discarded random variables through more elaborate updates during the search for the next initiation.

  • Algorithm requirements: Delayed-reaction simulation must calculate initiation times and store completion times because completions can occur between initiations.Stored delayed reactions can alter the state and propensity functions during simulation.
  • Rejection Method: The Rejection Method computes a timestep, discards it when a delayed completion occurs first, then restarts from the updated state.Its sample paths remain statistically exact despite the discarded timesteps because Poisson-process initiations have the loss-of-memory property.
  • Direct Method: Cai’s Direct Method uses one random number for each initiation and updates the state and propensities during the search across stored completion intervals.It maintains delayed completion times in an ordered d × 2 matrix Tstruct.
  • Trade-offs: The Direct Method avoids discarded random variables but adds machinery for processing delayed completions while finding the next initiation time.The authors note that this extra machinery can slow it relative to the Rejection Method, making the faster method system dependent.

B. The modified Next Reaction Method for systems with delays

The delayed-system extension stores completion times alongside the modified Next Reaction Method’s internal Poisson-process times, allowing initiations and completions to be simulated together. After initialization, it uses one random variable per initiation while retaining extensions to time-dependent and non-exponential systems.

  • Core representation: The delayed extension assigns each delayed reaction channel a vector of completion times stored in ascending order.The next state-changing event is then either a reaction initiation or a stored completion.
  • Algorithm: The algorithm initializes internal times Pk and Tk for every channel and sets each delayed channel’s completion vector to [∞].It then calculates propensities and uses the stored internal firing targets to select the next event.
  • Event handling: When a delayed completion is selected, the method updates the system, removes the earliest completion, and recalculates propensities.When a delayed reaction initiates, it updates the system and inserts its future completion time into the channel’s ordered vector.
  • Efficiency: After initialization, the delayed Next Reaction Method generates one random variable per initiation instead of the two used by the Direct Method.It uses every calculated random variable without the Direct Method’s more complicated machinery.
  • Efficiency: The authors expect the delayed Next Reaction Method to require less simulation time than either the Rejection or Direct Method.This expectation is based on its random-variable reuse and simpler event-update machinery.
  • Extensions: The delayed algorithm also extends easily to time-dependent rate constants and non-exponential waiting times between initiations.The initiation-time construction remains based on independent Poisson processes, with only the internal-time updates changing.

C. Numerical examples

The numerical examples compare Rejection, Direct, and modified Next Reaction methods for delayed chemical systems. Across the reported experiments, the Next Reaction Method is substantially more efficient, including a case where it takes 82.8 CPU seconds versus 179.5 and 167.2 seconds.

  • Example 1: The Direct Method was only 7.4% faster than Rejection in the reported implementation, despite Rejection discarding nearly half its generated random numbers.The comparison indicates that Direct-method calculation steps were not negligible relative to random-number generation.
  • Example 2: As the fraction of delayed completions increases, the Rejection Method becomes more wasteful because rejected time increments increase with delayed reaction initiations.The experiment varied k1 to control γ, the fraction of initiations with delayed completions.
  • Example 2: Across the k1 sweep, Rejection and Direct remained nearly equivalent, while the delayed Next Reaction Method was significantly more efficient for every k1.Rejection was slightly more efficient at smaller k1 and slightly less efficient at larger k1.

VII. CONCLUSION

The paper represents reaction initiation times through independent unit-rate Poisson processes with internal times determined by integrated propensities, yielding a modified Next Reaction Method. It extends this method to time-dependent propensities and delayed systems, where the reported simulations show computational efficiency over existing approaches.

  • Method: The modified Next Reaction Method represents reaction initiation times as firing times of independent unit-rate Poisson processes with integrated propensity functions as internal times.This representation is the basis for the modified algorithm.
  • Delayed systems: The method was extended to systems with delays and demonstrated computational efficiency over the Rejection and Direct Methods for such systems.The comparison names the Rejection Method of Barrio and Bratsun and the Direct Method of Cai.
  • Extensions: The method can also be extended to systems with non-exponential waiting times between reaction initiations.The paper states that this extension is straightforward.
  • Time-dependent propensities: The modified method is presented as a natural choice for systems whose propensities depend explicitly on time.The conclusion contrasts it with both the Gillespie Algorithm and the original Next Reaction Method in this setting.
  • Implications: The authors argue that one efficient method covering this broad range of chemical systems will be beneficial, including models with initiation-to-completion delays such as transcription and translation.They note that the Rejection Method is widely used for such delayed systems.

APPENDIX A: UNFINISHED CALCULATION

The appendix derives the distribution of the next reaction time when propensities depend explicitly on time. It then uses independence of the underlying Poisson processes to identify the reaction selected at that time.

  • Waiting-time distribution: When propensities depend explicitly on time, the absolute waiting time Δ until the next reaction is characterized by a distribution function involving the integrated propensity.The derivation considers each reaction channel's waiting time and the minimum across channels.
  • Reaction selection: Conditioning on F = Δ and using independence of the underlying Poisson processes yields the probability that reaction k fires at that time.The appendix begins this calculation after defining the minimum waiting time.
  • Waiting-time distribution: For each channel j, Δt_j denotes the time until that reaction fires, and F is defined as the minimum of these channel-specific waiting times.The next reaction occurs at F = Δ.
  • Derivation: The appendix combines equations (A1) and (A2), applying L’Hôpital’s rule to obtain the desired result.The cited passage describes the final algebraic step rather than restating the complete expression.

Figures

Figure 1 compares simulation speeds for three reaction methods with delays as k1 and the rejection rate increase. The delayed Next Reaction Method is significantly more efficient than the Rejection and Direct Methods.

  • Figure 1 compares the Rejection, Direct, and delayed Next Reaction Methods as the percentage of rejected timesteps increases with k1.Each method simulated system (18) 10^4 times, and CPU time was plotted against k1.
  • The Rejection and Direct Methods have nearly equivalent CPU times across the tested k1 values.
  • The delayed Next Reaction Method is significantly more efficient than both competing methods for all k1.
Loading 0708.0370v2…