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

^{1,b)}

### Abstract

Stochastic reaction-diffusion systems frequently exhibit behavior that is not predicted by deterministic simulation models. Stochastic simulation methods, however, are computationally expensive. We present a more efficient stochastic reaction-diffusion simulation algorithm that samples realizations from the exact solution of the reaction-diffusion master equation. The present algorithm, called partial-propensity stochastic reaction-diffusion (PSRD) method, uses an on-lattice discretization of the reaction-diffusion system and relies on partial-propensity methods for computational efficiency. We describe the algorithm in detail, provide a theoretical analysis of its computational cost, and demonstrate its computational performance in benchmarks. We then illustrate the application of PSRD to two- and three-dimensional pattern-forming Gray-Scott systems, highlighting the role of intrinsic noise in these systems.

R.R. was funded 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. ON-LATTICE STOCHASTIC REACTION-DIFFUSION

A. General concept

B. Discretization-corrected propensities

C. The next subvolume method (NSM) for on-lattice stochastic reaction-diffusion simulations

III. THE PARTIAL-PROPENSITY STOCHASTIC REACTION-DIFFUSION METHOD

A. General concept of PSRD

1. Composition-rejection sampling to select the subvolume

2. Partial propensities to sample the next reaction within a subvolume

B. Detailed description of the PSRD algorithm

1. Data structures

2. Algorithms

IV. BENCHMARKS

A. Colloidal aggregationmodel

B. Linear chain model

V. TWO- AND THREE-DIMENSIONAL SRD SIMULATIONS USING PSRD

VI. CONCLUSIONS AND DISCUSSION

### Key Topics

- Chemical reactions
- 42.0
- Diffusion
- 20.0
- Hydrogen reactions
- 11.0
- Boundary value problems
- 10.0
- Numerical modeling
- 9.0

## Figures

Partitioning of the computation domain into subvolumes. (A–C) Different possibilities of a subdiving a box-shaped computational domain in one, two, and three dimensions, respectively. *L* _{ x }, *L* _{ y }, and *L* _{ z } are the edge lengths of the computational domain in each direction. *K* _{ x }, *K* _{ y }, and *K* _{ z } are the numbers of subvolumes of edge length *h* in each direction. (D) Diffusion is modeled as jump “reactions” to face-connected subvolumes. The same chemical in different subvolumes is treated as a different species. Unimolecular “diffusion reactions” convert species as shown.

Partitioning of the computation domain into subvolumes. (A–C) Different possibilities of a subdiving a box-shaped computational domain in one, two, and three dimensions, respectively. *L* _{ x }, *L* _{ y }, and *L* _{ z } are the edge lengths of the computational domain in each direction. *K* _{ x }, *K* _{ y }, and *K* _{ z } are the numbers of subvolumes of edge length *h* in each direction. (D) Diffusion is modeled as jump “reactions” to face-connected subvolumes. The same chemical in different subvolumes is treated as a different species. Unimolecular “diffusion reactions” convert species as shown.

Illustration of the binning of the total propensities of the subvolumes used for composition-rejection sampling of the next subvolume. The illustration shows a computational domain divided into 4 subvolumes. Points A and B refer to the example in main text used to explain rejection sampling.

Illustration of the binning of the total propensities of the subvolumes used for composition-rejection sampling of the next subvolume. The illustration shows a computational domain divided into 4 subvolumes. Points A and B refer to the example in main text used to explain rejection sampling.

The data structures in PSRD. The contents of the data structures shown corresponds to the example reaction network in Eq. (1) with 3 species and 4 reactions. We assume that the computational domain is divided into 4 subvolumes. In the illustration, , , and are the lumped specific probability rates of the “diffusion reactions” of species 1, 2, and 3 respectively. See main text for details.

The data structures in PSRD. The contents of the data structures shown corresponds to the example reaction network in Eq. (1) with 3 species and 4 reactions. We assume that the computational domain is divided into 4 subvolumes. In the illustration, , , and are the lumped specific probability rates of the “diffusion reactions” of species 1, 2, and 3 respectively. See main text for details.

Computational cost of PSRD and NSM for the aggregation model (Eq. (15)). (A) Computational cost Θ of PSRD (squares) and NSM (circles) as a function of the number of subvolumes *N* _{v} with the size of the reaction network fixed to *N* = 10 (filled symbols) and *N* = 100 (empty symbols), respectively. The solid lines show the corresponding least-squares fits of the theoretical cost models: For *N* = 10, Θ^{PSRD} ≈ 0.02861 log *N* _{v} + 0.03925 *N*, Θ^{NSM} ≈ 0.1095 log *N* _{v} + 0.00581 *f* _{r} *M* + 0.00481(1 − *f* _{r})6*N*; for *N* = 100, Θ^{PSRD} ≈ 0.04401 log *N* _{v} + 0.003579 *N*, Θ^{NSM} ≈ 0.288 log *N* _{v} + 0.001375 *f* _{r} *M* + 0.001418(1 − *f* _{r})6*N*. We estimate for *N* = 10 and for *N* = 100. (B) Computational cost Θ of PSRD (squares) and NSM (circles) as a function of the number of species *N* in the reaction network with the number of subvolumes fixed to *N* _{v} = 512 (filled symbols) and *N* _{v} = 1000 (empty symbols), respectively. The solid lines show the corresponding least-squares fits of the theoretical cost models: For *N* _{v} = 512, Θ^{PSRD} ≈ 0.07559 log *N* _{v} + 0.002258 *N*, Θ^{NSM} ≈ 0.1356 log *N* _{v} + 0.002784 *f* _{r}(⌊*N* ^{2}/2⌋ + *N* + 1) + 0.002633(1 − *f* _{r})6*N*; for *N* _{v} = 1000, Θ^{PSRD} ≈ 0.07205 log *N* _{v} + 0.002777 *N*, Θ^{NSM} ≈ 0.1198 log *N* _{v} + 0.002762 *f* _{r}(⌊*N* ^{2}/2⌋ + *N* + 1) + 0.003163(1 − *f* _{r})6*N*. The fraction *f* _{r} = 0.04 for *N* _{v} = 512 and *f* _{r} = 0.02 for *N* _{v} = 1000.

Computational cost of PSRD and NSM for the aggregation model (Eq. (15)). (A) Computational cost Θ of PSRD (squares) and NSM (circles) as a function of the number of subvolumes *N* _{v} with the size of the reaction network fixed to *N* = 10 (filled symbols) and *N* = 100 (empty symbols), respectively. The solid lines show the corresponding least-squares fits of the theoretical cost models: For *N* = 10, Θ^{PSRD} ≈ 0.02861 log *N* _{v} + 0.03925 *N*, Θ^{NSM} ≈ 0.1095 log *N* _{v} + 0.00581 *f* _{r} *M* + 0.00481(1 − *f* _{r})6*N*; for *N* = 100, Θ^{PSRD} ≈ 0.04401 log *N* _{v} + 0.003579 *N*, Θ^{NSM} ≈ 0.288 log *N* _{v} + 0.001375 *f* _{r} *M* + 0.001418(1 − *f* _{r})6*N*. We estimate for *N* = 10 and for *N* = 100. (B) Computational cost Θ of PSRD (squares) and NSM (circles) as a function of the number of species *N* in the reaction network with the number of subvolumes fixed to *N* _{v} = 512 (filled symbols) and *N* _{v} = 1000 (empty symbols), respectively. The solid lines show the corresponding least-squares fits of the theoretical cost models: For *N* _{v} = 512, Θ^{PSRD} ≈ 0.07559 log *N* _{v} + 0.002258 *N*, Θ^{NSM} ≈ 0.1356 log *N* _{v} + 0.002784 *f* _{r}(⌊*N* ^{2}/2⌋ + *N* + 1) + 0.002633(1 − *f* _{r})6*N*; for *N* _{v} = 1000, Θ^{PSRD} ≈ 0.07205 log *N* _{v} + 0.002777 *N*, Θ^{NSM} ≈ 0.1198 log *N* _{v} + 0.002762 *f* _{r}(⌊*N* ^{2}/2⌋ + *N* + 1) + 0.003163(1 − *f* _{r})6*N*. The fraction *f* _{r} = 0.04 for *N* _{v} = 512 and *f* _{r} = 0.02 for *N* _{v} = 1000.

Computational cost of PSRD and NSM for the linear chain model (Eq. (17)). (A) Computational cost Θ of PSRD (squares) and NSM (circles) as a function of the number of subvolumes *N* _{v} with the size of the reaction network fixed to *N* = 10 (filled symbols) and *N* = 100 (empty symbols), respectively. The solid lines show the corresponding least-squares fits of the theoretical cost models: For *N* = 10, Θ^{PSRD} ≈ 0.03312 log *N* _{v} + 0.03703*N*, Θ^{NSM} ≈ 0.08256 log *N* _{v} + 0.02504 *f* _{r} *M* + 0.002615(1 − *f* _{r})6*N*; for *N* = 100, Θ^{PSRD} ≈ 0.04842 log *N* _{v} + 0.003786*N*, Θ^{NSM} ≈ 0.1428 log *N* _{v} + 0.002934 *f* _{r} *M* + 0.0001978(1 − *f* _{r})6*N* for *N* _{v} ⪅ 512 and Θ^{PSRD} ≈ 0.2923 log *N* _{v} − 0.01199*N*, Θ^{NSM} ≈ 0.5929 log *N* _{v} + 0.0000008 *f* _{r} *M* − 0.004924(1 − *f* _{r})6*N* for *N* _{v} ⪆ 512. We estimate for *N* = 10 and for *N* = 100. (B) Computational cost Θ of PSRD (squares) and NSM (circles) as a function of the number of species *N* in the reaction network with the number of subvolumes fixed to *N* _{v} = 512 (filled symbols) and *N* _{v} = 1728 (empty symbols), respectively. The solid lines show the corresponding least-squares fits of the theoretical cost models: For *N* _{v} = 512, Θ^{PSRD} ≈ 0.03051 log *N* + 0.5291, Θ^{NSM} ≈ 0.07885 log *N* + 0.5458; for *N* _{v} = 1728, Θ^{PSRD} ≈ 0.08479 log *N* + 0.5073, Θ^{NSM} ≈ 0.1642 log *N* + 0.5561. The fraction *f* _{r} = 0.06 for *N* _{v} = 512 and *f* _{r} = 0.03 for *N* _{v} = 1728.

Computational cost of PSRD and NSM for the linear chain model (Eq. (17)). (A) Computational cost Θ of PSRD (squares) and NSM (circles) as a function of the number of subvolumes *N* _{v} with the size of the reaction network fixed to *N* = 10 (filled symbols) and *N* = 100 (empty symbols), respectively. The solid lines show the corresponding least-squares fits of the theoretical cost models: For *N* = 10, Θ^{PSRD} ≈ 0.03312 log *N* _{v} + 0.03703*N*, Θ^{NSM} ≈ 0.08256 log *N* _{v} + 0.02504 *f* _{r} *M* + 0.002615(1 − *f* _{r})6*N*; for *N* = 100, Θ^{PSRD} ≈ 0.04842 log *N* _{v} + 0.003786*N*, Θ^{NSM} ≈ 0.1428 log *N* _{v} + 0.002934 *f* _{r} *M* + 0.0001978(1 − *f* _{r})6*N* for *N* _{v} ⪅ 512 and Θ^{PSRD} ≈ 0.2923 log *N* _{v} − 0.01199*N*, Θ^{NSM} ≈ 0.5929 log *N* _{v} + 0.0000008 *f* _{r} *M* − 0.004924(1 − *f* _{r})6*N* for *N* _{v} ⪆ 512. We estimate for *N* = 10 and for *N* = 100. (B) Computational cost Θ of PSRD (squares) and NSM (circles) as a function of the number of species *N* in the reaction network with the number of subvolumes fixed to *N* _{v} = 512 (filled symbols) and *N* _{v} = 1728 (empty symbols), respectively. The solid lines show the corresponding least-squares fits of the theoretical cost models: For *N* _{v} = 512, Θ^{PSRD} ≈ 0.03051 log *N* + 0.5291, Θ^{NSM} ≈ 0.07885 log *N* + 0.5458; for *N* _{v} = 1728, Θ^{PSRD} ≈ 0.08479 log *N* + 0.5073, Θ^{NSM} ≈ 0.1642 log *N* + 0.5561. The fraction *f* _{r} = 0.06 for *N* _{v} = 512 and *f* _{r} = 0.03 for *N* _{v} = 1728.

Normalized spatial concentration distribution of species S_{1} in the two-dimensional Gray-Scott reaction-diffusion system (Eq. (20)) for *F* = 0.04, *k* = 0.06, *k* _{1} = 1, and *D* _{1} = 2*D* _{2} = 2 × 10^{−5} in a square computational domain of area 0.64^{2}, divided into *N* _{v} = 64^{2} subvolumes (or subareas) of edge length *h* = 0.01. The concentration in each subvolume is shown as a color ranging from blue (concentration zero) to red (concentration one). (A, B) Concentration distributions, normalized by *u*, obtained using PSRD for *u* = 10^{6} (A) and *u* = 10^{7} (B), respectively. (C) Normalized concentration distribution obtained from a deterministic simulation using the same parameters, simulated using second-order finite differences. All snapshots are taken at final time *t* _{f} = 2000/(*k* _{1} *u* ^{2}).

Normalized spatial concentration distribution of species S_{1} in the two-dimensional Gray-Scott reaction-diffusion system (Eq. (20)) for *F* = 0.04, *k* = 0.06, *k* _{1} = 1, and *D* _{1} = 2*D* _{2} = 2 × 10^{−5} in a square computational domain of area 0.64^{2}, divided into *N* _{v} = 64^{2} subvolumes (or subareas) of edge length *h* = 0.01. The concentration in each subvolume is shown as a color ranging from blue (concentration zero) to red (concentration one). (A, B) Concentration distributions, normalized by *u*, obtained using PSRD for *u* = 10^{6} (A) and *u* = 10^{7} (B), respectively. (C) Normalized concentration distribution obtained from a deterministic simulation using the same parameters, simulated using second-order finite differences. All snapshots are taken at final time *t* _{f} = 2000/(*k* _{1} *u* ^{2}).

Normalized spatial concentration distribution of species S_{1} in the three-dimensional Gray-Scott reaction-diffusion system for *F* = 0.04, *k* = 0.06, *k* _{1} = 1, and *D* _{1} = 2*D* _{2} = 2 × 10^{−5} in a cubic computational domain of volume 0.64^{3}, divided into *N* _{v} = 64^{3} subvolumes of edge length *h* = 0.01. The concentration in each subvolume, normalized by *u* = 10^{8}, is shown as a color ranging from blue (concentration zero) to red (concentration one). The snapshot is taken at final time *t* _{f} = 2000/(*k* _{1} *u* ^{2}).

Normalized spatial concentration distribution of species S_{1} in the three-dimensional Gray-Scott reaction-diffusion system for *F* = 0.04, *k* = 0.06, *k* _{1} = 1, and *D* _{1} = 2*D* _{2} = 2 × 10^{−5} in a cubic computational domain of volume 0.64^{3}, divided into *N* _{v} = 64^{3} subvolumes of edge length *h* = 0.01. The concentration in each subvolume, normalized by *u* = 10^{8}, is shown as a color ranging from blue (concentration zero) to red (concentration one). The snapshot is taken at final time *t* _{f} = 2000/(*k* _{1} *u* ^{2}).

## Tables

The detailed algorithm of PSRD.

The detailed algorithm of PSRD.

Article metrics loading...

Full text loading...

Commenting has been disabled for this content