Stochastic simulation of reaction-diffusion systems: A fluctuating-hydrodynamics approach

Kim, Changho; Nonaka, Andy; Bell, John B.; Garcia, Alejandro L.; Donev, Aleksandar

doi:10.1063/1.4978775

We develop numerical methods for stochastic reaction-diffusion systems based on approaches used for fluctuating hydrodynamics (FHD). For hydrodynamic systems, the FHD formulation is formally described by stochastic partial differential equations (SPDEs). In the reaction-diffusion systems we consider, our model becomes similar to the reaction-diffusion master equation (RDME) description when our SPDEs are spatially discretized and reactions are modeled as a source term having Poisson fluctuations. However, unlike the RDME, which becomes prohibitively expensive for an increasing number of molecules, our FHD-based description naturally extends from the regime where fluctuations are strong, i.e., each mesoscopic cell has few (reactive) molecules, to regimes with moderate or weak fluctuations, and ultimately to the deterministic limit. By treating diffusion implicitly, we avoid the severe restriction on time step size that limits all methods based on explicit treatments of diffusion and construct numerical methods that are more efficient than RDME methods, without compromising accuracy. Guided by an analysis of the accuracy of the distribution of steady-state fluctuations for the linearized reaction-diffusion model, we construct several two-stage (predictor-corrector) schemes, where diffusion is treated using a stochastic Crank–Nicolson method, and reactions are handled by the stochastic simulation algorithm of Gillespie or a weakly second-order tau leaping method. We find that an implicit midpoint tau leaping scheme attains second-order weak accuracy in the linearized setting and gives an accurate and stable structure factor for a time step size of an order of magnitude larger than the hopping time scale of diffusing molecules. We study the numerical accuracy of our methods for the Schlögl reaction-diffusion model both in and out of thermodynamic equilibrium. We demonstrate and quantify the importance of thermodynamic fluctuations to the formation of a two-dimensional Turing-like pattern and examine the effect of fluctuations on three-dimensional chemical front propagation. By comparing stochastic simulations to deterministic reaction-diffusion simulations, we show that fluctuations accelerate pattern formation in spatially homogeneous systems and lead to a qualitatively different disordered pattern behind a traveling wave.

I. INTRODUCTION

While deterministic reaction-diffusion models have been successfully applied to explain various spatiotemporal phenomena such as pattern formation, and to gain insight into nonequilibrium transitions, it is now widely appreciated that spatiotemporal fluctuations in the concentration of chemical species play an essential role. Such internal or thermodynamic fluctuations, which arise from both reaction and diffusion processes, have molecular origin; microscopically, those processes occur through the movement and collision of individual molecules under thermal fluctuations. Hence, the deterministic macroscopic description eventually fails at smaller scales where the fluctuations are significant, and a stochastic mesoscopic description is needed. Examples include fluctuation-induced instabilities,¹ reversal of direction of front propagation,² violation of the law of mass action,³ long-time tails in kinetics,⁴ emergence of new steady states⁵ and patterns,⁶ acceleration of pattern formation,⁷ enhanced induction time for ignition,⁸ and the onset of homogeneous oscillations.⁹ Due to a small number of proteins involved in cellular functions,¹⁰ processes in cell biology are good examples^11–14 where the stochastic reaction-diffusion description provides an indispensable modeling tool.^15,16

A microscopic picture of reaction-diffusion, dating back to the work of Smoluchowski,¹⁷ assumes that molecules undergo independent Brownian motions and reactions can occur only when two molecules are close to each other. Based on this picture, the particle-based approach to simulate a reaction-diffusion system tracks the trajectories of diffusing molecules and uses the intermolecular distance to determine whether a reaction occurs. Exact sampling of the Smoluchowski model can be performed by first-passage kinetic Monte Carlo type algorithms;^18–20 approximate reactive Brownian dynamics (BD) using a fixed time step size forms another class of algorithms.^21,22 While molecular schemes, such as molecular dynamics (MD) and direct simulation Monte Carlo (DSMC), can be used for reaction-diffusion problems,^23,24 they are computationally even more expensive. Hybrid methods combining particle and coarse-grained descriptions, either using operator splitting²⁵ or domain decomposition,^26,27 have also been proposed.

For the mesoscopic description of a reactive system, the master equation approach is commonly used. For a well-mixed (i.e., spatially homogeneous) system, the time evolution of the system (i.e., the number of molecules of each chemical species) is described by the chemical master equation (CME). Exact sampling of the CME can be performed by the stochastic simulation algorithm (SSA) of Gillespie,²⁸ whereas the tau leaping method²⁹ can be employed as an approximate algorithm with a given time step size. Several variants of these methods have been proposed.^30,31 For a spatially inhomogeneous system, the time evolution of the system is commonly described by the reaction-diffusion master equation (RDME), which is also known as the multivariate master equation.^32,33 In this approach, the system is divided into homogeneous subsystems or cells and the number of molecules of each chemical species in each cell is tracked. Changes in the molecule numbers occur either through hopping events of a molecule between adjacent cells or through chemical reactions within a cell. Hopping events correspond to diffusive transport and are treated as first-order reactions. Since the RDME is a spatial extension of the CME, exact sampling of the RDME can be performed by SSA-type algorithms,^34,35 which are called inhomogeneous SSA (ISSA).

While conceptually simple and still widely used,^6,7,36 the traditional approach of solving the RDME by ISSA has the computational issue that the method becomes prohibitively slow as the number of molecules or cells increases. Since the cell volume should be chosen sufficiently small to ensure homogeneity over each cell, large, finely resolved grids are required for two- or three-dimensional problems. As the spatial resolution increases, the time interval between successive events becomes very short due to rapid diffusive transfer, and hopping events greatly outnumber reaction events, which slows down ISSA.³⁰ Several approaches have been proposed to improve the performance of stochastic sampling of the RDME, such as the next subvolume method³⁵ and its parallel simulation version.³⁷ Various implementations of the tau leaping method in a spatial context,^38–40 and the time-dependent propensity for diffusion method,⁴¹ have also been proposed. A more aggressive approach to reduce the computational cost by avoiding the sampling of the individual diffusion events is to split diffusion and reaction in each time step and to treat diffusion in a more efficient manner. Various sampling methods for diffusion have been proposed, including the Gillespie multi-particle method,⁴² the multinomial simulation algorithm,⁴³ the adaptive hybrid method on unstructured meshes,^44,45 and the diffusive finite-state projection algorithm.^46,47

In this paper, we propose a numerical algorithm for stochastic reaction-diffusion systems based on approaches used for fluctuating hydrodynamics (FHD). To incorporate the effects of thermal fluctuations in a fluid, in FHD one assumes that the dynamics of the fluid can be described by the usual hydrodynamic equations (e.g., the Navier–Stokes equations), augmenting each dissipative flux with a stochastic flux.⁴⁸ Those stochastic fluxes are modeled by spatiotemporal Gaussian white noise (GWN) and the resulting governing equations are written as stochastic partial differential equations (SPDEs). FHD was originally developed for equilibrium fluctuations by Landau and Lifshitz⁴⁹ and its validity has been justified for nonequilibrium systems⁵⁰ through the theory of coarse graining.⁵¹ For further discussion of FHD compared to MD, see Ref. 52. Various extensions and generalizations of FHD theory have been developed and successfully applied to fluctuation-induced phenomena; see Ref. 53 and references therein. Recent work by the authors has focused on FHD models of hydrodynamic transport^54,55 in binary fluid mixtures,^56,57 multiphase flows,⁵⁸ multispecies fluid mixtures,^59,60 multispecies reactive mixtures,⁶¹ and electrolytes.⁶²

Compared to our previous work,⁶¹ where the coupling effects of fluid hydrodynamic and chemical fluctuations have been investigated, here we focus on reaction and diffusion and neglect all other hydrodynamic processes (advection, viscous dissipation, thermal conduction, and cross term effects). Rather than using a Langevin description (i.e., based on Gaussian fluctuations) of chemistry, which is only valid in the limit of vanishing fluctuations,⁶¹ here we employ a more accurate description of reactions based on Poisson fluctuations. As pointed out in Ref. 51, even though a formal SPDE description is employed, the actual interpretation of FHD always requires the notion of a coarse-graining over a certain length scale. The FHD equations are discretized using a finite volume approach^63,64 that represents the solution in terms of the average over cells, which provides an effective coarse-graining. Therefore, reactions can be treated in a similar manner to the RDME approach when the SPDEs are spatially discretized and integrated in time using SSA or a weakly second-order tau leaping method.^65,66 Recent relevant work by others includes Ref. 67 in which the FHD approach has been applied to reaction-diffusion systems. However, only fluctuations arising from diffusion have been considered (i.e., no fluctuations from chemical reactions) and modeled as additive noise. The FHD approach has been also applied to concentration fluctuations in a ternary liquid mixture in equilibrium⁶⁸ and the Model H equations for binary mixtures.⁶⁹

The key difference between the FHD and RDME descriptions lies in the more efficient treatment of fast diffusion. A number of approximate numerical methods for the RDME^42–47 are based on operator splitting using first-order Lie or second-order Strang splitting.⁷⁰ In Appendix A we review and discuss in more detail a split scheme that uses multinomial diffusion sampling⁷¹ for diffusion and SSA for reactions. These RDME-based schemes use a time step size $Δ t$ comparable to the hopping time scale $τ_{h} = Δ x^{2} / (2 d D)$ with d being the spatial dimension, $Δ x$ being the grid spacing, and D being a typical diffusion coefficient. Even though $τ_{h}$ is much larger than the mean duration between successive events in ISSA, using $Δ t$ comparable to $τ_{h}$ is still very restrictive for large D or small $Δ x$ ⁠. In our FHD formulation, we treat diffusion implicitly using backward Euler or Crank–Nicolson, so that the time step size can be significantly larger (e.g., an order of magnitude larger for a given accuracy tolerance) than the hopping time scale. Since the time steps used in RDME simulations are already (usually an order of magnitude or more) larger than those in BD simulations, our approach allows even larger time step size compared to particle-based methods.

While the development of numerical schemes for stochastic reaction-diffusion systems described by spatiotemporal GWN dates back to the 1990s,^72,73 much of the prior work has not been guided by numerical analysis or extensive experience from deterministic computational fluid dynamics (CFD). With the help of well-established techniques for numerical solution of PDEs and SPDEs, we construct numerical schemes in a systematic manner to ensure that accuracy is maintained for a large time step size. To this end, we employ two-stage (i.e., predictor-corrector) Runge–Kutta temporal integrators.^64,74 Rather than using operator splitting, we treat reaction and diffusion together in each stage in a manner that is second-order weakly accurate for general linearized FHD equations. The construction of these schemes is guided by a stochastic accuracy analysis of the (static) structure factor for linearized FHD.^63,64 The structure factor is the steady-state spectrum of the concentration fluctuations, i.e., the covariance matrix in Fourier space, see Eq. (8). We apply the techniques in Refs. 63 and 64 to predict the discrete structure factors for our scheme and compare them to analytical predictions of the continuum structure factors for our model in the linearized setting.

The FHD approach inherently outperforms the RDME approach as the number of molecules per cell increases in exactly the same way that multinomial diffusion outperforms diffusion by hopping, or tau leaping outperforms SSA. In fact, the computational cost of FHD methods does not significantly change as the magnitude of the fluctuations changes. This is an obvious advantage of the FHD approach since the macroscopic limit cannot be efficiently simulated by the RDME approach. However, the validity of the FHD approach cannot be taken for granted when there are only a small number of molecules in each cell, since in FHD the number of molecules in each cell is a continuous real-valued variable, rather than a discrete non-negative integer variable as in the RDME. We investigate this issue carefully and propose techniques to improve the accuracy of the FHD description for the case of a small number of molecules per cell, making our numerical schemes robust even for large fluctuations. In particular, we develop a spatial discretization that significantly mitigates non-negativity of the species number densities and closely reproduces the Poisson thermodynamic equilibrium distribution for the number of molecules in a cell. For numerical examples considered in this paper, we show that the mean number of molecules in a cell can be as low as 10. However, if one is specifically interested in systems with only a small number of molecules per cell, one should use an integer-based description like RDME. Moreover, in very dilute cases, a particle-based description like BD is actually fastest since most cells will have essentially no molecules in them. However, for practical stochastic simulation of reaction-diffusion systems, where the populations of chemical species may have different orders of magnitude, this kind of robustness is required; even if there are a large number of molecules in a cell, some species may have a small number of molecules.

The rest of the paper is organized as follows. Section II presents the background for our approach, including the FHD description of reaction-diffusion systems and linearized analysis in a Gaussian approximation. Section III explains how the FHD reaction-diffusion equations can be spatially discretized using a finite-volume approach. Section IV presents temporal integrators for the spatially discretized equations that handle diffusion using existing FHD techniques and treat reactions using SSA or second-order tau leaping. Section V presents simulation results of several reaction-diffusion systems. In Section V A, for testing and validation of our numerical schemes, we use a one-species Schlögl model.^75,76 In Section V B, to compare our methods to each other and to RDME-based methods, we study two-dimensional Turing-like pattern formation in the three-species Baras–Pearson–Mansour (BPM) model.^77,78 In Section V C, to demonstrate the ability of our approach to scale to larger systems, we present numerical simulation results for three-dimensional front propagation in a two-species model.⁷ In Section VI, we offer some concluding remarks and suggest future research directions.

II. BACKGROUND

In Section II A, we present the continuous-time continuous-space FHD description of reaction-diffusion systems. Here, we assume that fluctuations in chemistry are described by GWN (i.e., Langevin type). A more accurate description of chemistry based on Poisson fluctuations is incorporated in the continuous-time discrete-space description in Section III. In Sections II B and II C, we introduce the structure factor and the Schlögl reaction-diffusion model, respectively. As one of the criteria for the development and analysis of numerical schemes, later in the paper, we investigate how accurately a numerical scheme produces the structure factor for the Schlögl model.

In this section, we introduce several GWN vector and scalar random fields and denote them by $𝓩 (𝒙, t) = (Z_{1} (𝒙, t), \dots, Z_{d} (𝒙, t))$ and $Z (𝒙, t)$ ⁠, respectively, with additional superscripts to distinguish the different fields. We assume that any two distinct processes are independent and that the noise intensity of each process is normalized: $⟨ Z_{j} (𝒙, t) Z_{j^{'}} (𝒙^{'}, t^{'}) ⟩ = δ_{j j^{'}} δ (𝒙 - 𝒙^{'}) δ (t - t^{'})$ and $⟨ Z (𝒙, t) Z (𝒙^{'}, t^{'}) ⟩$ = $δ (𝒙 - 𝒙^{'}) δ (t - t^{'})$ ⁠. Similarly, we denote GWN vector and scalar random processes by $W (t)$ and $W (t)$ ⁠, respectively, and assume $⟨ W_{j} (t) W_{j^{'}} (t^{'}) ⟩ = δ_{j j^{'}} δ (t - t^{'})$ and $⟨ W (t) W (t^{'}) ⟩$ $= δ (t - t^{'})$ ⁠.

A. FHD description

We consider a reaction-diffusion system having N_s species undergoing N_r reactions in d-dimensional space. By denoting the number density of species s by $n_{s} (𝒙, t)$ ⁠, the equations of FHD for $𝒏 (𝒙, t) = (n_{1} (𝒙, t), \dots, n_{N_{s}} (𝒙, t))$ are written formally as the SPDEs⁶¹

\frac{\partial}{\partial t} n_{s} = \nabla \cdot (D_{s} \nabla n_{s} + \sqrt{2 D_{s} n_{s}} (_{s}^{(D)}) + \sum_{r = 1}^{N_{r}} ν_{s r} (a_{r} (𝒏) + \sqrt{a_{r} (𝒏)} Z_{r}^{(R)}),

(1)

where D_s is the diffusion coefficient of species s, $a_{r} (𝒏)$ is the propensity function indicating the rate of reaction r, and $ν_{s r}$ is the stoichiometric coefficient of species s in reaction r. In the macroscopic limit of vanishing fluctuations, Eq. (1) approaches the deterministic reaction-diffusion PDE (law of large numbers),

\frac{\partial}{\partial t} n_{s} = D {}_{s}\nabla^{2} n_{s} + \sum_{r = 1}^{N_{r}} ν_{s r} a_{r} (𝒏) .

(2)

We explain below how the diffusion and reaction parts are obtained by considering the diffusion-only (i.e., no-reaction) and reaction-only (i.e., well-mixed) cases.

1. Diffusion

The diffusion-only SPDE

\frac{\partial}{\partial t} n_{s} = \nabla \cdot (D_{s} \nabla n_{s} + \sqrt{2 D_{s} n_{s}} _{s}^{(D)})

(3)

can be justified by considering a microscopic system where each molecule i undergoes independent Brownian motion,

{\dot{𝒙}}_{s, i} = \sqrt{2 D_{s}} W_{s, i} .

(4)

In this formal derivation,⁷⁹ one defines the instantaneous number density field $n_{s} (𝒙, t) = \sum_{i} δ (𝒙 - 𝒙_{s, i} (t))$ and uses Ito’s rule to obtain Eq. (3). This equation can also be obtained from the diffusion portion of the general multispecies FHD equations,^60,61 given by nonequilibrium statistical mechanics,⁸⁰ by assuming a dilute solution and considering only solute species. In addition, the linearized version of Eq. (3) can be obtained from the multivariate master equation model (i.e., diffusion by hopping) near the macroscopic limit.³² However, although relations to those equations reaffirm Eq. (3) near the macroscopic limit, it is important to note that Eq. (3) is formally exact even in the case where fluctuations are large, since it is simply a rewriting of Eq. (4), in a representation in which the particle numbering (identity) is lost.⁸¹

We note that the FHD equations (1) and (3) are not mathematically well-defined because the solution needs to be interpreted as a distribution (or a generalized function), and the square root of a distribution is not well-defined in general. The linearized FHD equation does not suffer from such an issue and is well-defined; the problems arise due to the multiplicative noise in Eq. (3). However, even though Eq. (3) is ill-defined, it is formally consistent with the law of large numbers (given by the deterministic reaction-diffusion equation (2)), the central limit theorem (given by the linearized FHD equation (12)), and large deviation theory for a collection of Brownian walkers. In this sense, Eq. (3) is a meaningful representation of the physical model that is useful in constructing well-defined mesoscopic descriptions via spatial discretization of the formal SPDEs. Compared to obtaining a mesoscopic model by directly coarse-graining the microscopic model, the spatial discretization of the SPDE is easier in general and can be done in a systematic manner.⁸²

2. Reaction

To see how the reaction part of Eq. (1) is obtained, consider a well-mixed system with volume $Δ V$ ⁠. By assuming that the time evolution of $𝑛 (t)$ follows the CME, we express the change over the infinitesimal time interval dt as follows:³⁰

d n_{s} = n_{s} (t + d t) - n_{s} (t) = \frac{1}{Δ V} \sum_{r = 1}^{N_{r}} ν_{s r} P (a_{r} (𝒏) Δ V d t),

(5)

where $P (m)$ denotes a Poisson random variable having mean m. Note that Eq. (5) is equivalent to the CME if interpreted in the Ito sense. The specific form of the chemical rate function $a_{r} (𝒏)$ that we use in this work is described in Section III C. Henceforth, we will formally write Eq. (5) in the differential form,

\frac{d}{d t} n_{s} = \sum_{r = 1}^{N_{r}} \frac{ν_{s r} P (a_{r} (𝒏) Δ V d t)}{Δ V d t} .

(6)

For a more mathematically precise representation, see Refs. 83 and 84.

The chemical Langevin equation (CLE)⁸⁵ is obtained under the assumption that the mean number of reaction occurrences is large.³⁰ That is, the assumption enables one to replace $P (m)$ by a Gaussian random variable having the same mean and variance, to give the CLE,

\frac{d}{d t} n_{s} = \sum_{r = 1}^{N_{r}} ν_{s r} [a_{r} (𝒏) + \sqrt{\frac{a_{r} (𝒏)}{Δ V}} W_{r}] .

(7)

Since reaction is assumed to be local, the reaction part of Eq. (1) is obtained from the spatial extension of Eq. (7).

One of the important conclusions of our previous work⁶¹ was that the Langevin description (7) is not consistent with equilibrium statistical mechanics. Alternative formulations based on a Langevin diffusion description^86,87 that are consistent at thermodynamic equilibrium fail to correctly model relaxation toward equilibrium.⁶¹ Instead, in order to correctly capture both small fluctuations and large deviations in equilibrium and non-equilibrium contexts, one must retain a description of chemical reactions as a Markov jump process. That is, one must describe reactions using a stochastic differential equation driven by Poisson rather than Gaussian noise.

B. Structure factor

The structure factor is the steady-state spectrum of the concentration fluctuations,

S_{s} (𝒌) = V ⟨ δ {\hat{n}}_{s, 𝒌} δ {\hat{n}}_{s, 𝒌}^{*} ⟩,

(8)

i.e., the variance of the Fourier mode of the number density of species s,

{\hat{n}}_{s, 𝒌} (t) = \frac{1}{V} \int n_{s} (𝒙, t) e^{- i 𝒌 \cdot 𝒙} d 𝒙 .

(9)

Here we have assumed a periodic domain of volume V and defined $δ {\hat{n}}_{s, 𝒌} = {\hat{n}}_{s, 𝒌} - ⟨ {\hat{n}}_{s, 𝒌} ⟩$ ⁠, where the brackets $⟨ ⟩$ denote the equilibrium average. Here we derive an analytic expression of the structure factor from the linearized FHD equation. We assume that there is only one species, N_s = 1, and drop the subscript s for species, to write Eq. (1) as

\frac{\partial}{\partial t} n = D \nabla^{2} n + \nabla \cdot (\sqrt{2 D n} ^{(D)}) + a (n) + \sqrt{2 Γ (n)} ^{(R)},

(10)

where

a (n) = \sum_{r = 1}^{N_{r}} ν_{r} a_{r} (n), Γ (n) = \frac{1}{2} \sum_{r = 1}^{N_{r}} ν_{r}^{2} a_{r} (n),

(11)

and we have expressed fluctuations arising from all reactions by a single GWN field $^{(R)}$ ⁠. At a spatially uniform stable steady state, $n (𝒙, t)$ fluctuates around mean number density $\bar{n} \equiv ⟨ n ⟩$ ⁠, where $a (\bar{n}) = 0$ and $a^{'} (\bar{n}) < 0$ ⁠. The linearization of Eq. (10) around this equilibrium state is given by the central limit theorem,

\frac{\partial}{\partial t} n = D \nabla^{2} n + \sqrt{2 D \bar{n}} \nabla \cdot ^{(D)} - r (n - \bar{n}) + \sqrt{2 \bar{Γ}} ^{(R)},

(12)

where $r = - a^{'} (\bar{n}) > 0$ is the effective reaction rate near equilibrium and $\bar{Γ} = Γ (\bar{n})$ ⁠.

The Fourier transform of Eq. (12) gives

\frac{d}{d t} δ {\hat{n}}_{𝒌} = - D k^{2} δ {\hat{n}}_{𝒌} + \sqrt{2 D \bar{n}} i 𝒌 \cdot {\hat{}}_{𝒌}^{(D)} - r δ {\hat{n}}_{𝒌} + \sqrt{2 \bar{Γ}} {\hat{}}_{𝒌}^{(R)} .

(13)

Since Eq. (13) has the form of the Ornstein–Uhlenbeck equation,³² the structure factor is easily obtained as

S (𝒌) = \frac{D \bar{n} k^{2} + \bar{Γ}}{D k^{2} + r} = \frac{\bar{n} k^{2} + \bar{Γ} / D}{k^{2} + ℓ^{- 2}},

(14)

where $ℓ = \sqrt{D / r}$ denotes the penetration depth. From Eq. (14), we observe that there are two limiting cases. In the small wave number limit $k ℓ ≪ 1$ ⁠, $S (𝒌)$ becomes $\bar{Γ} / r$ and does not depend on diffusion. In fact, the result $S (𝟎) = \bar{Γ} / r$ is also obtained from the CME assuming the whole system is well-mixed. On the other hand, in the large wave number limit $k ℓ ≫ 1$ ⁠, $S (𝒌)$ becomes $\bar{n}$ ⁠, which is the result for the diffusion-only system. Hence, fluctuations are reaction-dominated at a length scale larger than $ℓ$ and are diffusion-dominated at a length scale smaller than $ℓ$ ⁠.

We also observe that if the system is in detailed balance at its steady state, i.e., it is in thermodynamic equilibrium, then $\bar{Γ} = \bar{n} r$ and $S (𝒌) = \bar{n}$ ⁠, consistent with a product Poisson distribution with mean number density $\bar{n}$ ⁠. Therefore, in true thermodynamic equilibrium, the statistics of the fluctuations are independent of any kinetic parameters, as they must be according to equilibrium statistical mechanics.⁸⁸ In particular, the presence of the reactions does not change the Poisson statistics of the state of thermodynamic equilibrium. In Section V A, we use this property to judge the quality of numerical schemes.

C. Schlögl model

The Schlögl model^75,76 is given by the chemical reactions for species X,

2 X \underset{k_{2}}{\overset{k_{1}}{⇌}} 3 X, \emptyset \underset{k_{4}}{\overset{k_{3}}{⇌}} X .

(15)

Hence, we have N_s=1, N_r=4, $ν_{1} = ν_{3} = 1$ ⁠, $ν_{2} = ν_{4} = - 1$ ⁠, a(n)=k₁n²−k₂n³+k₃−k₄n, and $Γ (n) = \frac{1}{2} (k_{1} n^{2} + k_{2} n^{3} + k_{3} + k_{4} n)$ ⁠. Due to the cubic nonlinearity of a(n), the well-mixed system exhibits several kinds of distributions depending on the values of the rate constants. If detailed balance is satisfied, that is, $k_{1} n_{eq}^{2} = k_{2} n_{eq}^{3}$ and k₃= k₄n_eq, the system is in thermodynamic equilibrium and the distribution follows Poisson statistics with mean number density n_eq. Otherwise, depending on the number of real roots of a(n)= 0, the system exhibits a monostable distribution (for a single positive root) or a bistable distribution (for three positive roots).⁷⁶

The structure factor of the spatially extended Schlögl model can be calculated from Eq. (14). As expected from the fact that the equilibrium distribution of the system follows Poisson statistics, $S (𝒌) = n_{eq}$ in the case of thermodynamic equilibrium. Note, however, that having a monostable distribution does not imply thermodynamic equilibrium. The structure factor of the out-of-equilibrium monostable case is not flat but exhibits a transition near $k ℓ \sim 1$ ⁠. For the bistable case exhibiting metastability, the linearized theory is still applicable if one looks at fluctuations around one of the two peaks. However, in this work, we focus on the equilibrium and out-of-equilibrium cases where a(n) has a single positive root.

III. SPATIAL DISCRETIZATION

In this section, we discuss spatial discretization of the FHD equation using a finite-volume approach^63,64 that converts the SPDE into stochastic ordinary differential equations (SODEs) for the cell number density $n_{s, 𝒊} (t)$ ⁠. We develop numerical schemes to solve these SODEs in Section IV. In Section III A, we first discretize the diffusion-only SPDE (3). In Section III B, we add reactions and present the continuous-time discrete-space description of the reaction-diffusion system. In Section III C, we discuss techniques to handle a small number of molecules per cell.

For simplicity, in this paper we only consider periodic systems. However, our methods can be straightforwardly generalized to standard types of physical boundary conditions (Dirichlet, Neumann, or Robin). In particular, since chemistry is local and does not require boundary conditions, one can rely on methods we have developed in prior work without chemistry; see, for example, the discussion in Ref. 89.

A. Diffusion-only case

Due to the lack of regularity of $(_{s}^{(D)} (𝒙, t)$ in Eq. (3), pointwise values of $n_{s} (𝒙, t)$ are not physically meaningful. Hence, we consider instead the spatial average of $n_{s} (𝒙, t)$ over a cell. We partition the system domain $L_{1} \times L_{2} \times \dots \times L_{d}$ into cells of volume $Δ V = Δ x_{1} \dots Δ x_{d}$ and denote the cell number density of species s in cell $𝒊 = (i_{1}, \dots, i_{d})$ as

n_{s, 𝒊} (t) = \frac{1}{Δ V} \int_{cell 𝒊} n_{s} (𝒙, t) d 𝒙 .

(16)

We denote the face of a cell using the index f. If two contiguous cells have indices i and $𝒊 + 𝒆_{j}$ (with $𝒆_{j}$ being the unit vector along the j-axis), the face f shared by the cells is denoted by $𝒊 + \frac{1}{2} 𝒆_{j}$ ⁠.

To obtain a spatial discretization of Eq. (3) that ensures discrete fluctuation-dissipation balance,^63,64 we use the standard second-order discrete Laplacian operator for the deterministic diffusion part $D {}_{s}\nabla^{2} n_{s}$ and introduce a staggered grid for the stochastic diffusive flux term $\nabla \cdot (\sqrt{2 D_{s} n_{s}} (_{s}^{(D)})$ ⁠, see Fig. 1. For d = 1, a formally second-order spatial discretization of Eq. (3) is written as

\frac{d}{d t} n_{s, i} = D_{s} \frac{n_{s, i + 1} - 2 n_{s, i} + n_{s, i - 1}}{Δ x^{2}} + \sqrt{\frac{2 D_{s}}{Δ V}} \frac{\sqrt{{\tilde{n}}_{s, i + \frac{1}{2}}} W_{s, i + \frac{1}{2}} - \sqrt{{\tilde{n}}_{s, i - \frac{1}{2}}} W_{s, i - \frac{1}{2}}}{Δ x} .

(17)

The spatial average of $n_{s} (𝒙, t)$ over the interval of length $Δ x$ around face $i \pm \frac{1}{2}$ is approximated by ${\tilde{n}}_{s, i \pm \frac{1}{2}} (t)$ ⁠, whereas that of $_{s} (𝒙, t)$ is modeled by $\frac{1}{\sqrt{Δ V}} W_{s, i \pm \frac{1}{2}} (t)$ ⁠. To close the equation, ${\tilde{n}}_{s, i \pm \frac{1}{2}} (t)$ is approximated by an average of n_s,i(t) and $n_{s, i \pm 1} (t)$ ⁠, that is, ${\tilde{n}}_{s, i \pm \frac{1}{2}} = \tilde{n} (n_{s, i}, n_{s, i \pm 1})$ ⁠. Natural candidates for the averaging function $\tilde{n} (n_{1}, n_{2})$ would be the Pythagorean means: the arithmetic, geometric, and harmonic means. We choose a modified arithmetic average for $\tilde{n} (n_{1}, n_{2})$ described in Section III C, for reasons detailed in Appendix B.

FIG. 1.

View large Download slide

Finite-volume spatial discretization in two dimensions. The cell-averaged number density $n_{s, 𝒊} (t)$ is associated with the circles, and the face-averaged stochastic diffusive flux is associated with the crosses. The stochastic diffusive flux between the two cells having the red circles at the center is depicted by the blue arrow.

Generalization of the spatial discretization (17) to higher dimensions is straightforward. For each face, a GWN process $W_{s, 𝒇}$ is defined and ${\tilde{n}}_{s, 𝒇} (t)$ is calculated from the cell number densities of the two cells sharing the face by using the averaging function $\tilde{n} (n_{1}, n_{2})$ ⁠, see Fig. 1. By introducing notations $n_{s} (t) \equiv {n_{s, 𝒊} (t)}$ ⁠, $W_{s} (t) \equiv {W_{s, 𝒇} (t)}$ ⁠, and ${\tilde{n}}_{s} (t) \equiv \tilde{n} (n_{s} (t))$ ⁠, we express the resulting SODEs for ${n_{s, 𝒊} (t)}$ as

\frac{d}{d t} n_{s} = D {}_{s}\nabla_{d}^{2} n_{s} + \sqrt{\frac{2 D_{s}}{Δ V}} \nabla_{d} \cdot (\sqrt{{\tilde{n}}_{s}} W_{s}),

(18)

with the understanding that $\nabla_{d}^{2}$ denotes the standard (2d + 1)-point discrete Laplacian operator, and $\nabla_{d} \cdot$ denotes a discrete divergence operator.

B. Reaction-diffusion system

By combining Eqs. (6) and (18), we obtain the spatial discretization of the reaction-diffusion FHD equations as a system of Ito SODEs,

\frac{d}{d t} n_{s} = D {}_{s}\nabla_{d}^{2} n_{s} + \sqrt{\frac{2 D_{s}}{Δ V}} \nabla_{d} \cdot (\sqrt{{\tilde{n}}_{s}} W_{s}) + \sum_{r = 1}^{N_{r}} \frac{ν_{s r} P (a_{r} (𝒏) Δ V d t)}{Δ V d t} .

(19)

In Eq. (1), fluctuations in the reaction rate are modeled as GWN, while in Eq. (19), we assume Poisson fluctuations. Since the latter fluctuations are consistent with discrete nature of reactions based on the CME, the description in Eq. (19) is physically more accurate. In fact, it has been shown that the CLE description can give physically incorrect results since it is not consistent with a Gibbs–Boltzmann or Einstein equilibrium distribution, even for the case of a single well-mixed cell.⁶¹ As shown above, the inclusion of Poisson fluctuations for reaction, however, requires the notion of a mesoscopic cell and thus can be realized only after the SPDE is spatially discretized.

The choice of appropriate cell size is a delicate issue for the RDME and FHD descriptions. An upper bound on the cell size is given by the penetration depth due to the underlying assumption that each cell is homogeneous and reactions occur within a cell. In fact, there is not only an upper bound of the cell size for a valid description but also a lower bound. This can be seen by considering the fact that bimolecular reactions would become increasingly infrequent as the cell size decreases.^90,91 Several criteria for choosing the cell size have been proposed based on physical arguments^90,92,93 and mathematical analysis.⁹⁴ For a small value of the cell size, corrections in the rate constants of bimolecular reactions have been proposed.^90,95,96 However, these corrections do not fix the underlying problem which comes from the fact that reactions are treated as a purely local process with no associated spatial length scale. In microscopic (particle) models of reaction-diffusion such as the Smoluchowski¹⁷ model or the Doi model,⁹⁰ a microscopic reactive distance appears and controls the reaction rate for diffusion-limited reactions. By introducing a reactive distance into the model, and relaxing the restriction that a reaction should occur among the molecules in the same cell, a modified convergent RDME having well-defined limiting behavior for small cell size can be developed⁹⁷ and could be combined with our FHD description of diffusion. The dependence of stochastic Turing patterns on the grid size has been also investigated.⁹⁸

C. Maintaining non-negative densities

The spatially discretized FHD equation (18) or (19) is well defined but suffers from two issues that we now address. First, the number of molecules in a cell (i.e., $n_{s, 𝒊} Δ V$ ⁠) is not an integer. Second, the cell number density can become negative. When there are a small number of molecules per cell in the system, the behavior of the FHD description (19) depends sensitively on the averaging function $\tilde{n}$ and the propensity functions a_r that appear in the multiplicative noise terms. Hence, we carefully modify the form of $\tilde{n}$ and a_r for negative or very small densities in order to greatly reduce the chances of producing future negative densities.

In Section V A, we demonstrate that the arithmetic mean produces more accurate results for the equilibrium distribution than the other Pythagorean means. Based on the analysis given in Appendix B, we use the following modification to the arithmetic mean:

\tilde{n} (n_{1}, n_{2}) = \frac{n_{1} + n_{2}}{2} H_{0} (n_{1} Δ V) H_{0} (n_{2} Δ V),

(20)

where

H_{0} (x) = {\begin{matrix} 0 & (x \leq 0) \\ x & (0 < x < 1) \\ 1 & (x \geq 1) \end{matrix}

(21)

is a smoothed Heaviside function. The smoothed Heaviside function H₀ is introduced to ensure the continuity of $\tilde{n}$ at n₁ = 0 or n₂ = 0. As explained in Appendix B, this averaging function guarantees non-negativity for the diffusion-only system (18) in the continuous-time description. In our simulations, we find this modification greatly reduces the occurrence of negative density while closely matching the true equilibrium distribution, noting that in our formulation the stochastic diffusive flux is continuously turned off at $n_{1} \leq 0$ or $n_{2} \leq 0$ ⁠. We also note that the smoothing is based on the number of molecules in a cell and if both cells have at least one molecule (i.e., $n_{i} Δ V \geq 1$ ⁠), $\tilde{n}$ becomes exactly the arithmetic mean. As shown in Appendix B, the local modification near n = 0 does not cause any noticeable unphysical behavior for $n_{i} Δ V \geq 1$ ⁠. In Section V, we demonstrate that our numerical schemes based on Eq. (20) work very well even for a small number of molecules per cell.

For the propensity functions $a_{r} (𝒏)$ ⁠, we use the following correction to the law of mass action, which is usually included in the RDME description: if the deterministic rate expression contains $n_{s}^{2}$ (or $n_{s}^{3}$ ⁠, $\dots$ ⁠), then replace it by $n_{s} (n_{s} - \frac{1}{Δ V})$ (or $n_{s} (n_{s} - \frac{1}{Δ V}) (n_{s} - \frac{2}{Δ V})$ ⁠, $\dots$ ⁠). With this correction, at thermodynamic equilibrium, the mean reaction rate becomes equal to the one calculated from the deterministic rate expression with the mean number density. This can be seen from the fact that if $n_{s} Δ V$ follows Poisson statistics with mean ${\bar{n}}_{s} Δ V$ ⁠, $⟨ n_{s} (n_{s} - \frac{1}{Δ V}) ⟩ = {\bar{n}}_{s}^{2}$ and $⟨ n_{s} (n_{s} - \frac{1}{Δ V}) (n_{s} - \frac{2}{Δ V}) ⟩ = {\bar{n}}_{s}^{3}$ ⁠.

When reactions are combined with an FHD treatment of diffusion, number densities are no longer restricted to non-negative integers and special treatment is required to make reaction rates non-negative and physically sensible for small numbers of molecules. In this work, we evaluate the rate $a_{r} (𝒏)$ by using continuous-range number densities n (i.e., without trying to round $𝒏 Δ V$ to integers) and ensure that each term in the rate of each reaction is non-negative. For example, we take the rate expression of the Schlögl model (see Section II C) to be

a (n) = k_{1} n^{+} {(n - \frac{1}{Δ V})}^{+} - k_{2} n^{+} {(n - \frac{1}{Δ V})}^{+} {(n - \frac{2}{Δ V})}^{+} + k_{3} - k_{4} n^{+},

(22)

where $n^{+} = \max (n, 0)$ ⁠. We note that more mathematically justified algorithms have been proposed to handle reactions in regard to negative densities using operator splitting and exact solutions of reaction subproblems;^99,100 these methods cannot be used to address negative densities due to stochastic diffusive fluxes.

IV. TEMPORAL INTEGRATORS

In this section, we develop temporal integrators for the spatially discretized FHD equation (19). Our goal is twofold. First, we construct numerical methods that allow for a large time step size even in the presence of fast diffusion. By treating diffusion implicitly, the severe restriction on time step size can be bypassed. Second, we construct methods that maintain accuracy even if the time step size is much larger than the diffusive hopping time. Since it is quite difficult to achieve second-order weak accuracy for general multiplicative noise,¹⁰¹ our goal here is to ensure second-order accuracy wherever possible. In the limit in which the number of molecules per cell is very large and one can replace random numbers by their means, our schemes reduce to standard second-order schemes for deterministic reaction-diffusion PDEs. For linearized FHD, our midpoint tau leaping-based schemes are second-order weakly accurate, and all midpoint schemes reproduce at least second-order accurate static correlations, i.e., structure factors.

We build on previous work by some of us in Refs. 63, 64, and 74 and propose two (semi-) implicit schemes as an alternative numerical method to conventional RDME methods. We mainly consider the case where diffusion is much faster than reaction and molecules on average diffuse more than a cell length per time step (i.e., $2 d D Δ t ≫ Δ x^{2}$ ⁠). We focus here on unsplit schemes that do not rely on operator splitting. This is because we found that unsplit schemes give notably more accurate structure factors than corresponding split schemes in our case. In addition, including other transport processes (e.g., advection) and handling boundary conditions¹⁰² to second order is not straightforward for split schemes.

It is convenient to introduce dimensionless numbers, $α$ and $β$ ⁠, which measure how fast diffusion and reaction are relative to the given time step size $Δ t$ ⁠, respectively. For the single-species Equation (12), assuming $Δ x_{1} = \dots = Δ x_{d} = Δ x$ ⁠, we define

α = r Δ t, β = \frac{D Δ t}{Δ x^{2}},

(23)

where r is the chemical relaxation rate appearing in the linearized equation (12). Hence, we can express the well-mixed condition (i.e., the penetration depth $ℓ = \sqrt{D / r} ≫ Δ x$ ⁠) as $α ≪ β$ ⁠. In addition, the numerical stability condition of a scheme can also be given in terms of $α$ and $β$ ⁠. That is, if reaction and/or diffusion are treated explicitly in a scheme, values of $α$ and/or $β$ larger than a stability threshold cause numerical instability. For $α ≪ β$ ⁠, the stability limit is mainly determined by fast diffusion,

β \leq \frac{1}{2 d} \Leftrightarrow Δ t \leq Δ t_{\max} \equiv \frac{Δ x^{2}}{2 d D} .

(24)

Note that the stability limit becomes severe for large diffusion coefficients and small grid spacing and worsens with increasing dimension.

In Section IV A, we present several numerical schemes for the FHD equation (19), including two implicit schemes, and analyze the temporal orders of accuracy for the structure factors using the linearized analysis described in Appendix C. In Section IV B, we analyze the stochastic accuracy of the numerical schemes for large $Δ t$ by investigating the structure factor of the one-dimensional Schlögl model at different wavenumbers. Since analysis for the nonlinear equations is lacking at present, we numerically justify the handling of multiplicative noise in Section V A.

A. Schemes

The simplest method for integrating Eq. (19) in time is the Euler–Maruyama tau leaping (EMTau) scheme,

n_{s}^{k + 1} = n_{s}^{k} + D_{s} Δ t \nabla_{d}^{2} n_{s}^{k} + \sqrt{\frac{2 D_{s} Δ t}{Δ V}} \nabla_{d} \cdot (\sqrt{{\tilde{n}}_{s}^{k}} W_{s}^{k}) + \sum_{r = 1}^{N_{r}} \frac{ν_{s r} P (a_{r}^{k} Δ V Δ t)}{Δ V},

(25)

where superscripts denote the point in time at which quantities are evaluated, e.g., $n_{s}^{k} = n_{s} (k Δ t)$ and $a_{r}^{k} = a_{r} (𝒏 (k Δ t))$ ⁠, and we have used the compact notation for spatial discretization introduced in Section III B. Here, $\int_{k Δ t}^{(k + 1) Δ t} W (t^{'}) d t^{'}$ has been replaced by $\sqrt{Δ t} W^{k}$ ⁠, where W^k denotes a collection of standard random Gaussian variables sampled independently for each species on each grid face at each time step. That is, the stochastic diffusive flux of species s on face f at time step k is proportional to $W_{s, 𝒇}^{k}$ ⁠.

We also construct numerical schemes where reactions are treated by SSA, which is an exact (exponential) integrator for reactions. We denote by $ℜ_{s} (𝒏, τ)$ the (random) change in the number density of species s for a cell with initial state n obtained from SSA over at time interval $τ$ (in the absence of diffusion). We can then write the Euler–Maruyama SSA (EM-SSA) scheme as

n_{s}^{k + 1} = n_{s}^{k} + D_{s} Δ t \nabla_{d}^{2} n_{s}^{k} + \sqrt{\frac{2 D_{s} Δ t}{Δ V}} \nabla_{d} \cdot (\sqrt{{\tilde{n}}_{s}^{k}} W_{s}^{k}) + ℜ_{s} (𝒏^{k}, Δ t) .

(26)

The EMTau scheme is explicit in the sense that all terms on the right-hand side of Eq. (25) can be evaluated without knowing $n_{s}^{k + 1}$ ⁠. However, a by a stability condition (for derivation, see Eq. (C7))

β + \frac{α}{4 d} \leq \frac{1}{2 d},

(27)

which reduces to condition (24) for $α ≪ β$ ⁠. Since the EM-SSA scheme treats reactions using an exponential integrator, it is only subject to the stability limit (24) without a restriction on $α$ ⁠.

The stability limit imposed by fast diffusion can be overcome by using standard implicit methods such as the second-order implicit midpoint or Crank–Nicolson method, which gives a system of linear equations for $n_{s}^{k + 1}$ ⁠,

n_{s}^{k + 1} = n_{s}^{k} + D_{s} Δ t \nabla_{d}^{2} (\frac{n_{s}^{k} + n_{s}^{k + 1}}{2}) + \sqrt{\frac{2 D_{s} Δ t}{Δ V}} \nabla_{d} \cdot (\sqrt{{\tilde{n}}_{s}^{k}} W_{s}^{k}) + \sum_{r = 1}^{N_{r}} \frac{ν_{s r} P (a_{r}^{k} Δ V Δ t)}{Δ V} .

(28)

The linear system (28) can be solved efficiently iteratively using multigrid relaxation;¹⁰³ for $β ≲ 1$ ⁠, solving the linear system is not much more expensive than a step of a second-order explicit time stepping scheme. Note, however, that scheme (28) is only first order accurate overall for reaction-diffusion systems. Hence, in addition to improved stability, it is important to develop higher-order schemes to improve accuracy. Note that this is not as simple as replacing tau leaping in Eq. (28) with SSA as in Eq. (26); this would still be only first-order accurate even in the deterministic limit.

Here we construct numerical schemes based on the second-order temporal integrators for the linearized equations of FHD developed in Refs. 64 and 74. Those temporal integrators are second-order accurate in the weak sense for additive noise and are used here as the basis for handling diffusion. In order to add reactions into diffusion-only schemes, we consider two types of sampling methods: tau leaping and SSA. For tau leaping, we use the weakly second-order tau leaping method;^65,66 a similar two-stage scheme has been originally proposed for the CLE¹⁰⁴ to achieve second-order weak accuracy. Here we combine predictor-corrector midpoint schemes proposed in Ref. 64 (for diffusion) and the second-order tau leaping method (for reaction). Owing to the similar two-stage structures of those schemes, they fit together in a rather natural manner. In addition, since the resulting schemes still fit the framework of the implicit-explicit algorithms analyzed in Ref. 64, they are second-order weakly accurate for additive noise and an additional order of accuracy (i.e., third order) is gained for the structure factor.

We also develop midpoint schemes that use SSA instead of tau leaping for reactions. Unlike tau leaping-based schemes, the SSA-based schemes do not suffer from instability even in the presence of rapid reactions. The use of SSA may also help prevent the development of negative densities, which is one of the main numerical issues for large fluctuations. Hence, while SSA-based numerical schemes are computationally more expensive, they work better than tau leaping-based schemes when reactions are fast or when the number of molecules is small. The SSA-based schemes we propose here belong to a class of exponential Runge–Kutta schemes, and we construct them to ensure second-order deterministic accuracy, as well as second-order accuracy for the structure factor; a detailed analysis of their weak accuracy is at present missing even for linearized FHD.

1. Explicit midpoint schemes

As a prelude to constructing two-stage implicit methods, we first consider improving the accuracy of the explicit EMTau scheme (25) by using an explicit two-stage Runge–Kutta (predictor-corrector) approach. By combining the explicit midpoint predictor-corrector scheme from Refs. 64 and 74 (for diffusion) and the midpoint tau leaping scheme from Refs. 65 and 66 (for reaction), we obtain the explicit midpoint tau leaping (ExMidTau) scheme,

n_{s}^{⋆} = n_{s}^{k} + \frac{D_{s} Δ t}{2} \nabla_{d}^{2} n_{s}^{k} + \sqrt{\frac{D_{s} Δ t}{Δ V}} \nabla_{d} \cdot (\sqrt{{\tilde{n}}_{s}^{k}} W_{s}^{(1)}) + \sum_{r = 1}^{N_{r}} \frac{ν {}_{s r}P^{(1)} (a_{r}^{k} Δ V Δ t / 2)}{Δ V},

(29a)

n_{s}^{k + 1} = n_{s}^{k} + D_{s} Δ t \nabla_{d}^{2} n_{s}^{⋆} + \sqrt{\frac{D_{s} Δ t}{Δ V}} \nabla_{d} \cdot (\sqrt{{\tilde{n}}_{s}^{k}} W_{s}^{(1)}) + \sqrt{\frac{D_{s} Δ t}{Δ V}} \nabla_{d} \cdot (\sqrt{{\tilde{n}}_{s}^{∙}} W_{s}^{(2)}) + \sum_{r = 1}^{N_{r}} \frac{ν {}_{s r}P^{(1)} (a_{r}^{k} Δ V Δ t / 2)}{Δ V} + \sum_{r = 1}^{N_{r}} \frac{ν {}_{s r}P^{(2)} ({(2 a_{r}^{⋆} - a_{r}^{k})}^{+} Δ V Δ t / 2)}{Δ V},

(29b)

where the superscripts (1) and (2) indicate that the terms correspond to the first and second half of the time step, respectively. That is, $P^{(1)}$ (and similarly for W⁽¹⁾ and other random increments) denotes the same random number in both predictor and corrector stages and is only sampled once per time step. Following Refs. 65 and 66, the mean reaction rate for the second half step is corrected to ${(2 a_{r}^{⋆} - a_{r}^{k})}^{+}$ ⁠, where $a^{+} = \max (a, 0)$ ⁠.

For the magnitude of the stochastic diffusive fluxes over the second half of the time step, we consider the following three options for the face average value ${\tilde{n}}_{s}^{∙}$ ⁠:

{\tilde{n}}_{s}^{∙} = \tilde{n} (n_{s}^{k}),

(30a)

{\tilde{n}}_{s}^{∙} = \tilde{n} (n_{s}^{⋆}),

(30b)

{\tilde{n}}_{s}^{∙} = \tilde{n} ({(2 n_{s}^{⋆} - n_{s}^{k})}^{+}) . `

(30c)

While all options are consistent with the Ito interpretation, the effect of this choice on accuracy requires a nonlinear analysis that is not available at present. The option (30b) is used in Ref. 74 and shown to lead to second-order weak accuracy for FHD equations linearized around a time-dependent macroscopic state. The option (30c) is inspired by the midpoint tau leaping scheme.^65,66 However, it does not actually lead to second-order weak accuracy for multiplicative noise because the fluctuating diffusion equation does not have the simple noise structure that the CLE has.¹⁰⁴ For all our simulations, we use option (30c), as justified by numerical results in Section V A.

The reactions can also be treated using SSA, to give the explicit midpoint SSA (ExMidSSA) scheme,

n_{s}^{⋄} = n_{s}^{k} + \frac{D_{s} Δ t}{2} \nabla_{d}^{2} n_{s}^{k} + \sqrt{\frac{D_{s} Δ t}{Δ V}} \nabla_{d} \cdot (\sqrt{{\tilde{n}}_{s}^{k}} W_{s}^{(1)}),

(31a)

n_{s}^{⋆} = n_{s}^{⋄} + ℜ_{s}^{(1)} (𝒏^{⋄}, \frac{Δ t}{2}),

(31b)

n_{s}^{k + 1} = n_{s}^{k} + D_{s} Δ t \nabla_{d}^{2} n_{s}^{⋆} + \sqrt{\frac{D_{s} Δ t}{Δ V}} \nabla_{d} \cdot (\sqrt{{\tilde{n}}_{s}^{k}} W_{s}^{(1)}) + \sqrt{\frac{D_{s} Δ t}{Δ V}} \nabla_{d} \cdot (\sqrt{{\tilde{n}}_{s}^{∙}} W_{s}^{(2)}) + ℜ_{s}^{(1)} (𝒏^{⋄}, \frac{Δ t}{2}) + ℜ_{s}^{(2)} (𝒏^{⋆}, \frac{Δ t}{2}) .

(31c)

Here the predictor stage (31a) + (31b) is a split reaction-diffusion step, but the corrector is not split. Note that two $ℜ^{(1)}$ appearing in Eqs. (31b) and (31c) are the same random increment computed using SSA. In other words, the SSA algorithm is called once for each half of the time step; this has the same computational cost as calling SSA once to compute $ℜ_{s} (𝒏^{k}, Δ t)$ in the EM-SSA scheme (26).

Since both ExMidTau and ExMidSSA schemes treat diffusion explicitly, they are subject to stability limits. The ExMidTau scheme has the same stability limit (27) as the EMTau scheme, whereas the ExMidSSA scheme is subject to the same limit (24) as the EM-SSA scheme. The ExMidTau scheme with the option (30b) is an instance of the explicit midpoint scheme analyzed in Ref. 74 for weak noise (i.e., linearized FHD) and therefore achieves second-order weak accuracy for linearized FHD and gives third-order accurate equilibrium structure factors. On the other hand, the ExMidSSA scheme gives only second-order accurate structure factors.

2. Implicit midpoint schemes

Here we present two implicit midpoint schemes, where diffusion is treated implicitly based on the implicit midpoint predictor-corrector scheme.^64,74 By treating reactions using the second-order midpoint tau leaping scheme,^65,66 we obtain the implicit midpoint tau leaping (ImMidTau) scheme,

n_{s}^{⋆} = n_{s}^{k} + \frac{D_{s} Δ t}{2} \nabla_{d}^{2} n_{s}^{⋆} + \sqrt{\frac{D_{s} Δ t}{Δ V}} \nabla_{d} \cdot (\sqrt{{\tilde{n}}_{s}^{k}} W_{s}^{(1)}) + \sum_{r = 1}^{N_{r}} \frac{ν {}_{s r}P^{(1)} (a_{r}^{k} Δ V Δ t / 2)}{Δ V},

(32a)

n_{s}^{k + 1} = n_{s}^{k} + D_{s} Δ t \nabla_{d}^{2} (\frac{n_{s}^{k} + n_{s}^{k + 1}}{2}) + \sqrt{\frac{D_{s} Δ t}{Δ V}} \nabla_{d} \cdot (\sqrt{{\tilde{n}}_{s}^{k}} W_{s}^{(1)}) + \sqrt{\frac{D_{s} Δ t}{Δ V}} \nabla_{d} \cdot (\sqrt{{\tilde{n}}_{s}^{∙}} W_{s}^{(2)}) + \sum_{r = 1}^{N_{r}} \frac{ν {}_{s r}P^{(1)} (a_{r}^{k} Δ V Δ t / 2)}{Δ V} + \sum_{r = 1}^{N_{r}} \frac{ν {}_{s r}P^{(2)} ({(2 a_{r}^{⋆} - a_{r}^{k})}^{+} Δ V Δ t / 2)}{Δ V},

(32b)

where the three options for ${\tilde{n}}_{s}^{∙}$ are given in Eqs. (30). When SSA is used for the reactions, we obtain the implicit midpoint SSA (ImMidSSA) scheme,

n_{s}^{⋆} = n_{s}^{k} + \frac{D_{s} Δ t}{2} \nabla_{d}^{2} n_{s}^{⋆} + \sqrt{\frac{D_{s} Δ t}{Δ V}} \nabla_{d} \cdot (\sqrt{{\tilde{n}}_{s}^{k}} W_{s}^{(1)}),

(33a)

n_{s}^{k + 1} = n_{s}^{k} + D_{s} Δ t \nabla_{d}^{2} (\frac{n_{s}^{k} + n_{s}^{k + 1}}{2}) + \sqrt{\frac{D_{s} Δ t}{Δ V}} \nabla_{d} \cdot (\sqrt{{\tilde{n}}_{s}^{k}} W_{s}^{(1)}) + \sqrt{\frac{D_{s} Δ t}{Δ V}} \nabla_{d} \cdot (\sqrt{{\tilde{n}}_{s}^{∙}} W_{s}^{(2)}) + ℜ_{s} (𝒏^{⋆}, Δ t) .

(33b)

In the corrector stage, both schemes treat diffusion using the Crank–Nicolson method since this gives the most accurate structure factors for diffusion-only systems.⁶⁴ For the predictor step to the midpoint, we have chosen to use backward Euler for diffusion because this was found to be optimal using the structure factor analysis discussed in more detail in Section IV B.

We point out again that the difference in how reactions are included in the ImMidTau and ImMidSSA schemes stems from the fact that SSA is an exponential integrator whereas the midpoint tau leaping method is only a second-order integrator. This difference must be taken into account when analyzing the accuracy of SSA-based schemes both in the deterministic limit and in structure factor analysis. Like the explicit midpoint schemes in Section IV A 1, for linearized FHD, the ImMidTau scheme is second-order weakly accurate and gives a third-order accurate structure factor, whereas the ImMidSSA scheme gives only a second-order accurate structure factor. Since diffusion is treated implicitly in both schemes, they are not subject to a stability limit depending on $β$ ⁠. However, due to the explicit treatment of reactions, the ImMidTau scheme is subject to the stability condition $α \leq 2$ ⁠. The ImMidSSA scheme is unconditionally linearly stable and has no stability restrictions on $α$ and $β$ but can be considerably more expensive for systems with a large number of molecules.

A key element of this work that distinguishes it from our previous work based on the CLE⁶¹ is that here we replaced the GWN in the chemical noise with Poisson noise and used a weakly second-order tau leaping method^65,66 to account for the non-Gaussian nature of the chemical fluctuations. It is important to note that Poisson noise does not have a continuous range limit, i.e., the Poisson distribution remains integer-valued even as the number of molecules per cell becomes very large. Although it is tempting to replace the Poisson distribution with a Gaussian distribution, this changes the large deviation functional and therefore we recommend using tau leaping even in the case of weak fluctuations; we note that sampling from a Poisson distribution can be done with a cost essentially independent of the mean using well-designed rejection Monte Carlo methods. Because of the use of Poisson variates, which cannot be split into a mean and a fluctuation like a Gaussian variate can, there is no strict “deterministic limit” for our FHD discretizations. While the handling of diffusion degenerates to a standard second-order deterministic scheme in the absence of the noise, the chemical noise is always present and increments or decrements the number of molecules by integer numbers.

B. Structure factor analysis

Analyzing the accuracy of temporal integrators for stochastic differential equations is notably nontrivial, especially if driven by multiplicative noise. As mentioned above, because of the multiplicative noise, all of our midpoint schemes are formally only first-order weakly accurate. However, traditional weak-order accuracy is not the most important goal in FHD simulations. As first argued in Ref. 63 and then elaborated in Refs. 64 and 74, for FHD it is more important to attain discrete fluctuation-dissipation balance and higher-order accuracy for the spectrum of the equilibrium fluctuations.

Here, we analyze the accuracy of our numerical schemes by investigating the structure factor S(k) for the one-dimensional linearized FHD equation (12). The analytic expression for S(k) produced by a given scheme can be obtained as a function of $Δ x$ and $Δ t$ following the procedure described in Appendix C. Of specific interest to us is how accurately the implicit schemes reproduce S(k) at large wavenumbers corresponding to length scales comparable to $Δ x$ (i.e., $k ≳ {(D Δ t)}^{- 1 / 2}$ ⁠) when diffusion is the fastest process, $β ≫ 1 ≫ α$ ⁠. This is because incorrect diffusive dynamics at grid scales for $Δ t ≫ Δ x^{2} / D$ can lead to gross errors in the magnitude of the fluctuations at large wavenumbers.

Errors in the structure factor arise from two sources: spatial and temporal discretization. As explained in Appendix C, the predominant contribution of spatial discretization is to replace −k² with the symbol of the standard discrete Laplacian. In one dimension, this simply amounts to replacing k in the continuum expressions with the modified wavenumber $\tilde{k}$ defined by

\tilde{k} = \frac{\sin (\frac{k Δ x}{2})}{\frac{Δ x}{2}} .

(34)

Note that exactly the same expression applies to the RDME, where diffusion is simulated by lattice hops. In order to focus our attention on temporal integration errors, we will plot discrete structure factors as a function of $\tilde{k}$ instead of k, which effectively removes the spatial errors.

Figure 2 illustrates how S(k) deviates from the exact result (14) at different wavenumbers for large $Δ t$ ⁠. We compare S(k) obtained from the four schemes by using values $α = 0.025$ ⁠, $β = 10 α = 0.25$ for the explicit schemes and $α = 0.5$ ⁠, $β = 10 α = 5$ for the implicit schemes. Note that these values correspond to a case where the time step size is chosen as half of the stability limit $Δ t_{\max}$ for the explicit schemes, and it is increased by a factor of 20 (⁠ $Δ t = 10 Δ t_{\max}$ ⁠) for the implicit schemes. As described below, the accuracy at diffusion-dominated scales $k ℓ ≫ 1$ and reaction-dominated scales $k ℓ ≪ 1$ largely depends on how diffusion (i.e., explicit or implicit) and reaction (i.e., tau leaping or SSA) are treated, respectively.

FIG. 2.

View large Download slide

Discrete structure factors $S (\tilde{k})$ for the ExMidTau, ImMidTau, ExMidSSA, and ImMidSSA schemes for the one-dimensional linearized FHD equation (12) with $\bar{Γ} / \bar{n} r = 2$ (e.g., an out-of-equilibrium monostable Schlögl model). Note that different values of $β$ are chosen for the explicit schemes (⁠ $β = 0.25$ ⁠) and the implicit schemes (⁠ $β = 5$ ⁠) and $α$ is chosen as $α = 0.1 β$ The exact continuum result (14) is depicted by the dotted line.

As the time step size approaches the diffusive stability limit, $β \to 1 / 2$ in one dimension, the explicit schemes become inaccurate and eventually numerically unstable at the largest wavenumbers, as seen in the figure for $β = 1 / 4$ ⁠. Due to the small values of $α$ ⁠, both schemes give accurate results for reaction-dominated scales. Note, however, that the SSA-based schemes give exact S(0) and they are in general more accurate at reaction-dominated scales than the corresponding tau leaping-based schemes. Hence, we see that, for $β ≫ α$ ⁠, the accuracy of the explicit schemes is largely affected by numerical instability arising from diffusion.

On the other hand, both implicit schemes give fairly good S(k) in the overall range of k even though $Δ t$ is twenty times larger and $β = 5$ ⁠. As expected, the ImMidSSA scheme is more accurate for reaction-dominated scales $k ℓ ≪ 1$ ⁠, for which the ImMidTau scheme shows some errors because of the relatively large value of $α = 0.5$ ⁠. However, for intermediate scales $k ℓ \sim 1$ ⁠, the ImMidTau scheme is more accurate because it attains third-order accuracy for static covariances.

At diffusion-dominated scales $k ℓ ≫ 1$ ⁠, both implicit schemes give accurate results. This is not accidental, for we have selected these schemes from a family of schemes parameterized in Ref. 64 exactly for this reason. Specifically, the treatment of diffusion in both schemes is based on the implicit midpoint scheme (Crank–Nicolson), which gives the exact S(k) in the absence of reactions.⁶⁴ In addition, the reaction is incorporated in such a way that maintains fluctuation-dissipation balance for $k ℓ ≫ 1$ even for relatively large values of $α$ ⁠. For small $α$ ⁠, the time integration error of the ImMidTau scheme for the structure factor at the maximum wavenumber $k_{\max} ℓ = π \sqrt{β / α}$ is estimated as

\frac{S (k_{\max}) - S_{0} (k_{\max})}{S_{0} (k_{\max})} \approx - \frac{β}{2 {(1 + 2 β)}^{2}} α^{2},

(35)

where $S_{0} (k) = {lim}_{Δ t \to 0} S (k)$ is the structure factor in the absence of temporal integration errors (see Appendix C). Hence, for a given value of $α$ ⁠, S(k) gives accurate results at $k ℓ ≫ 1$ for large $β$ ⁠. For the ImMidSSA scheme, $β$ in the numerator is replaced by $β + (1 + 2 β) (\frac{\bar{Γ}}{\bar{n} r} - 1)$ and a similar stable behavior for large $β$ is observed. By expanding S(k_max) − S₀(k_max) for small $Δ t$ ⁠, we also see that the error is $O (α^{2} β) = O (Δ t^{3})$ for the ImMidTau scheme, whereas it is $O (Δ t^{2})$ for the ImMidSSA except at thermodynamic equilibrium (i.e., except when $\bar{Γ} = \bar{n} r$ ⁠), where it is third-order accurate.

V. NUMERICAL RESULTS

We perform numerical simulations for the following three stochastic reaction-diffusion systems. The simulation code is available at https://github.com/BoxLib-Codes/FHD_ReactDiff.git. In Section V A, we use the equilibrium Schlögl model in one, two, and three dimensions to validate our numerical methods. The analysis in Section IV B assumed additive noise, reflecting a large number of molecules per cell. Here we present numerical results demonstrating that the methodology continues to work when there are a small number of molecules per cell and the effects of multiplicative noise are significant. In particular, we show that our numerical methods, including the modified arithmetic-mean averaging function discussed in Section III C, accurately reproduce the Poisson statistics that characterize the thermodynamic equilibrium distribution. In Sections V B and V C, we study the effects of fluctuations on chemical pattern formation. In Section V B, we test our numerical methods on a time-dependent problem: two-dimensional Turing-like pattern formation in the three-species BPM model.^77,78 We investigate how accurately both time-transient and steady state behaviors are captured for the ImMidTau and ImMidSSA schemes when a large time step size is used. We consider the case where the populations of chemical species have different orders of magnitude, which is a frequently encountered situation where a conventional RDME-based method may not work efficiently. We demonstrate that the ImMidTau scheme scales very well with an increasing number of molecules per cell so that even the deterministic limit of vanishing fluctuations can be explored. In Section V C, we demonstrate the scalability to large systems and computational efficiency of our FHD approach by presenting a three-dimensional numerical simulation of chemical front propagation in a two-species model.⁷

As a reference method for comparison, we use an RDME-based method, as proposed in Appendix A, which is constructed via a standard operator splitting technique by combining multinomial diffusion sampling⁷¹ and SSA. Such a split scheme is notably more efficient than ISSA when there are a large number of molecules per cell and becomes an exact sampling method for the RDME in the limit $Δ t \to 0$ ⁠. This RDME-based scheme works with non-negative integer populations and reproduces correct fluctuations at thermodynamic equilibrium. However, diffusion imposes the same restriction (24) on $Δ t$ ⁠, and the split scheme produces only a first-order accurate structure factor in general.

In order to set a desired magnitude of fluctuations without changing any parameters for the macroscopic limit (e.g., penetration depth), we introduce a factor A, which scales the cell volume $Δ V = A Δ x_{1} \dots Δ x_{d}$ ⁠. It can be interpreted as the surface cross section in one dimension and the thickness of a system in two dimensions, and as a rescaling of the number density in three dimensions. Since the number of molecules in a cell is $n_{s, 𝒊} Δ V$ ⁠, the larger A is, the more molecules in a cell there are and the weaker the fluctuations become. However, the corresponding macroscopic system is unaffected by the value of A.

A. Schlögl model at thermodynamic equilibrium

In this section, we test the numerical schemes constructed in Sections III and IV on the Schlögl reaction-diffusion model, which is first introduced in Section II C. Simulation parameters are chosen to correspond to a system in thermodynamic equilibrium, so that the equilibrium fluctuations are Poisson and the structure factor S(k) = n_eq is constant, both with and without (i.e., diffusion only) chemical reactions. Specifically, we set the rate constants as k₁ = k₂ = k₃ = k₄ = 0.1 (see Eq. (15)), which gives n_eq = 1 and $α = 0.2 Δ t$ ⁠. We set the diffusion coefficient D = 1 and the grid spacing to unity for d = 1, 2, 3 and thus $β = Δ t$ ⁠. We consider the case where the mean number of molecules per cell is 10 by setting A = 10.

1. Continuous-time FHD equation

Prior to evaluating the different temporal integration strategies, we first focus our attention on the continuous-time discrete-space FHD equation (19) to establish a baseline for comparison and to evaluate the effectiveness of the choice of averaging function (20). To eliminate temporal integration errors, we use the EMTau scheme (25) with a very small time step size $Δ t = 10^{- 3}$ ⁠; results from the other FHD schemes are similar for sufficiently small $Δ t$ ⁠. The left panel of Fig. 3 shows the structure factors S(k) computed for N_c = 512 grid cells in one dimension for the arithmetic, geometric, and harmonic mean (AM, GM, HM) averaging functions $\tilde{n} (n_{1}, n_{2})$ ⁠. The correct flat spectrum is accurately reproduced by the modified AM averaging function (20). On the other hand, the GM and HM averaging functions give smaller S(k) at diffusion-dominated scales $k ℓ ≫ 1$ ⁠. This can be also observed from the diffusion-only system for all wavenumbers, as theoretically explained in Appendix B. Henceforth, we use the modified AM averaging function (20).

FIG. 3.

View large Download slide

Static fluctuations of the spatially discretized FHD equation (19) at thermodynamic equilibrium. (Left) Structure factors $S (\tilde{k})$ calculated by using different averaging functions $\tilde{n}$ ⁠. The solid lines depict the results from the equilibrium Schlögl model (⁠ $n_{eq} Δ V = 10$ ⁠), whereas the dotted lines are from the corresponding diffusion-only system. Note that the exact result for both systems is $S (\tilde{k}) = n_{eq}$ independent of $\tilde{k}$ ⁠, which corresponds to Poisson equilibrium fluctuations. (Right) Empirical histograms P(N) of the number of molecules per cell for the Schlögl model (red circles) and the diffusion-only system (blue crosses), computed using the arithmetic mean averaging function. For comparison, we show the correct Poisson distribution P_Poisson(N) and its Gaussian approximation P_Gauss(N). The inset shows the errors P(N) − P_Poisson(N) with error bars corresponding to two standard deviations.

The right panel of Fig. 3 shows that using the AM averaging function, the correct Poisson distribution for the number N of molecules in a cell is accurately reproduced for both reaction-diffusion and diffusion-only systems at thermodynamic equilibrium. From the equilibrium number density distribution $ρ (n)$ ⁠, we construct a discrete distribution for integer number of molecules N per cell,

P (N) = \int_{(N - \frac{1}{2}) / Δ V}^{(N + \frac{1}{2}) / Δ V} ρ (n) d n,

(36)

and compare it with a Poisson distribution P_Poisson(N) with mean $n_{eq} Δ V$ ⁠, as well as a Gaussian distribution P_Gauss(N) having the same mean and variance as P_Poisson(N). The agreement between P(N) and P_Poisson(N) is remarkable in the sense that FHD was originally proposed to account for only second moments of (small) Gaussian fluctuations. Since P_Poisson(N) is significantly different from P_Gauss(N) for $n^{eq} Δ V = 10$ ⁠, we confirm that our spatially discretized FHD equation (19) describes (large) Poisson fluctuations faithfully.

2. Time integration errors

In order to investigate time integration errors of our numerical schemes, we compare the numerical equilibrium distribution for a given time step size $Δ t$ with the target Poisson distribution P_Poisson(N) by using the following measures. First, we compute the Kullback–Leibler (KL) divergence (distance),

D_{KL} = \sum_{N = 0}^{\infty} P_{Poisson} (N) \log \frac{P_{Poisson} (N)}{P (N)} .

(37)

Second, we compute the probability of negative number densities,

P_{neg} = \int_{- \infty}^{0} ρ (n) d n .

(38)

Third, we compute the correlation coefficient $ζ$ between neighboring cells,

ζ = \frac{Cov [n_{𝒊}, n_{𝒊 \pm 𝒆_{j}}]}{Var [n]} .

(39)

Note that all three measures should be zero at thermodynamic equilibrium, as they would be for RDME.

We considered the three options for the stochastic flux amplitude ${\tilde{n}}_{s}^{∙}$ in Eqs. (30). For $ζ$ ⁠, the three options give similar values within standard errors of estimation. For D_KL and P_neg, option (30a) gives the largest values (i.e., least accurate) and option (30c) the smallest (not shown). Based on this result, we will adopt (30c) and use it for all of the simulations. Figure 4 shows how these measures deviate from zero as $Δ t$ increases for N_c = 64 cells in one dimension for the different schemes. As expected, for small values of $Δ t$ ⁠, all schemes give similar values. D_KL converges to a small value, which is consistent with the good agreement between P(N) and P_Poisson(N) seen in the right panel of Fig. 3. P_neg is observed to converge to zero as $Δ t \to 0$ ⁠, which demonstrates the effectiveness of the approach described in Section III C and agrees with the analysis in Appendix B. Also, no correlation between neighboring cells is observed for small $Δ t$ within statistical errors, which is consistent with the flat spectrum S(k) shown in the left panel of Fig. 3. As $Δ t$ approaches the explicit stability limit $Δ t_{\max}$ ⁠, rapid worsening is observed in both explicit schemes in all three measures. While P_neg and $ζ$ behave similarly for both schemes, D_KL remains small for larger values of $Δ t$ for the ExMidSSA scheme compared to the ExMidTau scheme.

FIG. 4.

View large Download slide

Deviations from the correct equilibrium distribution P_Poisson(N) as $Δ t$ increases for the four midpoint schemes applied to the one-dimensional Schlögl model at thermodynamic equilibrium. The left panel shows the KL divergence (37), the middle panel shows the probability of negative density (38), and the right panel shows the correlation coefficient (39) between neighboring cells. The red and green solid lines denote the ImMidTau and ExMidTau schemes, respectively, whereas the blue and purple dotted lines denote the ImMidSSA and ExMidSSA schemes, respectively. The arrows denote the stability limit $Δ t_{\max}$ of the explicit schemes, see Eq. (24)The error bars correspond to two standard deviations.

For both implicit schemes, it can be clearly seen that not only is the diffusion instability bypassed but also the accuracy is well maintained for large $Δ t$ ⁠. For comparable accuracy, an order of magnitude larger time step size than the explicit schemes can be chosen. The ImMidTau scheme gives smaller D_KL than the ImMidSSA scheme if $Δ t$ is smaller than a certain value. This is consistent with the observation that the former scheme has a higher temporal order of accuracy in S(k). However, due to inaccurate handling of reactions by tau leaping for large $Δ t$ ⁠, the ImMidTau scheme eventually gives larger D_KL. For P_neg and $ζ$ ⁠, a similar behavior is observed.

A similar behavior is observed for higher spatial dimensions d = 2 and d = 3 (not shown). However, for a given target accuracy tolerance, we find that a smaller time step size should be chosen, which is inversely proportional to the dimensionality d. This should not come as a surprise since the explicit stability limit $Δ t_{\max} \sim 1 / d$ ⁠, see Eq. (24). Therefore, we conclude that $Δ t$ can be chosen as an order of magnitude larger for the implicit schemes than for the explicit schemes independent of the spatial dimension. As mentioned, the computational overhead for solving linear systems can be reduced by an efficient iterative solver. Using multigrid relaxation,¹⁰³ the overall computational efficiency gain was roughly estimated to be a factor of 3. However, this factor largely depends on the problem as well as the implementation, especially on the linear solver used.

B. Turing-like pattern formation

In this section, we investigate pattern formation in the three-species Baras–Pearson–Mansour (BPM) model,^77,78

\begin{matrix} U + W \overset{k_{1}}{\to} V + W, & 2 V \underset{k_{3}}{\overset{k_{2}}{⇌}} W, \\ U \underset{k 5}{\overset{k_{4}}{⇌}} \emptyset, & V \underset{k 7}{\overset{k_{6}}{⇌}} \emptyset . \end{matrix}

(40)

We choose the rate constants so that the deterministic reaction-only system attains a limit cycle as its stable attractor, and we choose the diffusion coefficients so that a Turing-like pattern forms in the reaction-diffusion system.⁶¹ Specifically, we set $k_{1} = k_{2} = 2 \times 10^{- 4}$ ⁠, k₃ = 1, $k_{4} = 3.33 \times 10^{- 3}$ ⁠, k₅ = 16.7, $k_{6} = 3.67 \times 10^{- 2}$ ⁠, k₇ = 4.44, and D_U = 0.1, D_V = D_W = 0.01. We note that on the limit cycle, number densities of the three species oscillate in significantly different ranges: $n_{U} \in (999,2024)$ ⁠, $n_{V} \in (302,645)$ ⁠, and $n_{W} \in (18.2, 83.2)$ ⁠. For a physical domain with side lengths L_x = L_y = 32, we use three spatial resolutions with grid sizes N_c = 64², 128², and 256² cells. For the initial number densities, we choose a point on the limit cycle $(n_{U}^{0}, n_{V}^{0}, n_{W}^{0}) = (1686,534,56.4)$ and generate the initial number of molecules of each species s in each cell from a Poisson distribution with mean $n_{s}^{0} Δ V = n_{s}^{0} A Δ x Δ y$ ⁠. We use our implicit schemes in order to bypass the stiff stability limit imposed by the fast diffusion of U molecules. To obtain reference FHD results having minimal time integration errors, we use the ImMidTau scheme with $Δ t = 0.1$ ⁠. To test the importance of fluctuations, a deterministic version of the ImMidTau scheme is used with $Δ t = 0.1$ ⁠, with random initial conditions generated from a Poisson distribution corresponding to thickness A = 10.

Figure 5 shows snapshots of a final Turing pattern formed for A = 1 and 64² cells. While the pattern is qualitatively correct, the quantitative behavior of our FHD formulation may be questioned since the mean number of W molecules in a cell can be as low as 4.5 at small t in this case. To confirm that the FHD description applies even for relatively small numbers of molecules per cell, we compare the FHD results to reference RDME results obtained using the SSA/2+MN+SSA/2 scheme (A3) with $Δ t = 0.01$ ⁠. We find that the FHD reference simulations are qualitatively very similar to the RDME reference simulations over a wide range of thicknesses A, as we illustrate in Fig. 5. For our setup, after the initial formation of a disordered pattern of dots with low concentration of U molecules (blue dots in Fig. 5), the dots split and merge and diffuse to eventually form a stable regular pattern; note that the final patterns are nearly periodic lattice structures but their geometry is affected to some extent by the finite size of the domain. For A = 1, by $t = 5 \times 10^{4}$ ⁠, almost all samples had formed a steady pattern. Most samples formed a hexagonal (12 dots, see panels (a) and (c) in Fig. 5), and a few formed a monoclinic (11 dots, see panels (b) and (d) in Fig. 5) lattice of dots, for both FHD and RDME. Note that while FHD simulations using the ImMidTau scheme are equally efficient independent of A, RDME simulations become prohibitively expensive for large $A ≳ 100$ due to the very large number of U molecules (as many as $2 \times 10^{6} A$ ⁠) in the system. For weaker fluctuations, A = 1000, FHD simulations reveal that the annealing of the lattice defects takes much longer and we see several disordered or defective patterns even at $t = 5 \times 10^{4}$ (not shown). Therefore, not only do fluctuations accelerate the formation of the initial (disordered) pattern, but they also appear to accelerate the annealing of the defects.

FIG. 5.

View large Download slide

Two types of steady-state Turing-like patterns observed in the two-dimensional BPM model with 64² cells and A = 1. Snapshots of n_U are obtained at $t = 5 \times 10^{4}$ from FHD (left two panels) and RDME (right two panels) simulations with a small time step size. Panels (a) and (c) show a hexagonal structure (with 12 dots), whereas panels (b) and (d) exhibit a monoclinic structure (with 11 dots).

Since the formation of the pattern is driven by instability, it is itself a random process and a proper quantitative comparison between the different methods requires a careful statistical analysis of an ensemble of trajectories. In order to capture the time transient behavior of pattern formation, illustrated in the right panel of Fig. 6, we calculate the spatially averaged density ${\bar{n}}_{U} (t) = \frac{1}{N_{c}} \sum_{𝒊} n_{U, 𝒊} (t)$ ⁠. In the left panel, we compare sample trajectories of ${\bar{n}}_{U} (t)$ for A = 1 and A = 10 for RDME (similar results are obtained for FHD) and for deterministic reaction-diffusion. While ${\bar{n}}_{U} (t)$ initially oscillates as in the limit cycle, as the Turing-like pattern begins to form, the oscillation amplitude decays and ${\bar{n}}_{U} (t)$ eventually attains a steady value. By comparing A = 1 and 10, we see that larger fluctuations facilitate faster pattern formation, as observed in prior work⁶¹ by us and others. By comparing RDME results for thickness A = 10 with deterministic reaction-diffusion started from the same initial condition, we see that the effect of fluctuations on pattern formation is not just due to random initial conditions.

FIG. 6.

View large Download slide

Turing-like pattern formation in the two-dimensional BPM model with 64² cells. (Left) Spatially averaged density ${\bar{n}}_{U} (t)$ of species U for domain thickness A = 1 and A = 10 (RDME results), as well as deterministic reaction-diffusion started from random initial conditions corresponding to A = 10. (Right) Snapshots of n_U for A = 10 at four different times t at which ${\bar{n}}_{U} (t)$ attains a local minimum, indicated by circles in the left panel.

We generate 16 sample trajectories for each set of parameter values, fit each realization of ${\bar{n}}_{U} (t)$ using seven fitting parameters $a_{1}, \dots, a_{7}$ to

{\bar{n}}_{U} (t) = (1 - \tanh \frac{t - a_{1}}{a_{2}}) (a_{3} \sin (a_{4} t + a_{5}) + a_{6}) + a_{7},

(41)

and compare the distributions of the fitting parameters. Note that a₁ and a₇ correspond to the decay onset time and the steady spatial average density, respectively. In the left panel of Fig. 7, we compare the empirical distributions of (a₁, a₇) from the RDME and FHD results for different values of the thickness A and spatial resolutions. For each value of A, we observe that distributions obtained from different methods and/or resolutions coincide. For A = 1 and A = 10, we reconfirm that the RDME approach produces statistically very similar results for three resolutions and the FHD results are statistically indistinguishable from the RDME results. It is quite remarkable that FHD works even for A = 1 and 256² cells, which can have as low as 0.3 molecules per cell at small t, and this demonstrates the robustness of our treatment for a small number of molecules per cell. As the magnitude of the fluctuations increases (i.e., as A decreases), the pattern begins to form earlier (i.e., a₁ decreases), as already seen in Fig. 6, while the steady spatial average density a₇ becomes smaller. In addition, while the variance of a₁ does not change significantly as A varies, the variance of a₇ becomes larger for smaller A.

FIG. 7.

View large Download slide

(Left) Scatter plots of the decay onset time a₁ and the steady spatial average density a₇ for several values of the cell thickness A and two grid resolutions. RDME and FHD results are compared for A = 1 and A = 10, whereas only the FHD results (using the ImMidTau scheme) are shown for A = 100 and A = 1000 since RDME simulations are prohibitively expensive. (Right) Average values ${\bar{a}}_{7}$ of a₇ over 16 samples for the ImMidTau and ImMidSSA schemes with 64² cells as a function of time step size $Δ t$ ⁠, for cell thicknesses A of 0.1, 1, and 10. For comparison, the RDME results for $Δ t = 0.01$ are shown on the left with error bars corresponding to two standard deviations. Error bars are omitted for the implicit schemes; they are comparable to the RDME results. Note that a_i are normalized by the average values ${\bar{a}}_{i}^{\det}$ for deterministic reaction-diffusion started with random initial conditions corresponding to A = 10.

Finally, we investigate time integration errors of our implicit schemes for the Turing-like pattern formation. For the ImMidTau scheme, we increase $Δ t$ up to the stability limit arising from reactions $Δ t_{\max}^{(R)} \approx 1.3$ ⁠. The right panel of Fig. 7 shows the mean values ${\bar{a}}_{7}$ of the steady spatial average density a₇ over 16 samples versus $Δ t$ for 64² cells. While both schemes give similar values to the RDME results for small $Δ t$ ⁠, they show a different behavior for large $Δ t$ ⁠, which also depends on A. As expected, in the ImMidTau scheme, the value of ${\bar{a}}_{7}$ rapidly deviates from the RDME result as $Δ t$ approaches $Δ t_{\max}^{(R)}$ ⁠, especially if there are few molecules per cell, A = 0.1 or A = 1. On the other hand, in the ImMidSSA scheme, $Δ t$ can be increased beyond $Δ t_{\max}^{(R)}$ ⁠, and deviations from the RDME results remain small even for the smallest value of A. Hence, handling reactions by SSA not only removes the reaction stability constraint but also improves the accuracy for a small number of molecules per cell. However, it should be noted that this improvement comes at a significant computational cost, since the SSA scheme is much more expensive than tau leaping especially as the number of molecules per cell increases. Therefore, the SSA-based schemes are impractical in the regime of weak fluctuations due to poor scaling. We discuss some alternatives to SSA that may significantly improve the computational cost for weaker fluctuations in the Conclusions.

C. Front propagation

As a final example, we simulate three-dimensional front propagation in a two-species stochastic reaction-diffusion system having the following reaction network:

A \overset{k_{1}}{\to} \emptyset, 2 A + B \overset{k_{2}}{\to} 3 A, B \underset{k_{4}}{\overset{k_{3}}{⇌}} \emptyset .

(42)

This model has been proposed to reproduce axial segmentation in Ref. 7, where ISSA simulations have been performed for the one-dimensional case. Following Ref. 7, we set k₁ = 0.4, k₂ = 0.137, k₃ = 0.1, k₄ = 1 and D_A = 1, D_B = 10. For a physical domain with side lengths L_x = L_y = L_z = 512, we use 256³ cells. To initiate front propagation, we generate initial number densities as follows. We first assign to each cell i and species s a mean number density,

n_{s, 𝒊}^{0} = n_{s}^{(1)} + \frac{1}{2} (1 + \tanh \frac{r_{𝒊} - R}{ξ}) (n_{s}^{(2)} - n_{s}^{(1)}),

(43)

where $r_{𝒊}$ is the distance from the cell center to the center of the domain. This initializes a spherical region of radius R in the first uniform equilibrium state of the model, $(n_{A}^{(1)}, n_{B}^{(1)}) = (2.16, 1.35)$ ⁠, while the rest of the domain is initialized in the second uniform equilibrium state, $(n_{A}^{(2)}, n_{B}^{(2)}) = (0,10)$ ⁠, with a smooth transition region between the two states of width $\approx 2 ξ$ ⁠. Then, as in Section V B, we generate the initial number of molecules of each species in each cell from a Poisson distribution with mean $n_{s, 𝒊}^{0} Δ V = n_{s, 𝒊}^{0} A Δ x Δ y Δ z$ ⁠. We simulate the system for parameters A = 1000, R = 16, and $ξ = 4 = 2 Δ x$ using the ImMidTau scheme with $Δ t = 0.25$ ⁠. For comparison, we also simulate the corresponding deterministic system using a deterministic version of the ImMidTau scheme with the same time step size and (random) initial conditions. Simulations are performed using a parallel implementation of the algorithm using the BoxLib software framework.¹⁰⁵ We emphasize that a corresponding RDME system is too large to simulate with conventional RDME-based methods; while the total number of molecules in the system varies as the front propagates, it is of the order of 10¹² for A = 1000.

Figure 8 shows the growth of the spherical region as the more stable phase propagates into the less stable phase via a spherical traveling wave. While the phase boundary having a peak population of species A propagates, a Turing-like pattern develops behind the wave front; in one dimension this pattern is periodic and more pronounced in the presence of fluctuations.⁷ In two and three dimensions, fluctuations not only enhance the pattern but they also make it disordered, as seen in the figure by comparing the stochastic and deterministic cases. In addition, the phase boundary becomes more irregular under fluctuations. Note that the numerical results for the deterministic case are not perfectly radially symmetric not only due to the noisy initial conditions but also due to grid artifacts introduced by the standard discrete Laplacian, which is not perfectly isotropic¹⁰⁶ on length scales compared to the front width (i.e., the penetration depth); one would require an even finer grid to correct for this spatial discretization artifact.

FIG. 8.

View large Download slide

Three-dimensional FHD simulation of front propagation in a two-species stochastic reaction-diffusion model (42) using the ImMidTau scheme. The number density of species A is shown at four different times t for the stochastic reaction-diffusion system (in the top row) and the corresponding deterministic case (in the bottom row). The same initial conditions with Poisson fluctuations are used in both simulations.

VI. CONCLUSIONS

In this work, we have formulated a fluctuating hydrodynamics (FHD) model for reaction-diffusion systems and developed numerical schemes to solve the resulting stochastic ordinary differential equations (SODEs) (19) for the number densities $n_{s, 𝒊} (t)$ of chemical species in each cell. We obtained the diffusion part of the SODEs from an FHD description of a microscopic system consisting of molecules undergoing independent Brownian motions and added reactions in an equivalent manner to the reaction-diffusion master equation (RDME). We presented two implicit predictor-corrector schemes, the ImMidTau (32) and ImMidSSA (33) schemes, that treat reactions using tau leaping and SSA, respectively. In these schemes, diffusion is treated implicitly so that the stability limit imposed by fast diffusion can be bypassed and the time step size can be chosen to be significantly larger than the hopping time scale of diffusing molecules. In addition, two-stage Runge–Kutta temporal integrators are employed to improve the accuracy. To confirm the validity of our FHD formulation and demonstrate the performance of our numerical schemes, we numerically investigated not only a system at steady state (Schlögl reaction-diffusion model) but also time-dependent two-dimensional Turing-like pattern formation and three-dimensional front propagation.

Based on our analytical and numerical investigation, we conclude that the ImMidTau scheme is an efficient and robust alternative numerical method for reaction-diffusion simulations. The reason is threefold. First, the cost of the scheme does not increase for increasing number of molecules per cell (weaker fluctuations). For small numbers of molecules per cell (large fluctuations), the integer-valued RDME description is more appropriate than the continuous-range FHD description. However, by using the approach proposed in Section III C, we ensured that the FHD description remains robust and gives accurate results even for a small number of molecules per cell, as shown in Section V A. Hence, as shown in Section V B, our numerical methods can efficiently simulate reaction-diffusion systems over a broad range of relative magnitude of the fluctuations. Second, the scheme allows a significantly larger time step size without degrading accuracy compared to existing RDME-based numerical methods,^42–47 which use a fixed time step size for diffusion that is comparable to the hopping time scale. In particular, we found that the time step size could be chosen as an order of magnitude larger for the implicit schemes than for explicit methods, independent of the spatial dimension. Lastly, FHD can take advantage of the development of efficient parallel algorithms developed for computational fluid dynamics (CFD) that enable effective use of high-performance parallel architectures while providing the framework for treating more complex problems with additional physical phenomena. This enabled us to perform three-dimensional simulations of chemical front propagation involving as many as 10¹² molecules using the BoxLib CFD software framework.¹⁰⁵

The explicit tau leaping methods used here are quite simple to implement but are subject to a stability limit for fast reactions, and can lead to negative densities when fluctuations are large. While some implicit tau leaping methods have been developed, as an alternative we developed methods that use SSA for reactions. The ExMidSSA and ImMidSSA methods, however, do not scale as the fluctuations become weaker. This deficiency can be corrected by replacing SSA by a recently developed hybrid algorithm termed asynchronous tau leaping¹⁰⁷ that combines SSA and tau leaping in a dynamic manner by simulating multiple events with asynchronous time steps. Future work should develop FHD-based numerical schemes that are accurate and robust even for a small number of molecules per cell and also scale to the deterministic limit efficiently.

One of the advantages of the FHD approach for reaction-diffusion systems is its natural generalization to more complicated and realistic applications. Chemical reactions of interest usually occur in liquid solution, and often in a dense crowded environment such as the cytoplasm.^108–110 It is well-known that Brownian motion of liquid molecules or suspended macromolecules in liquids is dominated by hydrodynamic effects related to viscous dissipation.^81,110,111 This means that the diffusion model commonly used in reaction-diffusion models, including this work, which assumes that reactants are independent non-interacting Brownian walkers diffusing with a constant diffusion coefficient, does not apply in the majority of practical problems of interest. Notably, crowding or steric interactions affect the local hydrodynamic mobility of individual reactants, and hydrodynamic interactions (HI) among the diffusing reactants introduce strong correlations among the diffusive motions of the reacting particles (and also among reactants and passive crowding agents).¹¹⁰ Excluded volume due to steric repulsion introduces cross-diffusion effects, i.e., coupling between the diffusive fluxes for different species,¹¹² as well as concentration-dependent diffusion coefficients. Furthermore, it has been observed that cross-diffusion may lead to qualitatively different Turing instabilities.^113–115 Long-ranged contributions of hydrodynamic interactions can be captured by accounting for the advection of concentration fluctuations by the thermal velocity fluctuations, which follow a fluctuating Stokes equation.^81,111 Additional thermodynamic contributions to the diffusive fluxes such as cross-diffusion, barodiffusion, and thermodiffusion do not seem straightforward in the RDME but are easily included in our FHD formulation.⁶¹ In future work, we will investigate these hydrodynamic effects on stochastic reaction-diffusion phenomena.

ACKNOWLEDGMENTS

We would like to thank Jonathan Goodman for numerous informative discussions. This work was supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Applied Mathematics Program under Contract No. DE-AC02-05CH11231 and Award No. DE-SC0008271. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.

APPENDIX A: RDME-BASED SPLIT SCHEME

As a reference algorithm for stochastic reaction-diffusion simulation, we construct a numerical scheme for the RDME through a standard operator splitting technique, as done in a number of prior works.^25,44,45,47 This technique allows one to obtain a numerical scheme for the reaction-diffusion system by combining numerical methods for the diffusion-only and reaction-only systems. Here we combine multinomial diffusion sampling^43,71 for diffusion and SSA for reaction via Strang splitting.⁷⁰

One distinguishing feature of the resulting scheme, compared to exact sampling methods such as ISSA, is the use of a fixed time step size $Δ t$ for diffusion, which facilitates an efficient numerical simulation if diffusion is fast. As shown below, while there are several advantages of this scheme over ISSA, the choice of $Δ t$ is constrained as it is for explicit FHD schemes.

In Appendix A 1, after reviewing the multinomial diffusion sampling method, we present an RDME-based split scheme and discuss its advantages and disadvantages. In Appendix A 2, we present a stochastic accuracy analysis for this scheme.

1. Split scheme

In the multinomial diffusion (MN) sampling method, the number of molecules in each cell after time $Δ t$ is calculated by sampling how many molecules have moved from a cell to a neighboring cell. Here we follow the simple algorithm described by Balter and Tartakovsky⁷¹ and only allow particles to move to nearest-neighbor cells. More complicated but also more accurate algorithms that allow particles to jump to further than nearest-neighbor cells are described by Lampoudi et al.⁴³ We use the notation introduced in the body of the paper.

By denoting the number of molecules of species s diffusing from cell i to cell $𝒊^{'}$ over the time interval $Δ t$ as $N_{s, 𝒊 \to 𝒊^{'}}$ ⁠, the change in the number density $n_{s, 𝒊}$ can be expressed in terms of the sum of the inward and outward fluxes,

n_{s, 𝒊} (t + Δ t) = n_{s, 𝒊} (t) + \frac{1}{Δ V} \sum_{j = 1}^{d} (N_{s, 𝒊 + 𝒆_{j} \to 𝒊} + N_{s, 𝒊 - 𝒆_{j} \to 𝒊} - N_{s, 𝒊 \to 𝒊 + 𝒆_{j}} - N_{s, 𝒊 \to 𝒊 - 𝒆_{j}})

(A1a)

= n_{s, 𝒊} (t) + 𝔇_{𝒊} (n_{s} (t), Δ t),

(A1b)

where $n_{s} (t) = {n_{s, 𝒊} (t)}$ ⁠. For each cell i, the outward fluxes $(N_{s, 𝒊 \to 𝒊 + 𝒆_{1}}, N_{s, 𝒊 \to 𝒊 - 𝒆_{1}}, \dots, N_{s, 𝒊 \to 𝒊 + 𝒆_{d}},$ $N_{s, 𝒊 \to 𝒊 - 𝒆_{d}}, N_{s, 𝒊 \to 𝒊})$ are random variables sampled from the multinomial distribution with $\sum_{𝒊^{'}} N_{s, 𝒊 \to 𝒊^{'}} = n_{s, 𝒊} (t) Δ V$ total trials and probabilities $(p_{s}, p_{s}, \dots, p_{s}, p_{s}, 1 - 2 d p_{s})$ ⁠, where $p_{s} = D_{s} Δ t / Δ x^{2}$ ⁠, where we have assumed $Δ x_{1} = \dots = Δ x_{d} = Δ x$ ⁠.

For fast diffusion, this method becomes more efficient than treating hoping events one by one (as in ISSA). However, it is an approximate method since the actual distribution of the outward fluxes deviates from the multinomial distribution as $Δ t$ increases. In fact, $Δ t$ cannot be arbitrarily large and is limited by condition (24) because of the requirement $1 - 2 d p_{s} \geq 0$ ⁠. We also note that the number of molecules in a cell never becomes negative due to the constraint $\sum_{𝒊^{'}} N_{s, 𝒊 \to 𝒊^{'}} = n_{s, 𝒊} (t) Δ V$ ⁠. Hence, the fluxes on disjoint faces are correlated, which is different from the FHD description (18). In the deterministic limit, this scheme converges to a standard forward Euler scheme for the diffusion equation and is therefore only first-order accurate in time. In the stochastic setting, this scheme adds correlations between the fluxes through different faces of a given cell in such a way as to ensure discrete fluctuation-dissipation balance for any allowable time step size. In fact, for a system with diffusion only, this method ensures that the equilibrium fluctuations are strictly Poisson, as desired for independent Brownian walkers.

In order to handle chemical reactions, we use the SSA algorithm locally and independently in each cell, without any diffusive events. Let $ℜ_{s} (𝒏, τ)$ denote the (random) change in the number density of species s when a cell with initial state n is simulated using SSA over a time interval $τ$ ⁠. In the absence of diffusion, the SSA-based reaction scheme can be written as

n_{s, 𝒊} (t + Δ t) = n_{s, 𝒊} (t) + ℜ_{s} (𝒏_{𝒊} (t), Δ t) .

(A2)

If we combine the diffusion-only (A1) and reaction-only (A2) schemes using a Strang splitting approach,⁷⁰ we obtain the SSA/2+MN+SSA/2 scheme

n_{s, 𝐢}^{⋆} = n_{s, 𝐢}^{k} + ℜ_{s} (𝐧_{𝐢}^{k}, Δ t / 2),

(A3a)

n_{s, 𝐢}^{⋆ ⋆} = n_{s, 𝐢}^{⋆} + 𝔇_{𝒊} (n_{s}^{⋆}, Δ t),

(A3b)

n_{s, 𝐢}^{k + 1} = n_{s, 𝐢}^{⋆ ⋆} + ℜ_{s} (𝐧_{𝐢}^{⋆ ⋆}, Δ t / 2),

(A3c)

where superscripts denote time step or intermediate stage. (We note that a number of different splitting variants are possible; the version presented here gave the most accurate structure factor.) This split scheme has a number of advantages. It becomes an exact sampler (solver) for the RDME in the limit $Δ t \to 0$ ⁠, just like ISSA. It is notably more efficient than ISSA if there are many events per time interval $Δ t$ ⁠, and it can be parallelized in a straightforward manner using domain decomposition. Since the number of molecules is always a non-negative integer both in multinomial diffusion sampling and in SSA, this property is also preserved for this scheme. Moreover, since both sampling methods preserve the thermodynamic equilibrium distribution (i.e., the Poisson statistics), the split scheme also preserves it for any allowable time step size.

However, this scheme has some disadvantages. First, the time step size restriction (24) for diffusion becomes severe for fast diffusion. Second, SSA exhibits poor scalability with respect to the number of molecules in a cell. This can be resolved by replacing SSA with the tau leaping method, but the non-negativity is no longer guaranteed. Third, since the multinomial diffusion method used here is only first-order accurate, the accuracy of the scheme is first order even though Strang splitting is used. Constructing RDME-based diffusion methods that are more accurate is possible^43,46 but nontrivial and is not the focus of our work.

2. Structure factor analysis

In this section, we investigate the stochastic accuracy of the SSA/2+MN+SSA/2 scheme (A3). To this end, we consider the structure factor S(k) of an out-of-equilibrium monostable Schlögl model in one dimension (see Sections II B and II C). In the limit of many molecules per cell, an asymptotic expression of S(k) can be obtained for the scheme as a function of $Δ x$ and $Δ t$ ⁠. The multinomial fluxes can be approximated by correlated Gaussian ones and the type of analysis summarized in Appendix C can be applied; we do not give the details here for brevity. Note that we present a similar structure factor analysis for our FHD-based schemes with some background and details in Section IV B.

Figure 9 illustrates how S(k) deviates from the exact result (14) as $Δ t$ is increased to $Δ t_{\max}$ (equivalently, to $β = 0.5$ ⁠), for $α = 0.1 β$ ⁠, see Eq. (23). While S(k) is accurately reproduced at the reaction-dominated scales $k ℓ ≪ 1$ for all values of $β$ ⁠, it becomes inaccurate for smaller scales as $Δ t$ approaches $Δ t_{\max}$ ⁠. We recall that for a system at thermodynamic equilibrium, the split scheme exactly preserves the correct equilibrium distribution for any $Δ t < Δ t_{\max}$ ⁠. This property ensures that, for $α ≪ β$ ⁠, good structure factors are obtained even for systems outside of thermodynamic equilibrium, which exhibit a nonzero correlation length. For example, for $β = 0.25$ ⁠, the SSA/2+MN+SSA/2 scheme in Fig. 9 gives a notably more accurate S(k) than the FHD-based explicit schemes in Fig. 2. Hence, even though this split scheme is found to give only first-order accurate S(k), the resulting structure factor is relatively insusceptible to increasing $Δ t$ until the stability limit is approached.

FIG. 9.

View large Download slide

Structure factors $S (\tilde{k})$ obtained from the SSA/2+MN+SSA/2 scheme (A3) for the one-dimensional linearized FHD equation (12) with $\bar{Γ} / \bar{n} r = 2$ (e.g., an out-of-equilibrium monostable Schlögl model). Different values of $β$ are compared, with $α = 0.1 β$ ⁠, see Eq. (23)The exact result (14) is depicted by the dotted line. Note that the modified wavenumber $\tilde{k}$ is used, see Eq. (34)

APPENDIX B: AVERAGING FUNCTION $\tilde{𝑛}$

In this appendix, we show that the arithmetic mean should be chosen as the averaging function $\tilde{n} (n_{1}, n_{2})$ ⁠. Here we consider the diffusion-only case for a single species and investigate the equilibrium distribution of the spatially discretized FHD equation (18). Since the true equilibrium distribution for a bulk (infinite) system is known to be a product Poisson distribution from the corresponding microscopic system consisting of molecules undergoing independent Brownian motions, we choose $\tilde{n}$ so that the resulting equilibrium distribution is as close as possible to the true distribution. In addition, since the prevention of negative cell number densities is one of the essential issues for the development of a robust FHD numerical scheme, special care is taken to modify the form of $\tilde{n} (n_{1}, n_{2})$ for small values of n₁ and n₂. Here we focus on the continuous-time case and do not assume any specific temporal integrator.

The corresponding microscopic system has $N_{c} \bar{n} Δ V$ molecules, where $\bar{n}$ is the mean number density and N_c is the number of cells. The equilibrium distribution of the numbers of molecules in each cell follows the multinomial distribution with equal probabilities 1/N_c. Hence, it is straightforward to obtain the following second-order statistics of cell number density $n_{𝒊}$ ⁠:

Var [n_{𝒊}] = \frac{N_{c} - 1}{N_{c}} \frac{\bar{n}}{Δ V},

(B1)

Cov [n_{𝒊_{1}}, n_{𝒊_{2}}] = - \frac{1}{N_{c}} \frac{\bar{n}}{Δ V} for 𝒊_{1} \neq 𝒊_{2} .

(B2)

Equivalently, from Eqs. (B1) and (B2), the structure factor is also obtained as

S (𝒌) = \bar{n} for nonzero 𝒌 .

(B3)

If the FHD system (18) attains an equilibrium state, it can be shown that its second-order statistics are completely characterized by $⟨ \tilde{n} ⟩$ ⁠, which is the equilibrium average of $\tilde{n}$ over all faces,

Var [n_{𝒊}] = \frac{N_{c} - 1}{N_{c}} \frac{⟨ \tilde{n} ⟩}{Δ V},

(B4)

Cov [n_{𝒊_{1}}, n_{𝒊_{2}}] = - \frac{1}{N_{c}} \frac{⟨ \tilde{n} ⟩}{Δ V} for 𝒊_{1} \neq 𝒊_{2},

(B5)

S (𝒌) = ⟨ \tilde{n} ⟩ for nonzero 𝒌 .

(B6)

Comparing Eqs. (B4)–(B6) to the correct result Eqs. (B1)–(B3) suggests that one needs to choose $\tilde{n}$ so that $⟨ \tilde{n} ⟩$ is as close as possible to $\bar{n}$ ⁠. It is easy to see that the arithmetic mean (AM) would give the right answer: $⟨ {\tilde{n}}^{AM} ⟩ = ⟨ \frac{1}{2} (n_{1} + n_{2}) ⟩ = \bar{n}$ ⁠. On the other hand, to calculate $⟨ \tilde{n} ⟩$ for the geometric mean ${\tilde{n}}^{GM} = \sqrt{n_{1} n_{2}}$ or the harmonic mean ${\tilde{n}}^{HM} = 2 / (n_{1}^{- 1} + n_{2}^{- 1})$ ⁠, one needs to know the equilibrium distribution $ρ (n_{1}, n_{2})$ of two neighboring cells. However, under the reasonable assumption that all three averaging functions give similar distributions $ρ (n_{1}, n_{2})$ allowing only non-negative number densities, it can be easily shown that $⟨ {\tilde{n}}^{HM} ⟩ \leq ⟨ {\tilde{n}}^{GM} ⟩ \leq ⟨ {\tilde{n}}^{AM} ⟩$ from the well-known inequalities among the Pythagorean means. In fact, in Fig. 3, this ordering is observed from the structure factor of the diffusion-only system, see Eq. (B6). Hence, we conclude that the arithmetic mean is the right choice for $\tilde{n}$ ⁠.

However, if the arithmetic mean is employed without modification, Eq. (18) does not attain an equilibrium state. This is because almost surely at some point on some grid face we will have $\tilde{n} < 0$ so that the stochastic diffusive flux becomes undefined. The non-negativity of cell number densities is guaranteed if the stochastic diffusive flux through a face is turned off when the number density of either cell sharing the face becomes zero. Specifically, under some technical assumptions, it can be proven that if $\tilde{n} (n_{1}, n_{2})$ is a non-negative function that satisfies

\tilde{n} (n_{1}, n_{2}) = 0 for n_{1} \leq 0 or n_{2} \leq 0,

(B7)

then the number density of each cell never becomes negative. The non-negative arithmetic mean averaging function,

\tilde{n} (n_{1}, n_{2}) = {\begin{matrix} \frac{1}{2} (n_{1} + n_{2}) & if n_{1} > 0 and 𝑛_{2} > 0, \\ 0 & otherwise, \end{matrix}

(B8)

does not allow negative density (i.e., $ρ (n) = 0$ for $n < 0$ ⁠). However, due to the discontinuity of $\tilde{n}$ at n₁ = 0 or n₂ = 0, $ρ (n)$ does not decrease to zero as n becomes zero (i.e., ${lim}_{n \to 0^{+}} ρ (n) > 0$ ⁠) and a delta function is formed at n = 0, see Fig. 10.

FIG. 10.

View large Download slide

The equilibrium cell number density distributions $ρ (n)$ near n = 0 obtained from the arithmetic mean averaging function (B8) (depicted by the red solid line), which uses the discontinuous Heaviside function H, and Eq. (20) (depicted by the blue dashed line), which uses the smoothed Heaviside functions H₀The results are obtained from the one-dimensional diffusion-only FHD system (18) having $D = \bar{n} = Δ x = 1$ ⁠, A = 5, and N_c = 512 by using the EMTau scheme (25) with $Δ t = 10^{- 3}$

To understand this behavior, we consider numerically integrating Eq. (18) with a small time step size $Δ t > 0$ ⁠. Even if the density in a given cell n₁(t) has a small positive value, $\tilde{n}$ can be large if the density in the neighboring cell n₂(t) is large. In this case, $n_{1} (t + Δ t)$ can become negative due to the stochastic diffusive flux. But, once n₁ has become negative, the stochastic diffusive flux is turned off, and due to the deterministic diffusion, n₁ increases and becomes positive again. Hence, as shown in Fig. 10, $ρ (n)$ attains a peak in the negative density region near n = 0. The width and height of the peak are proportional to $Δ t^{1 / 2}$ and $Δ t^{- 1 / 2}$ ⁠, respectively, and the peak becomes a delta function in the limit $Δ t \to 0$ ⁠.

We note that the averaging in Eq. (B8) can be expressed as $\tilde{n} = \frac{1}{2} (n_{1} + n_{2}) H (n_{1} Δ V) H (n_{2} Δ V)$ ⁠, where H is the Heaviside function. To avoid a discontinuity in $\tilde{n}$ ⁠, we can use a smoothed Heaviside function to arrive at Eq. (20). Note that the smoothing is based on the number of molecules in a cell (⁠ $N = n Δ V$ ⁠) and the smoothing region $0 \leq N \leq 1$ is chosen so that the stochastic diffusive flux is modified only when there is less than one molecule in a cell. In Fig. 10, we show the distribution $ρ (n)$ near n = 0 obtained by the averaging function (20), for a rather small mean number of molecules per cell, $\bar{n} Δ V = 5$ ⁠. With the use of a smoothed Heaviside function, the spurious delta function at n = 0 as $Δ t \to 0$ is removed, and the probability of negative density is greatly reduced for small $Δ t$ ⁠.

APPENDIX C: LINEARIZED EQUATION ANALYSIS

In this appendix, we summarize how the discrete structure factor is obtained as a function of $Δ x$ and $Δ t$ when a given spatiotemporal discretization is applied to the linearized FHD equation (12), following the Fourier-space analysis developed in Ref. 63. For simplicity, we consider here the one-dimensional case.

Applying the spatial discretization given in Section III to Eq. (12) and taking a discrete Fourier transform, we obtain an Ornstein–Uhlenbeck equation for the Fourier coefficient $δ {\hat{n}}_{k} (t)$ ⁠,

\frac{d}{d t} δ {\hat{n}}_{k} = - D {\tilde{k}}^{2} δ {\hat{n}}_{k} + \sqrt{\frac{2 D \bar{n} {\tilde{k}}^{2}}{V}} W_{k}^{(D)} - r δ {\hat{n}}_{k} + \sqrt{\frac{2 \bar{Γ}}{V}} W_{k}^{(R)},

(C1)

where the modified wavenumber $\tilde{k}$ is defined in Eq. (34), and $W_{k}^{(D)} (t)$ and $W_{k}^{(R)} (t)$ are independent standard GWN processes. Compared to the continuous-space case (see Eq. (13)), k is replaced by $\tilde{k}$ due to the discrete Laplacian and divergence operators. Note that in this linearized analysis, Poisson processes have been replaced by Gaussian ones with the mean evaluated at the ensemble average (i.e., macroscopic) values of the density. For convenience, we have replaced complex-valued GWN processes by real-valued ones having the same noise intensities.

When the EMTau scheme (25) is used to solve Eq. (C1), we have the recursion

δ {\hat{n}}_{k} (t + Δ t) = [1 - (D {\tilde{k}}^{2} + r) Δ t] δ {\hat{n}}_{k} (t) + \sqrt{\frac{2 D \bar{n} {\tilde{k}}^{2} Δ t}{V}} W_{1} + \sqrt{\frac{2 \bar{Γ} Δ t}{V}} W_{2},

(C2)

where W₁ and W₂ are independent standard normal random variables. We write this as

δ {\hat{n}}_{k} (t + Δ t) \overset{d}{=} M_{k} δ {\hat{n}}_{k} (t) + N_{k} W,

(C3)

where $\overset{d}{=}$ denotes being equal in distribution. For any given temporal integrator, we can straightforwardly obtain analytic expressions for M_k and $N_{k} N_{k}^{*}$ ⁠. For example, for the EMTau scheme,

M_{k} = 1 - (D {\tilde{k}}^{2} + r) Δ t,

(C4a)

N_{k} N_{k}^{*} = \frac{2}{V} (D \bar{n} {\tilde{k}}^{2} + \bar{Γ}) Δ t .

(C4b)

The covariance of the noise $N_{k} N_{k}^{*}$ for multinomial diffusion (A1) can most easily be obtained from Eq. (C8) from the observation that, in the absence of reactions, the exact structure factor $S (k) = \bar{n}$ is obtained for any stable time step.

A similar procedure is applicable to numerical schemes having SSA for reactions. Since SSA is an exact integrator, the linearized reaction part in Eq. (C1) is exactly solved. For example, for the EM-SSA scheme (26), we have (cf. Eq. (C2))

δ {\hat{n}}_{k} (t + Δ t) = [- D {\tilde{k}}^{2} Δ t + e^{- r Δ t}] δ {\hat{n}}_{k} (t) + \sqrt{\frac{2 D \bar{n} {\tilde{k}}^{2} Δ t}{V}} W_{1} + \sqrt{\frac{\bar{Γ} (1 - e^{- 2 r Δ t})}{r V}} W_{2} .

(C5)

While the expressions of M_k and $N_{k} N_{k}^{*}$ become complicated for the predictor-corrector midpoint schemes, a theoretical analysis is still tractable with the help of symbolic algebra tools.

By calculating M_k and $N_{k} N_{k}^{*}$ for a given numerical scheme, the stability condition and the structure factor can be obtained as follows. From the condition that the amplification factor M_k should satisfy

| M_{k} | \leq 1 for all k,

(C6)

the stability condition is obtained. For the EMTau scheme, we obtain

\frac{D Δ t}{Δ x^{2}} + \frac{r Δ t}{4} \leq \frac{1}{2} .

(C7)

The analytic expression for $S (k) = V ⟨ δ {\hat{n}}_{k} δ {\hat{n}}_{k}^{*} ⟩$ can be calculated from

S (k) = \frac{V N_{k} N_{k}^{*}}{1 - M_{k}^{2}},

(C8)

which is obtained from the time invariance relation $⟨ δ {\hat{n}}_{k} (t) δ {\hat{n}}_{k}^{*} (t) ⟩ = ⟨ δ {\hat{n}}_{k} (t + Δ t) δ {\hat{n}}_{k}^{*} (t + Δ t) ⟩$ and Eq. (C3).⁶³ From the analytical expressions of S(k) and $S_{0} (k) = {lim}_{Δ t \to 0}$ $S (k)$ ⁠, [S(k_max) − S₀(k_max)]/S₀(k_max) is easily obtained, the series expansion of which for small $α$ (and fixed $β$ ⁠) gives Eq. (35) for the ImMidTau scheme. Similarly, series expansions for small $Δ t$ (i.e., fixed ratio $α / β$ ⁠) reveal the temporal order of accuracy of S(k).

REFERENCES

1.

D. A.

Kessler

and

H.

Levine

, “

Fluctuation-induced diffusive instabilities

,”

Nature

394

,

556

–

558

(

1998

).

https://doi.org/10.1038/29020

Google Scholar

Crossref

2.

E.

Khain

,

Y. T.

Lin

, and

L. M.

Sander

, “

Fluctuations and stability in front propagation

,”

EPL

93

,

28001

(

2011

).

https://doi.org/10.1209/0295-5075/93/28001

Google Scholar

Crossref

3.

A. A.

Winkler

and

E.

Frey

, “

Validity of the law of mass action in three-dimensional coagulation processes

,”

Phys. Rev. Lett.

108

,

108301

(

2012

).

https://doi.org/10.1103/physrevlett.108.108301

Google Scholar

Crossref

PubMed

4.

I. V.

Gopich

,

A. A.

Ovchinnikov

, and

A.

Szabo

, “

Long-time tails in the kinetics of reversible bimolecular reactions

,”

Phys. Rev. Lett.

86

,

922

–

925

(

2001

).

https://doi.org/10.1103/physrevlett.86.922

Google Scholar

Crossref

PubMed

5.

Y.

Togashi

and

K.

Kaneko

, “

Molecular discreteness in reaction-diffusion systems yields steady states not seen in the continuum limit.

,”

Phys. Rev. E

70

,

020901(R)

(

2004

).

https://doi.org/10.1103/physreve.70.020901

Google Scholar

Crossref

6.

H.

Wang

,

Z.

Fu

,

X.

Xu

, and

Q.

Ouyang

, “

Pattern formation induced by internal microscopic fluctuations

,”

J. Phys. Chem. A

111

,

1265

–

1270

(

2007

).

https://doi.org/10.1021/jp0674064

Google Scholar

Crossref

PubMed

7.

A.

Lemarchand

and

B.

Nowakowski

, “

Do the internal fluctuations blur or enhance axial segmentation?

,”

EPL

94

,

48004

(

2011

).

https://doi.org/10.1209/0295-5075/94/48004

Google Scholar

Crossref

8.

P.

Peeters

,

F.

Baras

, and

G.

Nicolis

, “

Sensitivity of the chlorite-thiosulfate system: A stochastic approach

,”

J. Chem. Phys.

93

,

7321

–

7328

(

1990

).

https://doi.org/10.1063/1.459405

Google Scholar

Crossref

9.

M.

Malek Mansour

,

J.

Dethier

, and

F.

Baras

, “

Onset of homogeneous oscillations in reactive systems

,”

J. Chem. Phys.

114

,

9265

–

9275

(

2001

).

https://doi.org/10.1063/1.1367389

Google Scholar

Crossref

10.

N.

Fedoroff

and

W.

Fontana

, “

Small numbers of big molecules

,”

Science

297

,

1129

–

1131

(

2002

).

https://doi.org/10.1126/science.1075988

Google Scholar

Crossref

PubMed

11.

K.

Doubrovinski

and

M.

Howard

, “

Stochastic model for Soj relocation dynamics in bacillus subtilis

,”

Proc. Natl. Acad. Sci. U. S. A.

102

,

9808

–

9813

(

2005

).

https://doi.org/10.1073/pnas.0500529102

Google Scholar

Crossref

PubMed

12.

D.

Fange

and

J.

Elf

, “

Noise-induced min phenotypes in E. coli

,”

PLoS Comput. Biol.

2

,

e80

(

2006

).

https://doi.org/10.1371/journal.pcbi.0020080

Google Scholar

Crossref

PubMed

13.

J. S.

van Zon

,

M. J.

Morelli

,

S.

Tănase-Nicola

, and

P. R.

ten Wolde

, “

Diffusion of transcription factors can drastically enhance the noise in gene expression

,”

Biophys. J.

91

,

4350

–

4367

(

2006

).

https://doi.org/10.1529/biophysj.106.086157

Google Scholar

Crossref

PubMed

14.

M. B.

Flegg

,

S.

Rüdiger

, and

R.

Erban

, “

Diffusive spatio-temporal noise in a first-passage time model for intracellular calcium release

,”

J. Chem. Phys.

138

,

154103

(

2013

).

https://doi.org/10.1063/1.4796417

Google Scholar

Crossref

PubMed

15.

K.

Burrage

,

P. M.

Burrage

,

A.

Leier

,

T.

Marquez-Lago

, and

D. V.

Nicolau

, Jr., in

Stochastic Simulation for Spatial Modelling of Dynamic Processes in a Living Cell

, Design and Analysis of Biomolecular Circuits, edited by

H.

Koeppl

,

G.

Setti

,

M.

di Bernardo

, and

D.

Densmore

(

Springer

,

New York

,

2011

), pp.

43

–

62

.

Google Scholar

16.

A.

Mahmutovic

,

D.

Fange

,

O. G.

Berge

, and

J.

Elf

, “

Lost in presumption: Stochastic reactions in spatial models

,”

Nat. Methods

9

,

1163

–

1166

(

2012

).

https://doi.org/10.1038/nmeth.2253

Google Scholar

Crossref

PubMed

17.

M.

von Smoluchowski

, “

Versuch einer mathematischen theorie der koagulationskinetik kolloider lösungen

,”

Z. Phys. Chem.

92

,

129

–

168

(

1917

).

Google Scholar

18.

A.

Donev

,

V. V.

Bulatov

,

T.

Oppelstrup

,

G. H.

Gilmer

,

B.

Sadigh

, and

M. H.

Kalos

, “

A first-passage kinetic Monte Carlo algorithm for complex diffusion-reaction systems

,”

J. Comput. Phys.

229

,

3214

–

3236

(

2010

).

https://doi.org/10.1016/j.jcp.2009.12.038

Google Scholar

Crossref

19.

K.

Takahashi

,

S.

Tănase-Nicola

, and

P. R.

ten Wolde

, “

Spatio-temporal correlations can drastically change the response of a MAPK pathway

,”

Proc. Natl. Acad. Sci. U. S. A.

107

,

2473

–

2478

(

2010

).

https://doi.org/10.1073/pnas.0906885107

Google Scholar

Crossref

PubMed

20.

A. J.

Mauro

,

J. K.

Sigurdsson

,

J.

Shrake

, and

P. J.

Atzberger

, “

A first-passage kinetic Monte Carlo method for reaction-drift-diffusion processes

,”

J. Comput. Phys.

259

,

536

–

567

(

2014

).

https://doi.org/10.1016/j.jcp.2013.12.023

Google Scholar

Crossref

21.

S. S.

Andrews

and

D.

Bray

, “

Stochastic simulation of chemical reactions with spatial resolution and single molecule detail

,”

Phys. Biol.

1

,

137

–

151

(

2004

).

https://doi.org/10.1088/1478-3967/1/3/001

Google Scholar

Crossref

PubMed

22.

J.

Lipková

,

K. C.

Zygalakis

,

S. J.

Chapman

, and

R.

Erban

, “

Analysis of Brownian dynamics simulations of reversible bimolecular reactions

,”

SIAM J. Appl. Math.

71

,

714

–

730

(

2011

).

https://doi.org/10.1137/100794213

Google Scholar

Crossref

23.

P.

Dziekan

,

A.

Lemarchand

, and

B.

Nowakowski

, “

Particle dynamics simulations of Turing patterns

,”

J. Chem. Phys.

137

,

074107

(

2012

).

https://doi.org/10.1063/1.4743983

Google Scholar

Crossref

PubMed

24.

P.

Dziekan

,

J. S.

Hansen

, and

B.

Nowakowski

, “

Nanoscale Turing structures

,”

J. Chem. Phys.

141

,

124106

(

2014

).

https://doi.org/10.1063/1.4895907

Google Scholar

Crossref

PubMed

25.

T.

Choi

,

M. R.

Maurya

,

D. M.

Tartakovsky

, and

S.

Subramaniam

, “

Stochastic operator-splitting method for reaction-diffusion systems

,”

J. Chem. Phys.

137

,

184102

(

2012

).

https://doi.org/10.1063/1.4764108

Google Scholar

Crossref

PubMed

26.

M. B.

Flegg

,

S. J.

Chapman

, and

R.

Erban

, “

The two-regime method for optimizing stochastic reaction-diffusion simulations

,”

J. R. Soc., Interface

9

,

859

–

868

(

2012

).

https://doi.org/10.1098/rsif.2011.0574

Google Scholar

Crossref

PubMed

27.

A.

Hellander

,

S.

Hellander

, and

P.

Lötstedt

, “

Coupled mesoscopic and microscopic simulation of stochastic reaction-diffusion processes in mixed dimensions

,”

Multiscale Model. Simul.

10

,

585

–

611

(

2012

).

https://doi.org/10.1137/110832148

Google Scholar

Crossref

28.

D. T.

Gillespie

, “

A general method for numerically simulating the stochastic time evolution of coupled chemical reactions

,”

J. Comput. Phys.

22

,

403

(

1976

).

https://doi.org/10.1016/0021-9991(76)90041-3

Google Scholar

Crossref

29.

D. T.

Gillespie

, “

Approximate accelerated stochastic simulation of chemically reacting systems

,”

J. Chem. Phys.

115

,

1716

(

2001

).

https://doi.org/10.1063/1.1378322

Google Scholar

Crossref

30.

D. T.

Gillespie

,

A.

Hellander

, and

L. R.

Petzold

, “

Perspective: Stochastic algorithms for chemical kinetics

,”

J. Chem. Phys.

138

,

170901

(

2013

).

https://doi.org/10.1063/1.4801941

Google Scholar

Crossref

PubMed

31.

T.

Székeley

, Jr. and

K.

Burrage

, “

Stochastic simulation in systems biology

,”

Comput. Struct. Biotechnol. J.

12

,

14

–

25

(

2014

).

https://doi.org/10.1016/j.csbj.2014.10.003

Google Scholar

Crossref

PubMed

32.

C. W.

Gardiner

,

Handbook of Stochastic Methods

, 3rd ed. (

Springer

,

New York

,

2004

).

Google Scholar

Crossref

33.

R.

Erban

,

S. J.

Chapman

, and

P. K.

Maini

, “

A practical guide to stochastic simulations of reaction-diffusion processes

,” eprint arXiv: 0704.1908.

34.

M.

Malek-Mansour

and

J.

Houard

, “

A new approximation scheme for the study of fluctuations in nonuniform nonequilibrium systems

,”

Phys. Lett. A

70

,

366

–

368

(

1979

).

https://doi.org/10.1016/0375-9601(79)90332-3

Google Scholar

Crossref

35.

J.

Elf

and

M.

Ehrenberg

, “

Spontaneous separation of bi-stable biochemical systems into spatial domains of opposite phases

,”

Syst. Biol.

1

,

230

–

236

(

2004

).

https://doi.org/10.1049/sb:20045021

Google Scholar

Crossref

36.

J.

Luis Hita

and

J. M.

Ortiz de Zárate

, “

Spatial correlations in nonequilibrium reaction-diffusion problems by the Gillespie algorithm

,”

Phys. Rev. E

87

,

052802

(

2013

).

https://doi.org/10.1103/physreve.87.052802

Google Scholar

Crossref

37.

B.

Wang

,

B.

Hou

,

F.

Xing

, and

Y.

Yao

, “

Abstract next subvolume method: A logical process-based approach for spatial stochastic simulation of chemical reactions

,”

Comput. Biol. Chem.

35

,

193

–

198

(

2011

).

https://doi.org/10.1016/j.compbiolchem.2011.05.001

Google Scholar

Crossref

PubMed

38.

T. T.

Marquez-Lago

and

K.

Burrage

, “

Binomial tau-leap spatial stochastic simulation algorithm for applications in chemical kinetics

,”

J. Chem. Phys.

127

,

104101

(

2007

).

https://doi.org/10.1063/1.2771548

Google Scholar

Crossref

PubMed

39.

D.

Rossinelli

,

B.

Bayati

, and

P.

Koumoutsakos

, “

Accelerated stochastic and hybrid methods for spatial simulations of reaction-diffusion systems

,”

Chem. Phys. Lett.

451

,

136

–

140

(

2008

).

https://doi.org/10.1016/j.cplett.2007.11.055

Google Scholar

Crossref

40.

K. A.

Iyengar

,

L. A.

Harris

, and

P.

Clancy

, “

Accurate implementation of leaping in space: The spatial partitioned-leaping algorithm

,”

J. Chem. Phys.

132

,

094101

(

2010

).

https://doi.org/10.1063/1.3310808

Google Scholar

Crossref

PubMed

41.

J.

Fu

,

S.

Wu

,

H.

Li

, and

L. R.

Petzold

, “

The time dependent propensity function for acceleration of spatial stochastic simulation of reaction-diffusion systems

,”

J. Comput. Phys.

274

,

524

–

549

(

2014

).

https://doi.org/10.1016/j.jcp.2014.06.025

Google Scholar

Crossref

PubMed

42.

J. V.

Rodríguez

,

J. A.

Kaandorp

,

M.

Dobrzyński

, and

J. G.

Blom

, “

Spatial stochastic modelling of the phosphoenolpyruvate-dependent phosphotransferase (PTS) pathway in Escherichia coli

,”

Bioinformatics

22

,

1895

–

1901

(

2006

).

https://doi.org/10.1093/bioinformatics/btl271

Google Scholar

Crossref

PubMed

43.

S.

Lampoudi

,

D. T.

Gillespie

, and

L. R.

Petzold

, “

The multinomial simulation algorithm for discrete stochastic simulation of reaction-diffusion systems

,”

J. Chem. Phys.

130

,

094104

(

2009

).

https://doi.org/10.1063/1.3074302

Google Scholar

Crossref

PubMed

44.

S.

Engblom

,

L.

Ferm

,

A.

Hellander

, and

P.

Lötstedt

, “

Simulation of stochastic reaction-diffusion processes on unstructured meshes

,”

SIAM J. Sci. Comput.

31

,

1774

–

1797

(

2009

).

https://doi.org/10.1137/080721388

Google Scholar

Crossref

45.

L.

Ferm

,

A.

Hellander

, and

P.

Lötstedt

, “

An adaptive algorithm for simulation of stochastic reaction-diffusion processes

,”

J. Comput. Phys.

229

,

343

–

360

(

2010

).

https://doi.org/10.1016/j.jcp.2009.09.030

Google Scholar

Crossref

46.

B.

Drawert

,

M. J.

Lawson

,

L.

Petzold

, and

M.

Khammash

, “

The diffusive finite state projection algorithm for efficient simulation of the stochastic reaction-diffusion master equation

,”

J. Chem. Phys.

132

,

074101

(

2010

).

https://doi.org/10.1063/1.3310809

Google Scholar

Crossref

PubMed

47.

A.

Hellander

,

M. J.

Lawson

,

B.

Drawert

, and

L.

Petzold

, “

Local error estimates for adaptive simulation of the reaction-diffusion master equation via operator splitting

,”

J. Comput. Phys.

266

,

89

–

100

(

2014

).

https://doi.org/10.1016/j.jcp.2014.02.004

Google Scholar

Crossref

PubMed

48.

H. C.

Öttinger

,

Beyond Equilibrium Thermodynamics

(

John Wiley & Sons

,

2005

).

Google Scholar

Crossref

49.

L. D.

Landau

and

E. M.

Lifshitz

, in

Fluid Mechanics

, Course of Theoretical Physics, 1st ed. (

Pergamon Press

,

1959

), Vol. 6, Chap. XVII, pp.

523

–

529

.

Google Scholar

Crossref

50.

P.

Español

, “

Stochastic differential equations for non-linear hydrodynamics

,”

Physica A

248

,

77

–

96

(

1998

).

https://doi.org/10.1016/s0378-4371(97)00461-5

Google Scholar

Crossref

51.

P.

Español

,

J. G.

Anero

, and

I.

Zúñga

, “

Microscopic derivation of discrete hydrodynamics

,”

J. Chem. Phys.

131

,

244117

(

2009

).

https://doi.org/10.1063/1.3274222

Google Scholar

Crossref

PubMed

52.

N. K.

Voulgarakis

and

J.-W.

Chu

, “

Bridging fluctuating hydrodynamics and molecular dynamics simulations of fluids

,”

J. Chem. Phys.

130

,

134111

(

2009

).

https://doi.org/10.1063/1.3106717

Google Scholar

Crossref

PubMed

53.

J. M.

Ortiz de Zárate

and

J. V.

Sengers

,

Hydrodynamic Fluctuations in Fluids and Fluid Mixtures

(

Elsevier

,

2006

).

Google Scholar

54.

A.

Donev

,

J. B.

Bell

,

A.

de la Fuente

, and

A. L.

Garcia

, “

Diffusive transport by thermal velocity fluctuations

,”

Phys. Rev. Lett.

106

,

204501

(

2011

).

https://doi.org/10.1103/physrevlett.106.204501

Google Scholar

Crossref

PubMed

55.

A.

Donev

,

J. B.

Bell

,

A.

de la Fuente

, and

A. L.

Garcia

, “

Enhancement of diffusive transport by non-equilibrium thermal fluctuations

,”

J. Stat. Mech.: Theory Exp.

2011

,

P06014

.

https://doi.org/10.1088/1742-5468/2011/06/p06014

Crossref

56.

A.

Donev

,

A.

Nonaka

,

Y.

Sun

,

T. G.

Fai

,

A. L.

Garcia

, and

J. B.

Bell

, “

Low Mach number fluctuating hydrodynamics of diffusively mixing fluids

,”

Commun. Appl. Math. Comput. Sci.

9

,

47

–

105

(

2014

).

https://doi.org/10.2140/camcos.2014.9.47

Google Scholar

Crossref

57.

A.

Nonaka

,

Y.

Sun

,

J. B.

Bell

, and

A.

Donev

, “

Low Mach number fluctuating hydrodynamics of binary liquid mixtures

,”

Commun. Appl. Math. Comput. Sci.

10

,

163

–

204

(

2015

).

https://doi.org/10.2140/camcos.2015.10.163

Google Scholar

Crossref

58.

A.

Chaudhri

,

J. B.

Bell

,

A. L.

Garica

, and

A.

Donev

, “

Modeling multiphase flow using fluctuating hydrodynamics

,”

Phys. Rev. E

90

,

033014

(

2014

).

https://doi.org/10.1103/physreve.90.033014

Google Scholar

Crossref

59.

K.

Balakrishnan

,

A. L.

Garcia

,

A.

Donev

, and

J. B.

Bell

, “

Fluctuating hydrodynamics of multispecies nonreactive mixtures

,”

Phys. Rev. E

89

,

013017

(

2014

).

https://doi.org/10.1103/physreve.89.013017

Google Scholar

Crossref

60.

A.

Donev

,

A.

Nonaka

,

A. K.

Bhattacharjee

,

A. L.

Garcia

, and

J. B.

Bell

, “

Low Mach number fluctuating hydrodynamics of multispecies liquid mixture

,”

Phys. Fluids

27

,

037103

(

2015

).

https://doi.org/10.1063/1.4913571

Google Scholar

Crossref

61.

A. K.

Bhattacharjee

,

K.

Balakrishnan

,

A. L.

Garcia

,

J. B.

Bell

, and

A.

Donev

, “

Fluctuating hydrodynamics of multi-species reactive mixtures

,”

J. Chem. Phys.

142

,

224107

(

2015

).

https://doi.org/10.1063/1.4922308

Google Scholar

Crossref

PubMed

62.

J.-P.

Péraud

,

A.

Nonaka

,

A.

Chaudhri

,

J. B.

Bell

,

A.

Donev

, and

A. L.

Garcia

, “

Low Mach number fluctuating hydrodynamics for electrolytes

,”

Phys. Rev. Fluids

1

,

074103

(

2016

).

https://doi.org/10.1103/physrevfluids.1.074103

Google Scholar

Crossref

63.

A.

Donev

,

E.

Vanden-Eijnden

,

A.

Garcia

, and

J.

Bell

, “

On the accuracy of finite-volume schemes for fluctuating hydrodynamics

,”

Commun. Appl. Math. Comput. Sci.

5

,

149

–

197

(

2010

).

https://doi.org/10.2140/camcos.2010.5.149

Google Scholar

Crossref

64.

S.

Delong

,

B. E.

Griffith

,

E.

Vanden-Eijnden

, and

A.

Donev

, “

Temporal integrators for fluctuating hydrodynamics

,”

Phys. Rev. E

87

,

033302

(

2013

).

https://doi.org/10.1103/physreve.87.033302

Google Scholar

Crossref

65.

Y.

Hu

,

T.

Li

, and

B.

Min

, “

A weak second order tau-leaping method for chemical kinetic systems

,”

J. Chem. Phys.

135

,

024113

(

2011

).

https://doi.org/10.1063/1.3609119

Google Scholar

Crossref

PubMed

66.

D. F.

Anderson

and

M.

Koyama

, “

Weak error analysis of numerical methods for stochastic models of population processes

,”

Multiscale Model. Simul.

10

,

1493

–

1524

(

2012

).

https://doi.org/10.1137/110849699

Google Scholar

Crossref

67.

P. J.

Atzberger

, “

Spatially adaptive stochastic numerical methods for intrinsic fluctuations in reaction-diffusion systems

,”

J. Comput. Phys.

229

,

3474

–

3501

(

2010

).

https://doi.org/10.1016/j.jcp.2010.01.012

Google Scholar

Crossref

68.

J. M.

Oritz de Zárate

,

J.

Luis Hita

, and

J. V.

Sengers

, “

Fluctuating hydrodynamics and concentration fluctuations in ternary mixtures

,”

C. R. Mec.

341

,

399

–

404

(

2013

).

https://doi.org/10.1016/j.crme.2013.01.018

Google Scholar

Crossref

69.

S. P.

Thampi

,

I.

Pagonabarraga

, and

R.

Adhikari

, “

Lattice-Boltzmann-Langevin simulations of binary mixtures

,”

Phys. Rev. E

84

,

046709

(

2011

).

https://doi.org/10.1103/physreve.84.046709

Google Scholar

Crossref

70.

G.

Strang

, “

On the construction and comparison of difference schemes

,”

SIAM J. Numer. Anal.

5

,

506

–

517

(

1968

).

https://doi.org/10.1137/0705041

Google Scholar

Crossref

71.

A.

Balter

and

A. M.

Tartakovsky

, “

Multinomial diffusion equation

,”

Phys. Rev. E

83

,

061143

(

2011

).

https://doi.org/10.1103/physreve.83.061143

Google Scholar

Crossref

72.

A.

Lemarchand

,

A.

Lesne

, and

M.

Mareshal

, “

Langevin approach to a chemical wave front: Selection of the propagation velocity in the presence of internal noise

,”

Phys. Rev. E

51

,

4457

–

4465

(

1995

).

https://doi.org/10.1103/physreve.51.4457

Google Scholar

Crossref

73.

M. A.

Karzazi

,

A.

Lemarchand

, and

M.

Mareschal

, “

Fluctuation effects on chemical wave fronts

,”

Phys. Rev. E

54

,

4888

–

4975

(

1996

).

https://doi.org/10.1103/physreve.54.4888

Google Scholar

Crossref

74.

S.

Delong

,

Y.

Sun

,

B. E.

Griffith

,

E.

Vanden-Eijnden

, and

A.

Donev

, “

Multiscale temporal integrators for fluctuating hydrodynamics

,”

Phys. Rev. E

90

,

063312

(

2014

).

https://doi.org/10.1103/physreve.90.063312

Google Scholar

Crossref

75.

F.

Schlögl

, “

Chemical reaction models for non-equilibrium phase transitions

,”

Z. Phys.

253

,

147

–

161

(

1972

).

https://doi.org/10.1007/bf01379769

Google Scholar

Crossref

76.

M.

Vellela

and

H.

Qian

, “

Stochastic dynamics and non-equilibrium thermodynamics of a bistable chemical system: The Schlögl model revisited

,”

J. R. Soc., Interface

6

,

925

–

940

(

2009

).

https://doi.org/10.1098/rsif.2008.0476

Google Scholar

Crossref

PubMed

77.

F.

Baras

,

J. E.

Pearson

, and

M.

Malek Mansour

, “

Microscopic simulation of chemical oscillations in homogeneous systems

,”

J. Chem. Phys.

93

,

5747

–

5750

(

1990

).

https://doi.org/10.1063/1.459568

Google Scholar

Crossref

78.

F.

Baras

,

M.

Malek Mansour

, and

J. E.

Pearson

, “

Microscopic simulation of chemical bistability in homogeneous systems

,”

J. Chem. Phys.

105

,

8257

–

8261

(

1996

).

https://doi.org/10.1063/1.472679

Google Scholar

Crossref

79.

D. S.

Dean

, “

Langevin equation for the density of a system of interacting Langevin processes

,”

J. Phys. A: Math. Gen.

29

,

L613

–

L617

(

1996

).

https://doi.org/10.1088/0305-4470/29/24/001

Google Scholar

Crossref

80.

G. D. C.

Kuiken

,

Thermodynamics of Irreversible Processes: Applications to Diffusion and Rheology

(

Wiley

,

1994

).

Google Scholar

81.

A.

Donev

and

E.

Vanden-Eijnden

, “

Dynamic density functional theory with hydrodynamic interactions and fluctuations

,”

J. Chem. Phys.

140

,

234115

(

2014

).

https://doi.org/10.1063/1.4883520

Google Scholar

Crossref

PubMed

82.

J. A.

de la Torre

,

P.

Español

, and

A.

Donev

, “

Finite element discretization of non-linear diffusion equations with thermal fluctuations

,”

J. Chem. Phys.

142

,

094115

(

2015

).

https://doi.org/10.1063/1.4913746

Google Scholar

Crossref

PubMed

83.

S.

Plyasunov

, “

On hybrid simulation schemes for stochastic reaction dynamics

,” arXiv: math/0504477.

84.

T.

Li

, “

Analysis of explicit tau-leaping schemes for simulating chemically reacting systems

,”

Multiscale Model. Simul.

6

,

417

–

436

(

2007

).

https://doi.org/10.1137/06066792x

Google Scholar

Crossref

85.

D. T.

Gillespie

, “

The chemical Langevin equation

,”

J. Chem. Phys.

113

,

297

–

306

(

2000

).

https://doi.org/10.1063/1.481811

Google Scholar

Crossref

86.

H. C.

Öttinger

and

M.

Grmela

, “

Dynamics and thermodynamics of complex fluids. II. Illustrations of a general formalism

,”

Phys. Rev. E

56

,

6633

(

1997

).

https://doi.org/10.1103/physreve.56.6633

Google Scholar

Crossref

87.

P.

Hänggi

,

H.

Grabert

,

P.

Talkner

, and

H.

Thomas

, “

Bistable systems: Master equation versus Fokker–Planck modeling

,”

Phys. Rev. A

29

,

371

–

378

(

1984

).

https://doi.org/10.1103/physreva.29.371

Google Scholar

Crossref

88.

J.

Keizer

,

Statistical Thermodynamics of Nonequilibrium Processes

(

Springer

,

Berlin

,

1987

).

Google Scholar

Crossref

89.

F. B.

Usabiaga

,

J. B.

Bell

,

R.

Delgado-Buscalioni

,

A.

Donev

,

T. G.

Fai

,

B. E.

Griffith

, and

G. S.

Peskin

, “

Staggered schemes for fluctuating hydrodynamics.

,”

Multiscale Model. Simul.

10

,

1369

–

1408

(

2012

).

https://doi.org/10.1137/120864520

Google Scholar

Crossref

90.

R.

Erban

and

S. J.

Chapman

, “

Stochastic modelling of reaction-diffusion processes: Algorithms for bimolecular reactions

,”

Phys. Biol.

6

,

046001

(

2009

).

https://doi.org/10.1088/1478-3975/6/4/046001

Google Scholar

Crossref

PubMed

91.

S. A.

Isaacson

, “

The reaction-diffusion master equation as an asymptotic approximation of diffusion to a small target

,”

SIAM J. Appl. Math.

70

,

77

–

111

(

2009

).

https://doi.org/10.1137/070705039

Google Scholar

Crossref

92.

D.

Bernstein

, “

Simulating mesoscopic reaction-diffusion systems using the Gillespie algorithm

,”

Phys. Rev. E

71

,

041103

(

2005

).

https://doi.org/10.1103/physreve.71.041103

Google Scholar

Crossref

93.

S. A.

Isaacson

and

C. S.

Peskin

, “

Incorporating diffusion in complex geometries into stochastic chemical kinetics simulations.

,”

SIAM J. Sci. Comput.

28

,

47

–

74

(

2006

).

https://doi.org/10.1137/040605060

Google Scholar

Crossref

94.

H.-W.

Kang

,

L.

Zheng

, and

H. G.

Othmer

, “

A new method for choosing the computational cell in stochastic reaction-diffusion systems

,”

J. Math. Biol.

65

,

1017

–

1099

(

2012

).

https://doi.org/10.1007/s00285-011-0469-6

Google Scholar

Crossref

PubMed

95.

D.

Fange

,

O. G.

Berg

,

P.

Sjöberg

, and

J.

Elf

, “

Stochastic reaction-diffusion kinetics in the microscopic limit

,”

Proc. Natl. Acad. Sci. U. S. A.

107

,

19820

–

19825

(

2010

).

https://doi.org/10.1073/pnas.1006565107

Google Scholar

Crossref

PubMed

96.

S.

Hellander

,

A.

Hellander

, and

L.

Petzold

, “

Reaction rates for mesoscopic reaction-diffusion kinetics

,”

Phys. Rev. E

91

,

023312

(

2015

).

https://doi.org/10.1103/physreve.91.023312

Google Scholar

Crossref

97.

S. A.

Isaacson

, “

A convergent reaction-diffusion master equation

,”

J. Chem. Phys.

139

,

054101

(

2013

).

https://doi.org/10.1063/1.4816377

Google Scholar

Crossref

PubMed

98.

Y.

Cao

and

R.

Erban

, “

Stochastic Turing patterns: Analysis of compartment-based approaches

,”

Bull. Math. Biol.

76

,

3051

–

3069

(

2014

).

https://doi.org/10.1007/s11538-014-0044-6

Google Scholar

Crossref

PubMed

99.

C. R.

Doering

,

K. V.

Sargsyan

, and

P.

Smereka

, “

A numerical method for some stochastic differential equations with multiplicative noise

,”

Phys. Lett. A

344

,

149

–

155

(

2005

).

https://doi.org/10.1016/j.physleta.2005.06.045

Google Scholar

Crossref

100.

E.

Moro

and

H.

Schurz

, “

Boundary preserving semianalytic numerical algorithms for stochastic differential equations

,”

SIAM J. Sci. Comput.

29

,

1525

–

1549

(

2007

).

https://doi.org/10.1137/05063725x

Google Scholar

Crossref

101.

A.

Abdulle

,

G.

Vilmart

, and

K. C.

Zygalakis

, “

Weak second order explicit stabilized methods for stiff stochastic differential equations.

,”

SIAM J. Sci. Comput.

35

,

A1792

–

A1814

(

2013

).

https://doi.org/10.1137/12088954x

Google Scholar

Crossref

102.

L.

Einkemmer

and

A.

Ostermann

, “

Overcoming order reduction in diffusion-reaction splitting. Part 1: Dirichlet boundary conditions.

,”

SIAM J. Sci. Comput.

37

,

A1577

–

A1592

(

2015

).

https://doi.org/10.1137/140994204

Google Scholar

Crossref

103.

W. L.

Briggs

,

V. E.

Henson

, and

S. F.

McCormick

,

A Multigrid Tutorial

, 2nd ed. (

SIAM

,

2000

).

Google Scholar

Crossref

104.

D. F.

Anderson

and

J. C.

Mattingly

, “

A weak trapezoidal method for a class of stochastic differential equations

,”

Commun. Math. Sci.

9

,

301

–

318

(

2011

).

https://doi.org/10.4310/cms.2011.v9.n1.a15

Google Scholar

Crossref

105.

C. A.

Rendleman

,

V. E.

Beckner

,

M.

Lijewski

,

W.

Crutchfield

, and

J. B.

Bell

, “

Parallelization of structured, hierarchical adaptive mesh refinement algorithms

,”

Comput. Visualization Sci.

3

,

147

–

157

(

2000

), software available at https://ccse.lbl.gov/BoxLib.

https://doi.org/10.1007/pl00013544

Google Scholar

Crossref

106.

R.

Ramadugu

,

S. P.

Thampi

,

R.

Adhikari

,

S.

Succi

, and

S.

Ansumali

, “

Lattice differential operators for computational physics

,”

EPL

101

,

50006

(

2013

).

https://doi.org/10.1209/0295-5075/101/50006

Google Scholar

Crossref

107.

Z.

Jdrzejewski-Szmek

and

K. T.

Blackwell

, “

Asynchronous

τ

-leaping

,”

J. Chem. Phys.

144

,

125104

(

2016

).

https://doi.org/10.1063/1.4944575

Google Scholar

Crossref

PubMed

108.

T.

Ando

and

J.

Skolnick

, “

Crowding and hydrodynamic interactions likely dominate in vivo macromolecular motion

,”

Proc. Natl. Acad. Sci. U. S. A.

107

,

18457

–

18462

(

2010

).

https://doi.org/10.1073/pnas.1011354107

Google Scholar

Crossref

PubMed

109.

J. A.

Dix

and

A. S.

Verkman

, “

Crowding effects on diffusion in solutions and cells

,”

Annu. Rev. Biophys.

37

,

247

–

263

(

2008

).

https://doi.org/10.1146/annurev.biophys.37.032807.125824

Google Scholar

Crossref

PubMed

110.

J.

Skolnick

, “

Perspective: On the importance of hydrodynamic interactions in the subcellular dynamics of macromolecules

,”

J. Chem. Phys.

145

,

100901

(

2016

).

https://doi.org/10.1063/1.4962258

Google Scholar

Crossref

PubMed

111.

A.

Donev

,

T. G.

Fai

, and

E.

Vanden-Eijnden

, “

A reversible mesoscopic model of diffusion in liquids: From giant fluctuations to Fick’s law

,”

J. Stat. Mech.: Theory Exp.

2014

,

P04004

.

https://doi.org/10.1088/1742-5468/2014/04/p04004

Crossref

112.

L.

Signon

,

B.

Nowakowski

, and

A. A.

Lemarchand

, “

Modeling somite scaling in small embryos in the framework of Turing patterns

,”

Phys. Rev. E

93

,

042402

(

2016

).

https://doi.org/10.1103/physreve.93.042402

Google Scholar

Crossref

PubMed

113.

T.

Biancalani

,

D.

Fanelli

, and

F.

Di Patti

, “

Stochastic Turing patterns in the brusselator model

,”

Phys. Rev. E

81

,

046215

(

2010

).

https://doi.org/10.1103/physreve.81.046215

Google Scholar

Crossref

114.

D.

Fanelli

,

C.

Cianci

, and

F.

Di Patti

, “

Turing instabilities in reaction-diffusion systems with cross diffusion.

,”

Eur. Phys. J. B

86

,

142

(

2013

).

https://doi.org/10.1140/epjb/e2013-30649-7

Google Scholar

Crossref

115.

N.

Kumar

and

W.

Horsthemke

, “

Effects of cross diffusion on Turing bifurcations in two-species reaction-transport systems

,”

Phys. Rev. E

83

,

036105

(

2011

).

https://doi.org/10.1103/physreve.83.036105

Google Scholar

Crossref

2017

Author(s)

Stochastic simulation of reaction-diffusion systems: A fluctuating-hydrodynamics approach

I. INTRODUCTION

II. BACKGROUND

A. FHD description

1. Diffusion

2. Reaction

B. Structure factor

C. Schlögl model

III. SPATIAL DISCRETIZATION

A. Diffusion-only case

B. Reaction-diffusion system

C. Maintaining non-negative densities

IV. TEMPORAL INTEGRATORS

A. Schemes

1. Explicit midpoint schemes

2. Implicit midpoint schemes

B. Structure factor analysis

V. NUMERICAL RESULTS

A. Schlögl model at thermodynamic equilibrium

1. Continuous-time FHD equation

2. Time integration errors

B. Turing-like pattern formation

C. Front propagation

VI. CONCLUSIONS

ACKNOWLEDGMENTS

APPENDIX A: RDME-BASED SPLIT SCHEME

1. Split scheme

2. Structure factor analysis

APPENDIX B: AVERAGING FUNCTION $\tilde{𝑛}$

APPENDIX C: LINEARIZED EQUATION ANALYSIS

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

Stochastic simulation of reaction-diffusion systems: A fluctuating-hydrodynamics approach

I. INTRODUCTION

II. BACKGROUND

A. FHD description

1. Diffusion

2. Reaction

B. Structure factor

C. Schlögl model

III. SPATIAL DISCRETIZATION

A. Diffusion-only case

B. Reaction-diffusion system

C. Maintaining non-negative densities

IV. TEMPORAL INTEGRATORS

A. Schemes

1. Explicit midpoint schemes

2. Implicit midpoint schemes

B. Structure factor analysis

V. NUMERICAL RESULTS

A. Schlögl model at thermodynamic equilibrium

1. Continuous-time FHD equation

2. Time integration errors

B. Turing-like pattern formation

C. Front propagation

VI. CONCLUSIONS

ACKNOWLEDGMENTS

APPENDIX A: RDME-BASED SPLIT SCHEME

1. Split scheme

2. Structure factor analysis

APPENDIX B: AVERAGING FUNCTION 𝑛∼

APPENDIX C: LINEARIZED EQUATION ANALYSIS

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only

APPENDIX B: AVERAGING FUNCTION $\tilde{𝑛}$