^{1}and Russell Schwartz

^{2,a)}

### Abstract

Models of reaction chemistry based on the stochastic simulation algorithm (SSA) have become a crucial tool for simulating complicated biological reactionnetworks due to their ability to handle extremely complicated networks and to represent noise in small-scale chemistry. These methods can, however, become highly inefficient for stiff reaction systems, those in which different reaction channels operate on widely varying time scales. In this paper, we develop two methods for accelerating sampling in SSA models: an exact method and a scheme allowing for sampling accuracy up to any arbitrary error bound. Both methods depend on the analysis of the eigenvalues of continuous time Markovmodels that define the behavior of the SSA. We show how each can be applied to accelerate sampling within known Markovmodels or to subgraphs discovered automatically during execution. We demonstrate these methods for two applications of sampling in stiff SSAs that are important for modelingself-assemblyreactions: sampling breakage times for multiply connected bondnetworks and sampling assembly times for multisubunit nucleation reactions. We show theoretically and empirically that our eigenvalue methods provide substantially reduced sampling times for a large class of models used in simulating self-assembly. These techniques are also likely to have broader use in accelerating SSA models so as to apply them to systems and parameter ranges that are currently computationally intractable.

This work was supported in part by U.S. National Science Foundation Award No. 0346981.

I. INTRODUCTION

II. THEORY

A. The chemical master equation and the SSA

B. Spectral sampling 1: Master equation approach

1. Spectral decomposition of the first-passagetime distribution

2. Exact sampling for the first-passage time distribution

C. Spectral sampling 2: Modified EMC method

D. AD of trapped subgraphs

III. APPLICATION WHEN THE SUBGRAPH IS KNOWN: FRACTURING BONDNETWORKS

A. Stiffness in SSA for bondnetworks

B. Master equation for bondnetworks

1. Spectral analysis for breaking the network

C. Simulation models used for bondnetworks

D. Experiments

E. Results

IV. APPLICATION WITH AD: NUCLEATION-LIMITED ASSEMBLY

A. Integer lattice models

B. AD for integer lattice

C. Experiments

D. Results

V. DISCUSSION

### Key Topics

- Eigenvalues
- 39.0
- Reaction kinetics modeling
- 24.0
- Spectral methods
- 20.0
- Networks
- 13.0
- Self assembly
- 13.0

## Figures

Illustration of trapped subgraphs in SSA models. (a) A simple CTMM on a three-path with transition rates and . (b) The probability landscape for the model. SSA is slow whenever the invariant density for the corresponding Markov chain is irregular. Here, the SSA takes steps to reach vertex 3. (c) CTMM model of a trimer assembly system with three subunits. Graph of possible configurations joined by reaction rates. States in which the trimer is broken are surrounded by solid lines and others by dashed lines.

Illustration of trapped subgraphs in SSA models. (a) A simple CTMM on a three-path with transition rates and . (b) The probability landscape for the model. SSA is slow whenever the invariant density for the corresponding Markov chain is irregular. Here, the SSA takes steps to reach vertex 3. (c) CTMM model of a trimer assembly system with three subunits. Graph of possible configurations joined by reaction rates. States in which the trimer is broken are surrounded by solid lines and others by dashed lines.

Pseudocode for spectral method 1.

Pseudocode for spectral method 1.

Schematic of the EMC-based spectral method. Filled vertices are the currently occupied nodes. Simulation advances the system state as a discrete mixture until such time as the state has relaxed to its slowest eigenstate . At each step, direct transitions to the absorbing vertex (gray) are computed according to the Kolmogorov matrix.

Schematic of the EMC-based spectral method. Filled vertices are the currently occupied nodes. Simulation advances the system state as a discrete mixture until such time as the state has relaxed to its slowest eigenstate . At each step, direct transitions to the absorbing vertex (gray) are computed according to the Kolmogorov matrix.

Pseudocode for spectral method 2.

Pseudocode for spectral method 2.

Pseudocode for AD.

Pseudocode for AD.

Number of SSA steps until first passage for the network generated by an -cycle . (a) Average number of steps . (b) Relative deviation .

Number of SSA steps until first passage for the network generated by an -cycle . (a) Average number of steps . (b) Relative deviation .

Number of SSA steps until first passage on an -dimensional unit hypercube . (a) A plot of the average number of steps vs rate ratio . (b) Relative deviation .

Number of SSA steps until first passage on an -dimensional unit hypercube . (a) A plot of the average number of steps vs rate ratio . (b) Relative deviation .

Number of rejection steps for the master equation method until first passage for the network generated by . (a) Average number of steps . (b) Standard deviation .

Number of rejection steps for the master equation method until first passage for the network generated by . (a) Average number of steps . (b) Standard deviation .

Number of rejection steps for the master equation method until first passage for . (a) Average number of steps . (b) Standard deviation .

Number of rejection steps for the master equation method until first passage for . (a) Average number of steps . (b) Standard deviation .

Number of EMC steps until first passage for the network generated by . (a) Average number of steps . (b) Standard deviation . (c) Fraction of times the trajectory escapes before relaxing to the slowest decay mode.

Number of EMC steps until first passage for the network generated by . (a) Average number of steps . (b) Standard deviation . (c) Fraction of times the trajectory escapes before relaxing to the slowest decay mode.

Number of EMC steps until first passage on . (a) Average number of steps . (b) Standard deviation . (c) Fraction of times the trajectory escapes before relaxing to the slowest decay mode.

Number of EMC steps until first passage on . (a) Average number of steps . (b) Standard deviation . (c) Fraction of times the trajectory escapes before relaxing to the slowest decay mode.

Comparative run times for the network generated by . (a) Ratio of SSA to master equation run times. (b) Ratio of SSA to EMC run times. (c) Region in two dimensional (2D) parameter space where each method is optimal.

Comparative run times for the network generated by . (a) Ratio of SSA to master equation run times. (b) Ratio of SSA to EMC run times. (c) Region in two dimensional (2D) parameter space where each method is optimal.

Comparative run times for first passage on . (a) Ratio of SSA to master equation run times. (b) Ratio of SSA to EMC run times. (c) Region in 2D parameter space where each method is optimal.

Comparative run times for first passage on . (a) Ratio of SSA to master equation run times. (b) Ratio of SSA to EMC run times. (c) Region in 2D parameter space where each method is optimal.

Pseudocode for AD for a simple path.

Pseudocode for AD for a simple path.

(a) Comparative run times for first passage for ten trimer counts. Ratio of SSA to ME-AD run times. (b) Comparative run times for first passage for 100 trimer counts. Ratio of SSA to EMC-AD run times.

(a) Comparative run times for first passage for ten trimer counts. Ratio of SSA to ME-AD run times. (b) Comparative run times for first passage for 100 trimer counts. Ratio of SSA to EMC-AD run times.

Article metrics loading...

Full text loading...

Commenting has been disabled for this content