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.

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,1 reversal of direction of front propagation,2 violation of the law of mass action,3 long-time tails in kinetics,4 emergence of new steady states5 and patterns,6 acceleration of pattern formation,7 enhanced induction time for ignition,8 and the onset of homogeneous oscillations.9 Due to a small number of proteins involved in cellular functions,10 processes in cell biology are good examples11–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,17 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 splitting25 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,28 whereas the tau leaping method29 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.30 Several approaches have been proposed to improve the performance of stochastic sampling of the RDME, such as the next subvolume method35 and its parallel simulation version.37 Various implementations of the tau leaping method in a spatial context,38–40 and the time-dependent propensity for diffusion method,41 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,42 the multinomial simulation algorithm,43 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.48 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 Lifshitz49 and its validity has been justified for nonequilibrium systems50 through the theory of coarse graining.51 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 transport54,55 in binary fluid mixtures,56,57 multiphase flows,58 multispecies fluid mixtures,59,60 multispecies reactive mixtures,61 and electrolytes.62 

Compared to our previous work,61 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,61 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 approach63,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 equilibrium68 and the Model H equations for binary mixtures.69 

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 RDME42–47 are based on operator splitting using first-order Lie or second-order Strang splitting.70 In  Appendix A we review and discuss in more detail a split scheme that uses multinomial diffusion sampling71 for diffusion and SSA for reactions. These RDME-based schemes use a time step size Δt comparable to the hopping time scale τh=Δx2/(2dD) 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.7 In Section VI, we offer some concluding remarks and suggest future research directions.

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)=(Z1(𝒙,t),,Zd(𝒙,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: Zj(𝒙,t)Zj(𝒙,t)=δjjδ(𝒙𝒙)δ(tt) and Z(𝒙,t)Z(𝒙,t) = δ(𝒙𝒙)δ(tt). Similarly, we denote GWN vector and scalar random processes by W(t) and W(t), respectively, and assume Wj(t)Wj(t)=δjjδ(tt) and W(t)W(t)=δ(tt).

We consider a reaction-diffusion system having Ns species undergoing Nr reactions in d-dimensional space. By denoting the number density of species s by ns(𝒙,t), the equations of FHD for 𝒏(𝒙,t)=(n1(𝒙,t),,nNs(𝒙,t)) are written formally as the SPDEs61 

tns=(Dsns+2Dsns(s(D))+r=1Nrνsr(ar(𝒏)+ar(𝒏)Zr(R)),
(1)

where Ds is the diffusion coefficient of species s, ar(𝒏) is the propensity function indicating the rate of reaction r, and νsr 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),

tns=D2sns+r=1Nrνsrar(𝒏).
(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

tns=(Dsns+2Dsnss(D))
(3)

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

𝒙˙s,i=2DsWs,i.
(4)

In this formal derivation,79 one defines the instantaneous number density field ns(𝒙,t)=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,80 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.32 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.81 

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.82 

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:30 

dns=ns(t+dt)ns(t)=1ΔVr=1NrνsrP(ar(𝒏)ΔVdt),
(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 ar(𝒏) that we use in this work is described in Section III C. Henceforth, we will formally write Eq. (5) in the differential form,

ddtns=r=1NrνsrP(ar(𝒏)ΔVdt)ΔVdt.
(6)

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

The chemical Langevin equation (CLE)85 is obtained under the assumption that the mean number of reaction occurrences is large.30 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,

ddtns=r=1Nrνsr[ar(𝒏)+ar(𝒏)ΔVWr].
(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 work61 was that the Langevin description (7) is not consistent with equilibrium statistical mechanics. Alternative formulations based on a Langevin diffusion description86,87 that are consistent at thermodynamic equilibrium fail to correctly model relaxation toward equilibrium.61 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.

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

Ss(𝒌)=Vδn^s,𝒌δn^s,𝒌*,
(8)

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

n^s,𝒌(t)=1Vns(𝒙,t)ei𝒌𝒙d𝒙.
(9)

Here we have assumed a periodic domain of volume V and defined δn^s,𝒌=n^s,𝒌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, Ns = 1, and drop the subscript s for species, to write Eq. (1) as

tn=D2n+(2Dn(D))+a(n)+2Γ(n)(R),
(10)

where

a(n)=r=1Nrνrar(n),Γ(n)=12r=1Nrνr2ar(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 n¯n, where a(n¯)=0 and a(n¯)<0. The linearization of Eq. (10) around this equilibrium state is given by the central limit theorem,

tn=D2n+2Dn¯(D)r(nn¯)+2Γ¯(R),
(12)

where r=a(n¯)>0 is the effective reaction rate near equilibrium and Γ¯=Γ(n¯).

The Fourier transform of Eq. (12) gives

ddtδn^𝒌=Dk2δn^𝒌+2Dn¯i𝒌^𝒌(D)rδn^𝒌+2Γ¯^𝒌(R).
(13)

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

S(𝒌)=Dn¯k2+Γ¯Dk2+r=n¯k2+Γ¯/Dk2+2,
(14)

where =D/r denotes the penetration depth. From Eq. (14), we observe that there are two limiting cases. In the small wave number limit k1, S(𝒌) becomes Γ¯/r and does not depend on diffusion. In fact, the result S(𝟎)=Γ¯/r is also obtained from the CME assuming the whole system is well-mixed. On the other hand, in the large wave number limit k1, S(𝒌) becomes 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 Γ¯=n¯r and S(𝒌)=n¯, consistent with a product Poisson distribution with mean number density 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.88 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.

The Schlögl model75,76 is given by the chemical reactions for species X,

2Xk1k23X,k3k4X.
(15)

Hence, we have Ns=1, Nr=4, ν1=ν3=1, ν2=ν4=1, a(n)=k1n2k2n3+k3k4n, and Γ(n)=12(k1n2+k2n3+k3+k4n). 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, k1neq2=k2neq3 and k3= k4neq, the system is in thermodynamic equilibrium and the distribution follows Poisson statistics with mean number density neq. 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).76 

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(𝒌)=neq 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 k1. 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.

In this section, we discuss spatial discretization of the FHD equation using a finite-volume approach63,64 that converts the SPDE into stochastic ordinary differential equations (SODEs) for the cell number density ns,𝒊(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.

Due to the lack of regularity of (s(D)(𝒙,t) in Eq. (3), pointwise values of ns(𝒙,t) are not physically meaningful. Hence, we consider instead the spatial average of ns(𝒙,t) over a cell. We partition the system domain L1×L2××Ld into cells of volume ΔV=Δx1Δxd and denote the cell number density of species s in cell 𝒊=(i1,,id) as

ns,𝒊(t)=1ΔVcell𝒊ns(𝒙,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 𝒊+12𝒆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 D2sns and introduce a staggered grid for the stochastic diffusive flux term (2Dsns(s(D)), see Fig. 1. For d = 1, a formally second-order spatial discretization of Eq. (3) is written as

ddtns,i=Dsns,i+12ns,i+ns,i1Δx2+2DsΔVns,i+12Ws,i+12ns,i12Ws,i12Δx.
(17)

The spatial average of ns(𝒙,t) over the interval of length Δx around face i±12 is approximated by ns,i±12(t), whereas that of s(𝒙,t) is modeled by 1ΔVWs,i±12(t). To close the equation, ns,i±12(t) is approximated by an average of ns,i(t) and ns,i±1(t), that is, ns,i±12=n(ns,i,ns,i±1). Natural candidates for the averaging function n(n1,n2) would be the Pythagorean means: the arithmetic, geometric, and harmonic means. We choose a modified arithmetic average for n(n1,n2) described in Section III C, for reasons detailed in  Appendix B.

FIG. 1.

Finite-volume spatial discretization in two dimensions. The cell-averaged number density ns,𝒊(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.

FIG. 1.

Finite-volume spatial discretization in two dimensions. The cell-averaged number density ns,𝒊(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.

Close modal

Generalization of the spatial discretization (17) to higher dimensions is straightforward. For each face, a GWN process Ws,𝒇 is defined and ns,𝒇(t) is calculated from the cell number densities of the two cells sharing the face by using the averaging function n(n1,n2), see Fig. 1. By introducing notations ns(t){ns,𝒊(t)}, Ws(t){Ws,𝒇(t)}, and ns(t)n(ns(t)), we express the resulting SODEs for {ns,𝒊(t)} as

ddtns=Dd2sns+2DsΔVd(nsWs),
(18)

with the understanding that d2 denotes the standard (2d + 1)-point discrete Laplacian operator, and d denotes a discrete divergence operator.

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

ddtns=Dd2sns+2DsΔVd(nsWs)+r=1NrνsrP(ar(𝒏)ΔVdt)ΔVdt.
(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.61 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 arguments90,92,93 and mathematical analysis.94 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 Smoluchowski17 model or the Doi model,90 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 developed97 and could be combined with our FHD description of diffusion. The dependence of stochastic Turing patterns on the grid size has been also investigated.98 

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., ns,𝒊Δ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 n and the propensity functions ar that appear in the multiplicative noise terms. Hence, we carefully modify the form of n and ar 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:

n(n1,n2)=n1+n22H0(n1ΔV)H0(n2ΔV),
(20)

where

H0(x)={0(x0)x(0<x<1)1(x1)
(21)

is a smoothed Heaviside function. The smoothed Heaviside function H0 is introduced to ensure the continuity of n at n1 = 0 or n2 = 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 n10 or n20. 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., niΔV1), 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 niΔV1. 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 ar(𝒏), 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 ns2 (or ns3, ), then replace it by ns(ns1ΔV) (or ns(ns1ΔV)(ns2ΔV), ). 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 nsΔV follows Poisson statistics with mean n¯sΔV, ns(ns1ΔV)=n¯s2 and ns(ns1ΔV)(ns2ΔV)=n¯s3.

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 ar(𝒏) 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)=k1n+(n1ΔV)+k2n+(n1ΔV)+(n2ΔV)++k3k4n+,
(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.

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,101 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., 2dDΔtΔx2). 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 conditions102 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 Δx1==Δxd=Δx, we define

α=rΔt,β=DΔtΔx2,
(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 =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,

β12dΔtΔtmaxΔx22dD.
(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.

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

nsk+1=nsk+DsΔtd2nsk+2DsΔtΔVd(nskWsk)+r=1NrνsrP(arkΔVΔt)ΔV,
(25)

where superscripts denote the point in time at which quantities are evaluated, e.g., nsk=ns(kΔt) and ark=ar(𝒏(kΔt)), and we have used the compact notation for spatial discretization introduced in Section III B. Here, kΔt(k+1)ΔtW(t)dt has been replaced by ΔtWk, where Wk 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 Ws,𝒇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

nsk+1=nsk+DsΔtd2nsk+2DsΔtΔVd(nskWsk)+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 nsk+1. However, a by a stability condition (for derivation, see Eq. (C7))

β+α4d12d,
(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 nsk+1,

nsk+1=nsk+DsΔtd2(nsk+nsk+12)+2DsΔtΔVd(nskWsk)+r=1NrνsrP(arkΔVΔt)ΔV.
(28)

The linear system (28) can be solved efficiently iteratively using multigrid relaxation;103 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 CLE104 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,

ns=nsk+DsΔt2d2nsk+DsΔtΔVd(nskWs(1))+r=1NrνP(1)sr(arkΔVΔt/2)ΔV,
(29a)
nsk+1=nsk+DsΔtd2ns+DsΔtΔVd(nskWs(1))+DsΔtΔVd(nsWs(2))+r=1NrνP(1)sr(arkΔVΔt/2)ΔV+r=1NrνP(2)sr((2arark)+Δ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(1) 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 (2arark)+, 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 ns:

ns=n(nsk),
(30a)
ns=n(ns),
(30b)
ns=n((2nsnsk)+).`
(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.104 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,

ns=nsk+DsΔt2d2nsk+DsΔtΔVd(nskWs(1)),
(31a)
ns=ns+s(1)(𝒏,Δt2),
(31b)
nsk+1=nsk+DsΔtd2ns+DsΔtΔVd(nskWs(1))+DsΔtΔVd(nsWs(2))+s(1)(𝒏,Δt2)+s(2)(𝒏,Δt2).
(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,

ns=nsk+DsΔt2d2ns+DsΔtΔVd(nskWs(1))+r=1NrνP(1)sr(arkΔVΔt/2)ΔV,
(32a)
nsk+1=nsk+DsΔtd2(nsk+nsk+12)+DsΔtΔVd(nskWs(1))+DsΔtΔVd(nsWs(2))+r=1NrνP(1)sr(arkΔVΔt/2)ΔV+r=1NrνP(2)sr((2arark)+ΔVΔt/2)ΔV,
(32b)

where the three options for ns are given in Eqs. (30). When SSA is used for the reactions, we obtain the implicit midpoint SSA (ImMidSSA) scheme,

ns=nsk+DsΔt2d2ns+DsΔtΔVd(nskWs(1)),
(33a)
nsk+1=nsk+DsΔtd2(nsk+nsk+12)+DsΔtΔVd(nskWs(1))+DsΔtΔVd(nsWs(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.64 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 α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 CLE61 is that here we replaced the GWN in the chemical noise with Poisson noise and used a weakly second-order tau leaping method65,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.

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Δx2/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 −k2 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 k defined by

k=sin(kΔx2)Δx2.
(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 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 Δtmax for the explicit schemes, and it is increased by a factor of 20 (Δt=10Δtmax) for the implicit schemes. As described below, the accuracy at diffusion-dominated scales k1 and reaction-dominated scales k1 largely depends on how diffusion (i.e., explicit or implicit) and reaction (i.e., tau leaping or SSA) are treated, respectively.

FIG. 2.

Discrete structure factors S(k) for the ExMidTau, ImMidTau, ExMidSSA, and ImMidSSA schemes for the one-dimensional linearized FHD equation (12) with Γ¯/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.

FIG. 2.

Discrete structure factors S(k) for the ExMidTau, ImMidTau, ExMidSSA, and ImMidSSA schemes for the one-dimensional linearized FHD equation (12) with Γ¯/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.

Close modal

As the time step size approaches the diffusive stability limit, β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 k1, for which the ImMidTau scheme shows some errors because of the relatively large value of α=0.5. However, for intermediate scales k1, the ImMidTau scheme is more accurate because it attains third-order accuracy for static covariances.

At diffusion-dominated scales k1, 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.64 In addition, the reaction is incorporated in such a way that maintains fluctuation-dissipation balance for k1 even for relatively large values of α. For small α, the time integration error of the ImMidTau scheme for the structure factor at the maximum wavenumber kmax=πβ/α is estimated as

S(kmax)S0(kmax)S0(kmax)β2(1+2β)2α2,
(35)

where S0(k)=limΔt0S(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 k1 for large β. For the ImMidSSA scheme, β in the numerator is replaced by β+(1+2β)(Γ¯n¯r1) and a similar stable behavior for large β is observed. By expanding S(kmax) − S0(kmax) for small Δt, we also see that the error is O(α2β)=O(Δt3) for the ImMidTau scheme, whereas it is O(Δt2) for the ImMidSSA except at thermodynamic equilibrium (i.e., except when Γ¯=n¯r), where it is third-order accurate.

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.7 

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 sampling71 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 Δt0. 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Δx1Δxd. 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 ns,𝒊Δ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.

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) = neq is constant, both with and without (i.e., diffusion only) chemical reactions. Specifically, we set the rate constants as k1 = k2 = k3 = k4 = 0.1 (see Eq. (15)), which gives neq = 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=103; 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 Nc = 512 grid cells in one dimension for the arithmetic, geometric, and harmonic mean (AM, GM, HM) averaging functions n(n1,n2). 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 k1. 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.

Static fluctuations of the spatially discretized FHD equation (19) at thermodynamic equilibrium. (Left) Structure factors S(k) calculated by using different averaging functions n. The solid lines depict the results from the equilibrium Schlögl model (neqΔV=10), whereas the dotted lines are from the corresponding diffusion-only system. Note that the exact result for both systems is S(k)=neq independent of 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 PPoisson(N) and its Gaussian approximation PGauss(N). The inset shows the errors P(N) − PPoisson(N) with error bars corresponding to two standard deviations.

FIG. 3.

Static fluctuations of the spatially discretized FHD equation (19) at thermodynamic equilibrium. (Left) Structure factors S(k) calculated by using different averaging functions n. The solid lines depict the results from the equilibrium Schlögl model (neqΔV=10), whereas the dotted lines are from the corresponding diffusion-only system. Note that the exact result for both systems is S(k)=neq independent of 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 PPoisson(N) and its Gaussian approximation PGauss(N). The inset shows the errors P(N) − PPoisson(N) with error bars corresponding to two standard deviations.

Close modal

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)=(N12)/ΔV(N+12)/ΔVρ(n)dn,
(36)

and compare it with a Poisson distribution PPoisson(N) with mean neqΔV, as well as a Gaussian distribution PGauss(N) having the same mean and variance as PPoisson(N). The agreement between P(N) and PPoisson(N) is remarkable in the sense that FHD was originally proposed to account for only second moments of (small) Gaussian fluctuations. Since PPoisson(N) is significantly different from PGauss(N) for neqΔ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 PPoisson(N) by using the following measures. First, we compute the Kullback–Leibler (KL) divergence (distance),

DKL=N=0PPoisson(N)logPPoisson(N)P(N).
(37)

Second, we compute the probability of negative number densities,

Pneg=0ρ(n)dn.
(38)

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

ζ=Cov[n𝒊,n𝒊±𝒆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 ns in Eqs. (30). For ζ, the three options give similar values within standard errors of estimation. For DKL and Pneg, 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 Nc = 64 cells in one dimension for the different schemes. As expected, for small values of Δt, all schemes give similar values. DKL converges to a small value, which is consistent with the good agreement between P(N) and PPoisson(N) seen in the right panel of Fig. 3. Pneg is observed to converge to zero as Δt0, 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 Δtmax, rapid worsening is observed in both explicit schemes in all three measures. While Pneg and ζ behave similarly for both schemes, DKL remains small for larger values of Δt for the ExMidSSA scheme compared to the ExMidTau scheme.

FIG. 4.

Deviations from the correct equilibrium distribution PPoisson(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 Δtmax of the explicit schemes, see Eq. (24)The error bars correspond to two standard deviations.

FIG. 4.

Deviations from the correct equilibrium distribution PPoisson(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 Δtmax of the explicit schemes, see Eq. (24)The error bars correspond to two standard deviations.

Close modal

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 DKL 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 DKL. For Pneg 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 Δtmax1/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,103 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.

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

U+Wk1V+W,2Vk2k3W,Uk4k5,Vk6k7.
(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.61 Specifically, we set k1=k2=2×104, k3 = 1, k4=3.33×103, k5 = 16.7, k6=3.67×102, k7 = 4.44, and DU = 0.1, DV = DW = 0.01. We note that on the limit cycle, number densities of the three species oscillate in significantly different ranges: nU(999,2024), nV(302,645), and nW(18.2,83.2). For a physical domain with side lengths Lx = Ly = 32, we use three spatial resolutions with grid sizes Nc = 642, 1282, and 2562 cells. For the initial number densities, we choose a point on the limit cycle (nU0,nV0,nW0)=(1686,534,56.4) and generate the initial number of molecules of each species s in each cell from a Poisson distribution with mean ns0ΔV=ns0AΔ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 642 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×104, 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 A100 due to the very large number of U molecules (as many as 2×106A) 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×104 (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.

Two types of steady-state Turing-like patterns observed in the two-dimensional BPM model with 642 cells and A = 1. Snapshots of nU are obtained at t=5×104 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).

FIG. 5.

Two types of steady-state Turing-like patterns observed in the two-dimensional BPM model with 642 cells and A = 1. Snapshots of nU are obtained at t=5×104 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).

Close modal

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 n¯U(t)=1Nc𝒊nU,𝒊(t). In the left panel, we compare sample trajectories of n¯U(t) for A = 1 and A = 10 for RDME (similar results are obtained for FHD) and for deterministic reaction-diffusion. While n¯U(t) initially oscillates as in the limit cycle, as the Turing-like pattern begins to form, the oscillation amplitude decays and 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 work61 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.

Turing-like pattern formation in the two-dimensional BPM model with 642 cells. (Left) Spatially averaged density 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 nU for A = 10 at four different times t at which n¯U(t) attains a local minimum, indicated by circles in the left panel.

FIG. 6.

Turing-like pattern formation in the two-dimensional BPM model with 642 cells. (Left) Spatially averaged density 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 nU for A = 10 at four different times t at which n¯U(t) attains a local minimum, indicated by circles in the left panel.

Close modal

We generate 16 sample trajectories for each set of parameter values, fit each realization of n¯U(t) using seven fitting parameters a1,,a7 to

n¯U(t)=(1tanhta1a2)(a3sin(a4t+a5)+a6)+a7,
(41)

and compare the distributions of the fitting parameters. Note that a1 and a7 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 (a1, a7) 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 2562 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., a1 decreases), as already seen in Fig. 6, while the steady spatial average density a7 becomes smaller. In addition, while the variance of a1 does not change significantly as A varies, the variance of a7 becomes larger for smaller A.

FIG. 7.

(Left) Scatter plots of the decay onset time a1 and the steady spatial average density a7 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 a¯7 of a7 over 16 samples for the ImMidTau and ImMidSSA schemes with 642 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 ai are normalized by the average values a¯idet for deterministic reaction-diffusion started with random initial conditions corresponding to A = 10.

FIG. 7.

(Left) Scatter plots of the decay onset time a1 and the steady spatial average density a7 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 a¯7 of a7 over 16 samples for the ImMidTau and ImMidSSA schemes with 642 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 ai are normalized by the average values a¯idet for deterministic reaction-diffusion started with random initial conditions corresponding to A = 10.

Close modal

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 Δtmax(R)1.3. The right panel of Fig. 7 shows the mean values a¯7 of the steady spatial average density a7 over 16 samples versus Δt for 642 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 a¯7 rapidly deviates from the RDME result as Δt approaches Δtmax(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 Δtmax(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.

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

Ak1,2A+Bk23A,Bk3k4.
(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 k1 = 0.4, k2 = 0.137, k3 = 0.1, k4 = 1 and DA = 1, DB = 10. For a physical domain with side lengths Lx = Ly = Lz = 512, we use 2563 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,

ns,𝒊0=ns(1)+12(1+tanhr𝒊Rξ)(ns(2)ns(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, (nA(1),nB(1))=(2.16,1.35), while the rest of the domain is initialized in the second uniform equilibrium state, (nA(2),nB(2))=(0,10), with a smooth transition region between the two states of width 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 ns,𝒊0ΔV=ns,𝒊0AΔ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.105 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 1012 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.7 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 isotropic106 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.

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.

FIG. 8.

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.

Close modal

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 ns,𝒊(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 1012 molecules using the BoxLib CFD software framework.105 

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 leaping107 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).110 Excluded volume due to steric repulsion introduces cross-diffusion effects, i.e., coupling between the diffusive fluxes for different species,112 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.61 In future work, we will investigate these hydrodynamic effects on stochastic reaction-diffusion phenomena.

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.

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 sampling43,71 for diffusion and SSA for reaction via Strang splitting.70 

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 Tartakovsky71 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.43 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 Ns,𝒊𝒊, the change in the number density ns,𝒊 can be expressed in terms of the sum of the inward and outward fluxes,

ns,𝒊(t+Δt)=ns,𝒊(t)+1ΔVj=1d(Ns,𝒊+𝒆j𝒊+Ns,𝒊𝒆j𝒊Ns,𝒊𝒊+𝒆jNs,𝒊𝒊𝒆j)
(A1a)
=ns,𝒊(t)+𝔇𝒊(ns(t),Δt),
(A1b)

where ns(t)={ns,𝒊(t)}. For each cell i, the outward fluxes (Ns,𝒊𝒊+𝒆1,Ns,𝒊𝒊𝒆1,,Ns,𝒊𝒊+𝒆d,Ns,𝒊𝒊𝒆d,Ns,𝒊𝒊) are random variables sampled from the multinomial distribution with 𝒊Ns,𝒊𝒊=ns,𝒊(t)ΔV total trials and probabilities (ps,ps,,ps,ps,12dps), where ps=DsΔt/Δx2, where we have assumed Δx1==Δxd=Δ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 12dps0. We also note that the number of molecules in a cell never becomes negative due to the constraint 𝒊Ns,𝒊𝒊=ns,𝒊(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

ns,𝒊(t+Δt)=ns,𝒊(t)+s(𝒏𝒊(t),Δt).
(A2)

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

ns,𝐢=ns,𝐢k+s(𝐧𝐢k,Δt/2),
(A3a)
ns,𝐢=ns,𝐢+𝔇𝒊(ns,Δt),
(A3b)
ns,𝐢k+1=ns,𝐢+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 Δt0, 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 possible43,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 Δtmax (equivalently, to β=0.5), for α=0.1β, see Eq. (23). While S(k) is accurately reproduced at the reaction-dominated scales k1 for all values of β, it becomes inaccurate for smaller scales as Δt approaches Δtmax. We recall that for a system at thermodynamic equilibrium, the split scheme exactly preserves the correct equilibrium distribution for any Δt<Δtmax. 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.

Structure factors S(k) obtained from the SSA/2+MN+SSA/2 scheme (A3) for the one-dimensional linearized FHD equation (12) with Γ¯/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 k is used, see Eq. (34)

FIG. 9.

Structure factors S(k) obtained from the SSA/2+MN+SSA/2 scheme (A3) for the one-dimensional linearized FHD equation (12) with Γ¯/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 k is used, see Eq. (34)

Close modal

In this appendix, we show that the arithmetic mean should be chosen as the averaging function n(n1,n2). 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 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 n(n1,n2) for small values of n1 and n2. Here we focus on the continuous-time case and do not assume any specific temporal integrator.

The corresponding microscopic system has Ncn¯ΔV molecules, where n¯ is the mean number density and Nc is the number of cells. The equilibrium distribution of the numbers of molecules in each cell follows the multinomial distribution with equal probabilities 1/Nc. Hence, it is straightforward to obtain the following second-order statistics of cell number density n𝒊:

Var[n𝒊]=Nc1Ncn¯ΔV,
(B1)
Cov[n𝒊1,n𝒊2]=1Ncn¯ΔVfor𝒊1𝒊2.
(B2)

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

S(𝒌)=n¯fornonzero𝒌.
(B3)

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

Var[n𝒊]=Nc1NcnΔV,
(B4)
Cov[n𝒊1,n𝒊2]=1NcnΔVfor𝒊1𝒊2,
(B5)
S(𝒌)=nfornonzero𝒌.
(B6)

Comparing Eqs. (B4)–(B6) to the correct result Eqs. (B1)–(B3) suggests that one needs to choose n so that n is as close as possible to n¯. It is easy to see that the arithmetic mean (AM) would give the right answer: nAM=12(n1+n2)=n¯. On the other hand, to calculate n for the geometric mean nGM=n1n2 or the harmonic mean nHM=2/(n11+n21), one needs to know the equilibrium distribution ρ(n1,n2) of two neighboring cells. However, under the reasonable assumption that all three averaging functions give similar distributions ρ(n1,n2) allowing only non-negative number densities, it can be easily shown that nHMnGMnAM 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 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 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 n(n1,n2) is a non-negative function that satisfies

n(n1,n2)=0forn10orn20,
(B7)

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

n(n1,n2)={12(n1+n2)ifn1>0and𝑛2>0,0otherwise,
(B8)

does not allow negative density (i.e., ρ(n)=0 for n<0). However, due to the discontinuity of n at n1 = 0 or n2 = 0, ρ(n) does not decrease to zero as n becomes zero (i.e., limn0+ρ(n)>0) and a delta function is formed at n = 0, see Fig. 10.

FIG. 10.

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 H0The results are obtained from the one-dimensional diffusion-only FHD system (18) having D=n¯=Δx=1, A = 5, and Nc = 512 by using the EMTau scheme (25) with Δt=103

FIG. 10.

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 H0The results are obtained from the one-dimensional diffusion-only FHD system (18) having D=n¯=Δx=1, A = 5, and Nc = 512 by using the EMTau scheme (25) with Δt=103

Close modal

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 n1(t) has a small positive value, n can be large if the density in the neighboring cell n2(t) is large. In this case, n1(t+Δt) can become negative due to the stochastic diffusive flux. But, once n1 has become negative, the stochastic diffusive flux is turned off, and due to the deterministic diffusion, n1 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 Δt1/2 and Δt1/2, respectively, and the peak becomes a delta function in the limit Δt0.

We note that the averaging in Eq. (B8) can be expressed as n=12(n1+n2)H(n1ΔV)H(n2ΔV), where H is the Heaviside function. To avoid a discontinuity in 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 0N1 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, n¯ΔV=5. With the use of a smoothed Heaviside function, the spurious delta function at n = 0 as Δt0 is removed, and the probability of negative density is greatly reduced for small Δt.

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 δn^k(t),

ddtδn^k=Dk2δn^k+2Dn¯k2VWk(D)rδn^k+2Γ¯VWk(R),
(C1)

where the modified wavenumber k is defined in Eq. (34), and Wk(D)(t) and Wk(R)(t) are independent standard GWN processes. Compared to the continuous-space case (see Eq. (13)), k is replaced by 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

δn^k(t+Δt)=[1(Dk2+r)Δt]δn^k(t)+2Dn¯k2ΔtVW1+2Γ¯ΔtVW2,
(C2)

where W1 and W2 are independent standard normal random variables. We write this as

δn^k(t+Δt)=dMkδn^k(t)+NkW,
(C3)

where =d denotes being equal in distribution. For any given temporal integrator, we can straightforwardly obtain analytic expressions for Mk and NkNk. For example, for the EMTau scheme,

Mk=1(Dk2+r)Δt,
(C4a)
NkNk=2V(Dn¯k2+Γ¯)Δt.
(C4b)

The covariance of the noise NkNk 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)=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))

δn^k(t+Δt)=[Dk2Δt+erΔt]δn^k(t)+2Dn¯k2ΔtVW1+Γ¯(1e2rΔt)rVW2.
(C5)

While the expressions of Mk and NkNk* become complicated for the predictor-corrector midpoint schemes, a theoretical analysis is still tractable with the help of symbolic algebra tools.

By calculating Mk and NkNk* for a given numerical scheme, the stability condition and the structure factor can be obtained as follows. From the condition that the amplification factor Mk should satisfy

|Mk|1forallk,
(C6)

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

DΔtΔx2+rΔt412.
(C7)

The analytic expression for S(k)=Vδn^kδn^k* can be calculated from

S(k)=VNkNk*1Mk2,
(C8)

which is obtained from the time invariance relation δn^k(t)δn^k*(t)=δn^k(t+Δt)δn^k*(t+Δt) and Eq. (C3).63 From the analytical expressions of S(k) and S0(k)=limΔt0S(k), [S(kmax) − S0(kmax)]/S0(kmax) 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).

1.
D. A.
Kessler
and
H.
Levine
, “
Fluctuation-induced diffusive instabilities
,”
Nature
394
,
556
558
(
1998
).
2.
E.
Khain
,
Y. T.
Lin
, and
L. M.
Sander
, “
Fluctuations and stability in front propagation
,”
EPL
93
,
28001
(
2011
).
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
).
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
).
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
).
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
).
7.
A.
Lemarchand
and
B.
Nowakowski
, “
Do the internal fluctuations blur or enhance axial segmentation?
,”
EPL
94
,
48004
(
2011
).
8.
P.
Peeters
,
F.
Baras
, and
G.
Nicolis
, “
Sensitivity of the chlorite-thiosulfate system: A stochastic approach
,”
J. Chem. Phys.
93
,
7321
7328
(
1990
).
9.
M.
Malek Mansour
,
J.
Dethier
, and
F.
Baras
, “
Onset of homogeneous oscillations in reactive systems
,”
J. Chem. Phys.
114
,
9265
9275
(
2001
).
10.
N.
Fedoroff
and
W.
Fontana
, “
Small numbers of big molecules
,”
Science
297
,
1129
1131
(
2002
).
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
).
12.
D.
Fange
and
J.
Elf
, “
Noise-induced min phenotypes in E. coli
,”
PLoS Comput. Biol.
2
,
e80
(
2006
).
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
).
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
).
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
.
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
).
17.
M.
von Smoluchowski
, “
Versuch einer mathematischen theorie der koagulationskinetik kolloider lösungen
,”
Z. Phys. Chem.
92
,
129
168
(
1917
).
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
).
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
).
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
).
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
).
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
).
23.
P.
Dziekan
,
A.
Lemarchand
, and
B.
Nowakowski
, “
Particle dynamics simulations of Turing patterns
,”
J. Chem. Phys.
137
,
074107
(
2012
).
24.
P.
Dziekan
,
J. S.
Hansen
, and
B.
Nowakowski
, “
Nanoscale Turing structures
,”
J. Chem. Phys.
141
,
124106
(
2014
).
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
).
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
).
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
).
28.
D. T.
Gillespie
, “
A general method for numerically simulating the stochastic time evolution of coupled chemical reactions
,”
J. Comput. Phys.
22
,
403
(
1976
).
29.
D. T.
Gillespie
, “
Approximate accelerated stochastic simulation of chemically reacting systems
,”
J. Chem. Phys.
115
,
1716
(
2001
).
30.
D. T.
Gillespie
,
A.
Hellander
, and
L. R.
Petzold
, “
Perspective: Stochastic algorithms for chemical kinetics
,”
J. Chem. Phys.
138
,
170901
(
2013
).
31.
T.
Székeley
, Jr.
and
K.
Burrage
, “
Stochastic simulation in systems biology
,”
Comput. Struct. Biotechnol. J.
12
,
14
25
(
2014
).
32.
C. W.
Gardiner
,
Handbook of Stochastic Methods
, 3rd ed. (
Springer
,
New York
,
2004
).
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
).
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
).
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
).
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
).
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
).
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
).
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
).
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
).
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
).
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
).
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
).
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
).
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
).
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
).
48.
H. C.
Öttinger
,
Beyond Equilibrium Thermodynamics
(
John Wiley & Sons
,
2005
).
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
.
50.
P.
Español
, “
Stochastic differential equations for non-linear hydrodynamics
,”
Physica A
248
,
77
96
(
1998
).
51.
P.
Español
,
J. G.
Anero
, and
I.
Zúñga
, “
Microscopic derivation of discrete hydrodynamics
,”
J. Chem. Phys.
131
,
244117
(
2009
).
52.
N. K.
Voulgarakis
and
J.-W.
Chu
, “
Bridging fluctuating hydrodynamics and molecular dynamics simulations of fluids
,”
J. Chem. Phys.
130
,
134111
(
2009
).
53.
J. M.
Ortiz de Zárate
and
J. V.
Sengers
,
Hydrodynamic Fluctuations in Fluids and Fluid Mixtures
(
Elsevier
,
2006
).
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
).
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
.
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
).
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
).
58.
A.
Chaudhri
,
J. B.
Bell
,
A. L.
Garica
, and
A.
Donev
, “
Modeling multiphase flow using fluctuating hydrodynamics
,”
Phys. Rev. E
90
,
033014
(
2014
).
59.
K.
Balakrishnan
,
A. L.
Garcia
,
A.
Donev
, and
J. B.
Bell
, “
Fluctuating hydrodynamics of multispecies nonreactive mixtures
,”
Phys. Rev. E
89
,
013017
(
2014
).
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
).
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
).
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
).
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
).
64.
S.
Delong
,
B. E.
Griffith
,
E.
Vanden-Eijnden
, and
A.
Donev
, “
Temporal integrators for fluctuating hydrodynamics
,”
Phys. Rev. E
87
,
033302
(
2013
).
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
).
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
).
67.
P. J.
Atzberger
, “
Spatially adaptive stochastic numerical methods for intrinsic fluctuations in reaction-diffusion systems
,”
J. Comput. Phys.
229
,
3474
3501
(
2010
).
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
).
69.
S. P.
Thampi
,
I.
Pagonabarraga
, and
R.
Adhikari
, “
Lattice-Boltzmann-Langevin simulations of binary mixtures
,”
Phys. Rev. E
84
,
046709
(
2011
).
70.
G.
Strang
, “
On the construction and comparison of difference schemes
,”
SIAM J. Numer. Anal.
5
,
506
517
(
1968
).
71.
A.
Balter
and
A. M.
Tartakovsky
, “
Multinomial diffusion equation
,”
Phys. Rev. E
83
,
061143
(
2011
).
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
).
73.
M. A.
Karzazi
,
A.
Lemarchand
, and
M.
Mareschal
, “
Fluctuation effects on chemical wave fronts
,”
Phys. Rev. E
54
,
4888
4975
(
1996
).
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
).
75.
F.
Schlögl
, “
Chemical reaction models for non-equilibrium phase transitions
,”
Z. Phys.
253
,
147
161
(
1972
).
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
).
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
).
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
).
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
).
80.
G. D. C.
Kuiken
,
Thermodynamics of Irreversible Processes: Applications to Diffusion and Rheology
(
Wiley
,
1994
).
81.
A.
Donev
and
E.
Vanden-Eijnden
, “
Dynamic density functional theory with hydrodynamic interactions and fluctuations
,”
J. Chem. Phys.
140
,
234115
(
2014
).
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
).
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
).
85.
D. T.
Gillespie
, “
The chemical Langevin equation
,”
J. Chem. Phys.
113
,
297
306
(
2000
).
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
).
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
).
88.
J.
Keizer
,
Statistical Thermodynamics of Nonequilibrium Processes
(
Springer
,
Berlin
,
1987
).
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
).
90.
R.
Erban
and
S. J.
Chapman
, “
Stochastic modelling of reaction-diffusion processes: Algorithms for bimolecular reactions
,”
Phys. Biol.
6
,
046001
(
2009
).
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
).
92.
D.
Bernstein
, “
Simulating mesoscopic reaction-diffusion systems using the Gillespie algorithm
,”
Phys. Rev. E
71
,
041103
(
2005
).
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
).
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
).
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
).
96.
S.
Hellander
,
A.
Hellander
, and
L.
Petzold
, “
Reaction rates for mesoscopic reaction-diffusion kinetics
,”
Phys. Rev. E
91
,
023312
(
2015
).
97.
S. A.
Isaacson
, “
A convergent reaction-diffusion master equation
,”
J. Chem. Phys.
139
,
054101
(
2013
).
98.
Y.
Cao
and
R.
Erban
, “
Stochastic Turing patterns: Analysis of compartment-based approaches
,”
Bull. Math. Biol.
76
,
3051
3069
(
2014
).
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
).
100.
E.
Moro
and
H.
Schurz
, “
Boundary preserving semianalytic numerical algorithms for stochastic differential equations
,”
SIAM J. Sci. Comput.
29
,
1525
1549
(
2007
).
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
).
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
).
103.
W. L.
Briggs
,
V. E.
Henson
, and
S. F.
McCormick
,
A Multigrid Tutorial
, 2nd ed. (
SIAM
,
2000
).
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
).
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.
106.
R.
Ramadugu
,
S. P.
Thampi
,
R.
Adhikari
,
S.
Succi
, and
S.
Ansumali
, “
Lattice differential operators for computational physics
,”
EPL
101
,
50006
(
2013
).
107.
Z.
Jdrzejewski-Szmek
and
K. T.
Blackwell
, “
Asynchronous τ-leaping
,”
J. Chem. Phys.
144
,
125104
(
2016
).
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
).
109.
J. A.
Dix
and
A. S.
Verkman
, “
Crowding effects on diffusion in solutions and cells
,”
Annu. Rev. Biophys.
37
,
247
263
(
2008
).
110.
J.
Skolnick
, “
Perspective: On the importance of hydrodynamic interactions in the subcellular dynamics of macromolecules
,”
J. Chem. Phys.
145
,
100901
(
2016
).
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
.
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
).
113.
T.
Biancalani
,
D.
Fanelli
, and
F.
Di Patti
, “
Stochastic Turing patterns in the brusselator model
,”
Phys. Rev. E
81
,
046215
(
2010
).
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
).
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
).