A constant-time kinetic Monte Carlo algorithm for simulation of large biochemical reaction networks

Abstract

The time evolution of species concentrations in biochemical reaction networks is often modeled using the stochastic simulation algorithm (SSA) [Gillespie, J. Phys. Chem. 81, 2340 (1977)]. The computational cost of the original SSA scaled linearly with the number of reactions in the network. Gibson and Bruck developed a logarithmic scaling version of the SSA which uses a priority queue or binary tree for more efficient reaction selection [Gibson and Bruck, J. Phys. Chem. A 104, 1876 (2000)]. More generally, this problem is one of dynamic discrete random variate generation which finds many uses in kinetic Monte Carlo and discrete event simulation. We present here a constant-time algorithm, whose cost is independent of the number of reactions, enabled by a slightly more complex underlying data structure. While applicable to kinetic Monte Carlo simulations in general, we describe the algorithm in the context of biochemical simulations and demonstrate its competitive performance on small- and medium-size networks, as well as its superior constant-time performance on very large networks, which are becoming necessary to represent the increasing complexity of biochemical data for pathways that mediate cell function.

An exact accelerated stochastic simulation algorithm

“…Further developments include multiscale SSAs such as "nested stochastic simulation," 6 the multiscale methods, 7,8 and the "slow-scale stochastic simulation" algorithm. 9 Another acceleration method 10 uses rejection sampling to achieve constant time scaling with the number of reaction channels; this differs from the present work which uses rejection sampling to improve scaling with respect to the number of reaction events.…”

A partial-propensity variant of the composition-rejection stochastic simulation algorithm for chemical reaction networks

“…For small networks, PSSA-CR is outperformed by other methods due to the additional overhead involved in the composition-rejection sampling. SSA formulations such as SDM, 12 NRM, 4 SSA-CR, 5 PDM, or SPDM 6 might be more efficient here. In addition, PSSA-CR only achieves the O(1) scaling for weakly coupled networks for which ratio of maximum to minimum non-zero reaction propensity is bounded by a constant throughout a simulation.…”

“…5 It is limited to chemical reaction networks composed of elementary reactions involving at most two reactants. Non-elementary reactions can be treated by decomposing them into elementary reactions.…”

“…We apply the composition-rejection sampling strategy 8,5 to obtain the group index I and the element index J of the next reaction µ. The group index I is sampled in two steps: (1) the composition step to find the bin b I and (2) the rejection step to find Σ I inside that bin.…”

“…The partial-propensity SSA with composition-rejection sampling (PSSA-CR) is based on the idea of factorizing the reaction propensities as described below, 6 grouping and binning them, and using composition-rejection sampling 8,5 to determine the index of the next reaction.…”