Using Equation-Free Computation to Accelerate Network-Free Stochastic Simulation of Chemical Kinetics

Abstract

The chemical kinetics of many complex systems can be concisely represented by reaction rules, which can be used to generate reaction events via a kinetic Monte Carlo method that has been termed network-free simulation. Here, we demonstrate accelerated network-free simulation through a novel approach to equation-free computation. In this process, variables are introduced that approximately capture system state. Derivatives of these variables are estimated using short bursts of exact stochastic simulation and finite differencing. The variables are then projected forward in time via a numerical integration scheme, after which a new exact stochastic simulation is initialized and the whole process repeats. The projection step increases efficiency by bypassing the firing of numerous individual reaction events. As we show, the projected variables may be defined as populations of building blocks of chemical species. The maximal number of connected molecules included in these building blocks determines the degree of approximation. Equation-free acceleration of network-free simulation is found to be both accurate and efficient.

Figure 1

Illustration of the basic concepts of equation-free computation through simulation of stochastic logistic growth. We consider a model wherein the possible reactions are birth events N → N + 1, with a rate βN, and death events N → N − 1, with rate δN², where N is the population. We set β = 2 s⁻¹ and δ = 0.001 s⁻¹. (a) An equation-free trajectory, which is derived through a combination of stochastic simulation runs (black points) and extrapolation (red line), and the deterministic solution (green broken line). The periods over which the exact simulations were performed are indicated by gray boxes. (b) A magnification of the sample path illustrates the procedure. (c) A schematic diagram of the equation-free method and its four stages. 1. Simulation. Exact stochastic simulation is executed. 2. Restriction. Coarse state variables are calculated based on the system state at the end of exact simulation. 3. Projection. The coarse variables are projected forward in time. 4. Lifting. A detailed system configuration is generated from the coarse state variables for the next exact simulation.

Figure 2

Illustration of a model for the assembly of a linear aggregate. (a) The contact map of the model, which consists of six molecules A, B, C, D, and E considered in the model and their interactions. Each molecule has specific binding site(s), which are responsible for interaction with sites in other molecules. There are 21 possible complexes in this model. (b) Shorthand notation to denote the configurations of a single molecule. (c) First-degree moieties using the short-hand notation. The shaded boxes identify redundant information because populations of these moieties can be deduced from the populations of the unshaded moieties. (d) Second-degree moieties.

Figure 3

Moiety variables for a model for early events in EGFR signaling. (a) The contact map of the model. (b) The first-degree moieties. The moiety variables can be defined as model outputs (i.e., observables) using the BioNetGen language (BNGL). Redundant information involving EGF and SOS molecules is removed, similar to the linear assembly model. (c) The second-degree moieties.

Figure 4

Stochastic simulation of EGFR signaling according to the model of Blinov et al. with equation-free acceleration. For this model, the reaction network and the full list of reacting species can be generated from the model’s rules using BioNetGen. To perform stochastic simulations, we implemented Gillespie’s direct method with equation-free acceleration in problem-specific code. Each burst of stochastic simulation consisted of 4096 reaction events. At each event time, we calculated and recorded the values of the coarse-grained moiety variables. To estimate derivatives for these variables, we performed linear fits to their stochastic time courses. The resulting slopes were taken to be the time derivatives of the variables, which were used in projection via Euler’s method. The time increment used in projection was chosen to be three times the duration of the prior burst of stochastic simulation.

Figure 5

The efficiency of equation-free acceleration of stochastic simulation for the EGFR model with K = 1, 2, 3, and 4 relative to exact stochastic simulation (Gillespie’s direct method). The simulations were independently executed in parallel (32 simultaneous threads) on a machine with 32 Intel Xeon CPUs (E5-2698 v3, 2.30 GHz).