^{1,a)}and Ivo F. Sbalzarini

^{1,b)}

### Abstract

Several real-world systems, such as gene expression networks in biological cells, contain coupled chemical reactions with a time delay between reaction initiation and completion. The non-Markovian kinetics of such reactionnetworks can be exactly simulated using the delay stochastic simulation algorithm (dSSA). The computational cost of dSSA scales with the total number of reactions in the network. We reduce this cost to scale at most with the smaller number of species by using the concept of partial reaction propensities. The resulting delay partial-propensity direct method (dPDM) is an exact dSSA formulation for well-stirred systems of coupled chemical reactions with delays. We detail dPDM and present a theoretical analysis of its computational cost. Furthermore, we demonstrate the implications of the theoretical cost analysis in two prototypical benchmark applications. The dPDM formulation is shown to be particularly efficient for strongly coupled reactionnetworks, where the number of reactions is much larger than the number of species.

R.R. was financed by a grant from the Swiss SystemsX.ch initiative, grant WingX, evaluated by the Swiss National Science Foundation. This project was also supported with a grant from the Swiss SystemsX.ch initiative, grant LipidX-2008/011 to I.F.S.

I. INTRODUCTION

II. THE DELAY STOCHASTIC SIMULATION ALGORITHM (dSSA)

A. The delay direct method

III. THE DELAY PARTIAL-PROPENSITY DIRECT METHOD (dPDM)

A. Prerequisites for dPDM

1. Partial propensities

2. Partial-propensity SSAs

B. Detailed description of the dPDM algorithm

IV. BENCHMARKS

A. A strongly coupled reactionnetwork:Colloidalaggregation model

B. A weakly coupled reactionnetwork: Cyclic chain model

V. CONCLUSIONS AND DISCUSSION

### Key Topics

- Chemical reactions
- 31.0
- Chemical reaction theory
- 7.0
- Chemical kinetics
- 6.0
- Reaction kinetics modeling
- 5.0
- Aggregation
- 4.0

## Figures

Illustration of the main steps in dPDM. (a) Illustration of the linear search to find the interval *p* such that the global time of firing (initiation) of the next reaction . In this figure, the number of pending reactions Δ = 6. (b) Illustration of the partial-propensity structure and the grouping based on the index of the common factored-out reactant. The group index *I* of the next reaction is sampled using linear search over the total propensities of the groups, . The element index *J* within the selected group is found using linear search over the partial propensities stored in group *I*.

Illustration of the main steps in dPDM. (a) Illustration of the linear search to find the interval *p* such that the global time of firing (initiation) of the next reaction . In this figure, the number of pending reactions Δ = 6. (b) Illustration of the partial-propensity structure and the grouping based on the index of the common factored-out reactant. The group index *I* of the next reaction is sampled using linear search over the total propensities of the groups, . The element index *J* within the selected group is found using linear search over the partial propensities stored in group *I*.

Computational cost of dPDM (squares) and dDM (circles). The average (over 100 independent runs) CPU time Θ per reaction initiation (i.e., the average time to execute steps 2–7 in Table II for dPDM, and Table I for dDM) is shown as a function of the number of species *N* in the reaction network. (a) Logarithmic plot of Θ(*N*) for the strongly coupled colloidal aggregation model, considering systems of size up to *N* = 320. Θ is *O*(*N*) for dPDM and for dDM. (b) Linear plot of Θ(*N*) for the strongly coupled colloidal aggregation model, considering systems of size up to *N* = 2000 (2 million reactions). While the scaling of the computational cost remains linear for all system sizes tested, the slope increases around *N* = 1000. This is the system size beyond which the partial-propensity structure does not fit into the computer's cache memory any more. (c) Linear plot of Θ(*N*) for the weakly coupled cyclic chain model. The solid lines are linear least square fits. Θ is *O*(*N*) for both dPDM and dDM, but with a smaller slope for dPDM.

Computational cost of dPDM (squares) and dDM (circles). The average (over 100 independent runs) CPU time Θ per reaction initiation (i.e., the average time to execute steps 2–7 in Table II for dPDM, and Table I for dDM) is shown as a function of the number of species *N* in the reaction network. (a) Logarithmic plot of Θ(*N*) for the strongly coupled colloidal aggregation model, considering systems of size up to *N* = 320. Θ is *O*(*N*) for dPDM and for dDM. (b) Linear plot of Θ(*N*) for the strongly coupled colloidal aggregation model, considering systems of size up to *N* = 2000 (2 million reactions). While the scaling of the computational cost remains linear for all system sizes tested, the slope increases around *N* = 1000. This is the system size beyond which the partial-propensity structure does not fit into the computer's cache memory any more. (c) Linear plot of Θ(*N*) for the weakly coupled cyclic chain model. The solid lines are linear least square fits. Θ is *O*(*N*) for both dPDM and dDM, but with a smaller slope for dPDM.

## Tables

Outline of the algorithm for the delay direct method (dDM) with global times. We only give the main steps and refer to the original publication for the detailed substeps (Ref. 13).

Outline of the algorithm for the delay direct method (dDM) with global times. We only give the main steps and refer to the original publication for the detailed substeps (Ref. 13).

Detailed algorithm for the delay partial-propensity direct method (dPDM), explicitly describing all substeps.

Detailed algorithm for the delay partial-propensity direct method (dPDM), explicitly describing all substeps.

Article metrics loading...

Full text loading...

Commenting has been disabled for this content