 Research
 Open Access
 Published:
The Spatial Chemical Langevin Equation and Reaction Diffusion Master Equations: moments and qualitative solutions
Theoretical Biology and Medical Modellingvolume 12, Article number: 5 (2015)
Abstract
Background
It has been established that stochastic effects play an important role in spatiotemporal biochemical networks. A popular method of representing such stochastic systems is the Reaction Diffusion Master Equation (RDME). However, simulating sample paths from the RDME can be computationally expensive, particularly at large populations. Here we investigate an uncommon, but much faster alternative: the Spatial Chemical Langevin Equation (SCLE).
Methods
We investigate moment equations and correlation functions analytically, then we compare sample paths and moments of the SCLE to the RDME and associated deterministic solutions. Sample paths are generated computationally by the Next Subvolume method (RDME) and the EulerMaruyama method (SCLE), while a deterministic solution is obtained with an Euler method. We consider the GrayScott model, a wellknown pattern generating system, and a predator–prey system with spatially inhomogeneous parameters as sample applications.
Results
For linear reaction networks, it is well known that the first order moments of all three approaches match, that the RDME and SCLE match to the second moment, and that all approaches diverge at third order moments. For nonlinear reaction networks, differential equations governing moments do not form a closed system, but a general moment equation can be compared term wise. All approaches match at the leading order, and the RDME and SCLE match at the second leading order. As expected, the SCLE captures many dynamics of the RDME where deterministic methods fail to represent them. However, areas of the parameter space in the GrayScott model exist where either the SCLE and RDME give qualitatively different predictions, or the RDME predicts patterns, while the SCLE does not.
Conclusions
The SCLE provides a fast alternative to existing methods for simulation of spatial stochastic biochemical networks, capturing many aspects of dynamics represented by the RDME. This becomes very useful in search of quantitative parameters yielding desired qualitative solutions. However, there exist parameter sets where both the qualitative and quantitative behaviour of the SCLE can differ when compared to the RDME, so care should be taken in its use for applications demanding greater accuracy.
Background
Spatialstochastic effects are increasingly found to play important roles throughout a range of biological scales, from intracellular and intercellular processes, to ecological and epidemiological scales [13]. A widely used approach to study stochastic spatial dynamics is the Reaction Diffusion Master Equation (RDME), in which space is partitioned discretely into a number of voxels. Diffusion can then occur between different voxels, and reactions can occur within voxels on the assumption that reactants are wellmixed.
The RDME is generally analytically intractable. However, there do exist some closed form solutions for systems involving monomolecular reactions [4]. Meanfield approaches provide some analytical tools to help understand systems with bimolecular reactions [5,6], but these do not provide exact solutions. It is possible to generate sample paths consistent with the RDME using a variety of spatial Stochastic Simulation Algorithms (SSSAs) [79]. While widely used, SSSAs suffer from several drawbacks. For example, there are spatial resolution limits under which artefacts in particle interactions might occur [10], and also some effects at boundaries might not be accurately captured [11]. However, a more significant drawback is the fact that SSSAs are eventdriven algorithms. Thus, at large numbers of particles, the number of events per time step can become very large, and SSSAs become prohibitively slow. While one might argue that deterministic approaches might suffice in such regimes, it has been shown that stochastic effects can give rise to important effects here, such as noiseinduced oscillations and patterns [3,12]. There exist alternative algorithms based on the RDME which are faster in such scenarios, but these sacrifice the exactness of the SSSA [13], and thus are not guaranteed to faithfully represent the behaviour of the RDME. In this paper, we investigate one such alternative: the Spatial Chemical Langevin Equation (SCLE). In the nonspatial setting, the Chemical Langevin Equation (CLE), can be derived from the Chemical Master Equation (CME) [14], which in turn can be derived from a microscopic description of chemical processes [15]. The CLE and CME can then be extended to the SCLE and RDME, respectively, by introducing diffusion analogously to linear reactions. The SCLE consists of a family of stochastic differential equations (SDEs), which have the advantage that they can be simulated with fixed time steps, thus shedding the computational overhead associated with the eventdriven SSSAs. Furthermore, there are simple schemes available to simulate SDEs. However, very little work has been done on investigating the SCLE in detail (see [16,17]). Moreover, these studies do not incorporate diffusion terms in a manner consistent with linear reactions in the nonspatial case. Since diffusion might be viewed as a linear reaction, it is important to maintain a consistent formulation between diffusion and linear reactions when introducing diffusion into existing nonspatial models.
Differences between deterministic and Master Equation approaches have been wellstudied (e.g. [1820]). However, previous work on comparing Langevin and Master equations has concentrated on nonspatial settings using the CLE and the CME [21], and also in nonspatial systems with delay [22]. In particular, it has been found that for linear reaction networks, that the first two moments of CME and CLE match. It was also demonstrated numerically that the CME and CLE can give similar moments in a nonlinear reaction with a population of the order of 100 [21], but no formal proof extrapolating to different scenarios or populations was provided. While the method of adding spatial interactions is a straightforward extension involving adding linear interactions to the nonspatial case, there remain some open questions as to the consequences of such an implementation. For example, in a spatial model, it is possible for populations of some species to be high in some voxels, and low in others. Thus, it remains uncertain as to whether inaccuracies arising from voxels with small populations will propagate through the system. Furthermore, there is the interesting question as to how spatial correlation functions might behave in the spatial setting, as there is no such analogue in the nonspatial case.
In this work, we address the applicability of the SCLE as a substitute for the RDME. Since a key aspect of spatial models includes spatial correlations of different species, we investigate moments and correlation functions of the RDME and SCLE, and compare these with deterministic solutions for general reaction networks. In general, the moment equations for systems with nonlinear reactions do not form a closed set of equations, and cannot be solved without further assumptions on the underlying distribution of particle populations [23,24]. However, we can draw conclusions on the moments of the RDME and SCLE without closing this set of equations, namely by comparing each corresponding equation termwise. We provide a thorough numerical investigation of the GrayScott model to investigate whether the SCLE is capable of capturing pattern formation driven by intrinsic noise.
Results and discussion
Formulation
We construct a space Ξ, partitioned into N disjoint voxels, ξ _{1}, ξ _{2}, …, ξ _{ N }. Consider a number, S, of distinct chemical species: X _{1}, X _{2}, …, X _{ S }. Each of these can react according to some reaction, of which there are R in total. We denote the population of X _{ s } (s ∈ {1, …, S}) in each ξ _{ n } ∈ Ξ by \( {u}_n^{(s)} \). We can use these quantities to define the total population for each ξ _{ n } ∈ Ξ as \( {\mathbf{u}}_n={\left({u}_n^{(1)},{u}_n^{(2)},\dots, {u}_n^{(S)}\right)}^T \) (with T denoting the transpose of the vector) and also the total state of the system as an S x N matrix U = (u _{1}, u _{2}, …, u _{ N }). We allow the state of the system to change by reactions r ∈ {1, …, R}:
where α _{ sr } and β _{ sr } are natural numbers defining the stoichiometries of the reaction r, and \( {a}_r^{(n)}\left(\mathbf{U}\right) \) is the reaction propensity of reaction r occurring within voxel ξ _{ n }. Note that we set the propensity functions to be dependent on voxel attributes to allow for flexibility in our approach, e.g. to consider different sized voxels, or to introduce spatiallyinhomogeneous reaction rates. Using the reactions defined in (1), we can define a stoichiometric matrix V ^{(r,n)}, with the same dimensionality as U, such that the occurrence of reaction r in voxel ξ _{ n } causes a change of state from U to U + V ^{(r,n)}. We denote the value of V ^{(r,n)} at the entry corresponding to species s in voxel x by \( {v}_{s,x}^{\left(r,n\right)} \). Note that for typical systems, it will be the case that V ^{(r,n)} is sparse, i.e. \( {v}_{s,x}^{\left(r,n\right)}=0 \) whenever x ≠ n, which can be for most values of x. However, it is helpful to define the stoichiometric matrix on a larger space since it allows for more compact expressions in our subsequent analysis. We allow species s to diffuse from ξ _{ m } to ξ _{ n } with propensity \( {d}_{m,n}^{(s)}{u}_m^{(s)} \), where \( {d}_{m,n}^{(s)} \) represents the microscopic diffusion rate of species s from voxel m to voxel n. The change of total state as a result of such a diffusion event is given by a matrix \( {\mathbf{W}}_{m,n}^{(s)} \), which is of the same size as U _{,} has a value of 1 at the position corresponding to \( {u}_n^{(s)} \) and −1 at the entry corresponding to \( {u}_m^{(s)} \), and is zero everywhere else. We refer to the element of \( {\mathbf{W}}_{m,n}^{(s)} \) corresponding to species s’ in voxel x’ by \( {w}_{s\hbox{'},x\hbox{'}}^{\left(s,m,n\right)} \). In the proceeding analysis, for some collection of n (n ∈ ℕ) random variables, X _{1}, X _{2}, …, X _{ n }, we denote the expectation of \( {\displaystyle \prod_{i=1}^n{X}_i} \) by \( \left\langle {\displaystyle \prod_{i=1}^n{X}_i}\right\rangle \) and refer to this quantity as an n ^{th} order moment. Furthermore, for random variables X and Y, we define the correlation of X and Y to be 〈XY〉, which can also be alternatively viewed as a second order moment. We continue our analysis by considering the RDME, then the SCLE and deterministic approaches. Immediately after each approach, we present the moment equation implied by the corresponding approach. For the reader’s convenience, definitions of the RDME, SCLE and deterministic approaches are found in equations (2), (4) and (6), respectively. The corresponding correlation functions are in (3), (5) and (7), respectively.
Correlation functions from master equations
Using the definitions developed in the previous section, we can write out a corresponding RDME for the probability that the system is at some state u at a time t, i.e. P(U = u, t). We denote this by P(u, t) for brevity. This leaves us with the following form for the RDME:
where the first and second lines represent contributions of reaction and diffusion events, respectively. To investigate various spatial phenomena, we consider the correlation function of the populations of species x in voxels ξ _{ a } and ξ _{ b }, \( \left\langle {u}_a^{(x)}{u}_b^{(x)}\right\rangle \). In situations with nonlinear propensities, ODEs describing these correlation functions depend on higher order correlation functions, which in turn are related to correlation functions of everincreasing order. Thus, we consider a general correlation function of order O represented by \( \left\langle {\displaystyle \prod_{o=1}^O{u}_{n_o}^{\left({s}_o\right)}}\right\rangle \) where for each o we choose some n _{ o } ∈ {1, …, N} and s _{ o } ∈ {1, …, S}. We apply the operator \( {\displaystyle \sum_{\mathbf{u}}{\displaystyle \prod_{o=1}^O{u}_{n_o}^{\left({s}_o\right)}}} \) to both sides of (2) to give the following ODE (see Additional file 1 for derivation):
where the expression is structured as in (2), and the first and second lines represent contributions from reaction and diffusion events, respectively. Observe that the leading order of both \( {\displaystyle \prod_{o=1}^O\left({u}_{n_o}^{\left({s}_o\right)}+{v}_{s_o,{n}_o}^{\left(r,n\right)}\right)} \) and \( {\displaystyle \prod_{o=1}^O\left({u}_{n_o}^{\left({s}_o\right)}+{w}_{s_o,{n}_o}^{\left(s,m,n\right)}\right)} \) is \( {\displaystyle \prod_{o=1}^O{u}_{n_o}^{\left({s}_o\right)}} \). Consequently, the moments of the diffusion terms depend only on moments of order O and below. However, the order of the moments arising from reaction terms depends on the specific form of \( {a}_r^{(n)}\left(\mathbf{U}\right) \). If all such reaction propensities are linear or constant, then (3) forms a closed set of equations. Otherwise, the moments of order O are dependent on moments of order O + 1 or higher, and thus the set of ODEs in (3) is dependent on an infinite set of higher order ODEs, and cannot be solved without some extra assumptions on the process. For our purposes, we do not seek to solve (3), but only compare each ODE generated by (3) termwise to corresponding ODEs obtained from the SCLE. This then provides some insight into how the spatial correlations compare between these approaches.
Correlation functions from Langevin and deterministic equations
We consider the Langevin and deterministic regimes together in this section, since their analytical treatment is similar. Using the same established definitions from the case of the RDME, we can write out a set of SDEs for each chemical species in each voxel. Taken together, these provide a Langevin representation, which we take to be a counterpart to the RDME of the previous section. To arrive at such representations, one can either start from the Chemical Langevin Equation in the nonspatial setting (as in [21]) and introduce diffusion between voxels analogously to linear reactions. An alternative way to arrive at Langevin representations is to proceed directly by an expansion of the RDME [25]. Note that such representations assume that the population of various quantities are continuous variables, as opposed to the RDME approaches which preserve discreteness. The SCLE consists of the following coupled SDEs for the population of molecules of each species s _{ o } in each voxel \( {\xi}_{n_o} \):
where each \( {W}_r^{(n)} \) is a Wiener process used to represent noise in the occurrence of reaction r in voxel ξ _{ n } . Similarly, each \( {W}_{mn}^{(s)} \) incorporates the stochastic nature of diffusion of species s from voxel m to voxel n. We refer to the collection of SDEs represented by (4) as the SCLE.
An expression for an ODE describing \( \left\langle {\displaystyle \prod_{o=1}^O{u}_{n_o}^{\left({s}_o\right)}}\right\rangle \) follows by applying Ito’s Lemma to (4) (see Additional file 1):
where we have structured the expression in the same way as in (3), with the contribution from reaction and diffusion appearing in the first and second lines, respectively. The expression becomes clearer when seen in terms of binomial expansions of various terms in (3). The reaction terms in (5) are identical to the leading two orders of the reaction terms appearing in the expansion of (3). The same holds for the diffusion terms. Thus, each of the infinite set of ODEs representing correlation functions from the RDME and SCLE match up to the second leading order.
For the deterministic model, we consider the SCLE without any contribution from Wiener terms. That is, for each s _{ o } in voxel \( {\xi}_{n_o} \), we consider the following deterministic representation:
The ODEs describing correlation functions and moment equations from (6) can be easily calculated, since in the deterministic scenario, expectation and products commute, i.e., \( \left\langle {\displaystyle \prod_{o=1}^O{u}_{n_o}^{\left({s}_o\right)}}\right\rangle ={\displaystyle \prod_{o=1}^O\left\langle {u}_{n_o}^{\left({s}_o\right)}\right\rangle } \). However, to facilitate comparison with the RDME and SCLE, we consider ODEs of the same form as (3) and (5). The result follows directly from the derivation of (5) (see Additional file 1 for details):
which lends itself to the same interpretation as (5). The deterministic representation only accounts for leading order terms of the correlation functions and moment equations.
Implications of moment equations
From equations (3), (5) and (7), we can make some analytical and also some heuristic predictions on the behaviour of the RDME, SCLE and deterministic approaches. First note that all moment equations match at the leading order, and the RDME and SCLE match at the second leading order.
In the case of linear networks (i.e. \( {a}_r^{(n)}\left(\mathbf{U}\right) \) is linear in U for all r and n), we have that all three approaches give rise to a set of moment equations, which are closed, and can thus be solved. Since all three approaches agree at the leading order, the mean population of all three approaches match. Similarly, since the RDME and SCLE match at second leading order, we can conclude that the variance of populations between the RDME and SCLE agree for linear networks, while the same cannot be said for deterministic approaches. Since all three approaches deviate at the third leading order, we anticipate that the three approaches will not agree at the level of third order moments and higher. These results have previously been considered in the nonspatial case [21], and their extension to the spatial case in this work is to be expected, since diffusion in space can be considered as a linear reaction in itself.
In the case that there are second order or higher terms of species populations appearing in \( {a}_r^{(n)}\left(\mathbf{U}\right) \) for some r and n, then the equations (3), (5) and (7) do not form a closed set, since each equation of order O is dependent on orders of O + 1 or higher. Thus, in the system of (infinite) ODEs describing these moments, there will be terms where all three approaches disagree, since they match at most at the leading two orders. Thus, moments of populations with nonlinear reaction propensities cannot be expected to match. These arguments also hold when applied to the spatial correlation functions within the system, since spatial correlation functions can be defined in terms of moments.
In the limit of large populations, each ODE is dominated by its leading order terms, and thus one can expect moments and correlation functions to match for systems with very large populations. At smaller populations we anticipate that the SCLE should provide a better approximation to the RDME than a comparable deterministic approach, since the RDME and SCLE match at two leading order terms, whereas a deterministic approach only matches at the leading order term. However it is as yet unclear at what numbers the SCLE and RDME might give similar moments. This problem is compounded by the fact that in spatial scenarios, it is possible for some areas of space to have large populations, with others having small populations, such as the case in Turing patterns. Thus, we proceed with a numerical investigation.
The GrayScott model
On comparing the ODEs governing spatial correlation functions for the RDME, SCLE and deterministic approaches, it is clear that the same spatial correlations should not be expected from the three approaches, but it is uncertain when one approach should wellapproximate the other. To provide an application with which we can illustrate the potential and limitations of each approach, we consider the GrayScott model [26], a widely used pattern generating system. While it was originally intended as a model of glycolysis [27], it has been shown that the model can generate several different patterns within a narrow parameter range [28], and that intrinsic noise can drive these patterns [12]. As such, it provides a good framework with which to probe differences between the three approaches.
For this we consider two reacting species, U and V. Denote the population of U and V in voxel ξ _{ i } by U _{ i } and V _{ i }, respectively. The GrayScott model is characterised by the following reactions occurring within each ξ _{ i } ∈ Ξ:
with reaction propensities within each voxel being k _{1} U _{ i } V _{ i }(V _{ i } − 1)/Ω^{2}, k _{2} V _{ i }, k _{3} U _{ i } and k _{4}Ω, respectively. The Ω term serves as a parameter to vary the system size. We conduct numerical experiments in a 2D Cartesian square. Each side of the square consists of L (L ∈ ℕ) square voxels, each with length h, giving the system a total length of Lh per side. We refer to the diffusion constants of U and V as D _{ U } and D _{ V }, respectively. Thus, the diffusion propensities of U and V from voxel ξ _{ i } to neighbouring voxels are given by D _{ U } U _{ i }/h ^{2} and D _{ V } V _{ i }/h ^{2}, respectively. In order to conduct numerical experiments, it is helpful to parameterise the reaction rates. Thus, we define the reaction rates in terms of two parameters k and F, and define k _{1} = 1, k _{2} = F + k, k _{3} = F and k _{4} = F. This parameterisation has been chosen so as to be consistent with previous studies of the GrayScott Model [12,28]. For all numerical investigations, we take Ω = 250, D _{ U } = 2 × 10^{− 5}, D _{ V } = 1 × 10^{− 5}, h = 0.01 and L = 50. For initial conditions, we use U = 250 and V = 0 everywhere, except in a centred box of 3×3 voxels, where we use U = 0 and V = 250. We operate in arbitrary length and time units in order to maintain consistency with previously published studies, which have followed this strategy. In order to simulate sample paths from the RDME, we implemented the next subvolume method [29]. To generate comparable paths from the SCLE and the deterministic approaches, we implemented the EulerMaruyama method [30] and an Euler method, respectively. As expected, the SCLE simulations were executed much faster than corresponding RDME simulations. The typical computational time to simulate 2500 time units from the RDME was on the order of two days, whereas the corresponding time from the SCLE was on the order of half an hour (with a time step of 0.1 time units), as computed using MATLAB R2013a running on a 3.2 GHz Intel Xeon processor. Reducing the time step on the EulerMaruyama and Euler methods to 0.025 time units was found to have no impact on the results.
In the area of parameter space where all three approaches generated patterns, typical patterns from the SCLE resembled something between the RDME and deterministic approaches (see Figure 1). That the qualitative patterns are not necessarily the same is not surprising, given one cannot expect the same correlation functions from all three approaches.
However, it is not always the case that the three approaches predict the existence of patterns for the same parameters (see Figure 2). We observe three different scenarios: 1) All three approaches predict patterns, 2) Only the RDME and SCLE predict patterns and 3) Only the RDME predicts patterns. Such a result might be expected from the analytical expressions themselves. The region of parameter space where the RDME and the SCLE match is sizeable, and only in a few of the investigated parameters did the RDME predict patterns where the SCLE did not (see Figure 3). Our numerical results indicate that the RDME and SCLE give qualitatively similar results in the vast majority of the parameter space, however, there are cases where it can give qualitatively different results.
A predator–prey system with a safe haven
To investigate first and second moments of a nonlinear reaction system numerically, we consider a modified predator prey system in a space of 5 by 5 voxels, where prey are safe from predation in a small area. Such models can provide insight into the wider impact of implementing conservation schemes in local areas. We denote prey and predator species by A and B, respectively, and write their populations in voxel i as A _{ i } and B _{ i }, respectively. The following reactions occur in every voxel:
with corresponding propensities of k _{1} A _{ i }, k _{2} A _{ i } B _{ i }, k _{3} A _{ i } and k _{4} B _{ i }, respectively. In voxel (1,1), we set k _{ 2 } = 0, signifying that prey are safe from predators in this voxel. Furthermore, both A and B diffuse between voxels with the same diffusion rate k _{ diff }.
Sample paths from the RDME were generated by using the next subvolume method. Corresponding paths from the SCLE and deterministic approaches were generated by using the EulerMaruyama and Euler methods, respectively. The population of B in voxels (1,1), (3,3) and (5,5) are shown in Figure 4, which shows first and second moments estimated from 10,000 realizations. The following parameters were used in the presented simulations: k _{ 1 } = 1, k _{ 2 } = 0.02, k _{ 3 } = 0.25, k _{ 4 } = 1, and k _{ diff } = 0. Initial conditions of 100 individuals of A and 100 individuals of B were set in every voxel. At all times, values of the SCLE were found to interpolate those of the RDME and deterministic approaches. For example, at 10 seconds, the values of the second moment of B in voxel (1,1) from the RDME, SCLE and deterministic approaches was 85100, 84200 and 81800 rounded to three significant figures, respectively.
Discussion
We conducted an investigation of the SCLE both analytically and numerically, with an emphasis on comparing moments, correlation functions and qualitative behaviour to the RDME and deterministic approaches. We show that for systems with linear reaction networks, the SCLE provides correct descriptions of first and second order moments, but not for third moments and higher. For nonlinear cases, it cannot be guaranteed that the moments match.
In such nonlinear scenarios, ODEs describing moments and correlation functions do not form a closed system of equations, but depend on higher order moment equations. Thus, these systems are not solvable without further assumptions being imposed on the system [23,24]. However, by comparing each equation termwise, we showed that the RDME and SCLE match at the leading and secondleading orders, whereas deterministic approaches match only at the leading order. The implications of this depend on the specific reaction propensities in the network of interest. For linear networks, the RDME and SCLE match at the first and second moments, whereas deterministic approaches only represent first moments accurately. For nonlinear networks, little can be said conclusively. These results are summarised in Table 1. Ascertaining the population size where this would lead to qualitative differences could not be performed analytically, so a numerical investigation was performed instead.
Spatial studies of a predator–prey system showed that the SCLE can provide first and second order moments which closely match those of the RDME for populations on the order of 100 individuals (see Figure 4). These findings reinforce what was previously investigated in nonspatial settings [21]. The moments of the SCLE were found to interpolate between that of deterministic approaches and the RDME in all numerical investigations.
To demonstrate the applicability of the SCLE in capturing phenomena driven by intrinsic noise, we considered simulations of the GrayScott model [26,31]. While the qualitative solutions between the RDME and SCLE were typically similar, there were cases where they differed (see Figures 1 and 2). In particular, there were even some regions of parameter space where patterns might be observed in the RDME, but not in the SCLE (see Figures 2 and 3). Where the RDME, SCLE and deterministic approaches all predicted patterns, it was interesting to note that the resulting patterns obtained from the SCLE appeared to be an intermediate between patterns associated with the RDME and the PDE solutions.
Such results were conducted for populations of the order of a few hundred particles. In smaller populations, it is clear from equations (3), (5) and (7) that we cannot expect the moments nor correlation functions to match. It has recently been shown in the nonspatial case the CLE can be interpreted to give complex values for nonlinear reactions at small populations [32], thus making it problematic to compare sample paths from master equations and Langevin equation at small populations, since it would involve comparing real and complex numbers with one another. In the spatial case, another artefact might appear in the SCLE: the fact that the SCLE admits continuous valued population sizes means that the notion of a particle being contained entirely within one voxel might be lost.
A key advantage the SCLE has over the RDME approach is computational efficiency. In executing the simulations in making Figures 1 and 3, the computational time required was on the order of days for simulations from the RDME, whereas analogous simulations were completed in the order of hours from the SCLE. The computational savings can be even more significant for larger systems, since SSSAs scale in computational speed according to a polynomial of the population size, whereas the speed of the EulerMaruyama method is independent of the population size. As such, the main benefit of the SCLE is reached in the region where it gives the best accuracy, making it especially suitable for stochastic simulation on the macroscopic scales and parameter sweeps.
Conclusions
We demonstrated that the SCLE captures many qualitative and quantitative characteristics of the RDME, which deterministic models fail to represent. The SCLE was found to be significantly faster to simulate than the RDME. Qualitative differences in behaviour of the RDME to the SCLE occurred for specific parameter settings, but such occurrences were uncommon. Analytical results demonstrate that the SCLE matches the RDME most at large population sizes. Thus, we anticipate that the SCLE should provide a strong framework for the simulation of reaction–diffusion systems with many particles.
Abbreviations
 CLE:

Chemical Langevin Equation
 CME:

Chemical Master Equation
 RDME:

Reaction Diffusion Master Equation
 SCLE:

Spatial Chemical Langevin Equation
References
 1.
McAdams HH, Arkin A. Stochastic mechanisms in gene expression. Proc Natl Acad Sci U S A. 1997;94:814–9.
 2.
Maini PK, Painter KJ, Chau HNP. Spatial pattern formation in chemical and biological systems. J Chem SocFaraday Transac. 1997;93:3601–10.
 3.
Hempel H, SchimanskyGeier L, GarciaOjalvo J. Noisesustained pulsating patterns and global oscillations in subexcitable media. Phys Rev Lett. 1999;82:3713–6.
 4.
Jahnke T, Huisinga W. Solving the chemical master equation for monomolecular reaction systems analytically. J Math Biol. 2007;54:1–26.
 5.
Cardy JL, Tauber UC. Field theory of branching and annihilating random walks. J Stat Phys. 1998;90:1–56.
 6.
Dodd PJ, Ferguson NM. A manybody field theory approach to stochastic models in population biology. PLoS One. 2009;4:e6855.
 7.
Baras F, Mansour MM. Reaction–diffusion master equation: a comparison with microscopic simulations. Phys R E. 1996;54:6139–48.
 8.
Elf J, Doncic A, Ehrenberg M, Frauenfelder H, Moss F. Mesoscopic reaction–diffusion in intracellular signaling. In: Society of PhotoOptical Instrumentation Engineers (SPIE) Conference Series. Bellingham, Washington, USA: SPIE; 2003.
 9.
Erban R, Chapman J, Maini P. A practical guide to stochastic simulations of reaction–diffusion processes. 2007.
 10.
Isaacson SA. The reaction–diffusion master equation as an asymptotic approximation of diffusion to a small target. Siam J Appl Math. 2009;70:77–111.
 11.
Leier A, MarquezLago TT. Correction factors for boundary diffusion in reaction–diffusion master equations. J Chem Phys. 2011;135(13):134109. http://scitation.aip.org/content/aip/journal/jcp/135/13/10.1063/1.3634003.
 12.
Wang HL, Fu ZP, Xu XH, Qi OY. Pattern formation induced by internal microscopic fluctuations. J Phys Chem A. 2007;111:1265–70.
 13.
MarquezLago TT, Burrage K. Binomial tauleap spatial stochastic simulation algorithm for applications in chemical kinetics. J Chem Phys. 2007;127(10):104101. http://scitation.aip.org/content/aip/journal/jcp/127/10/10.1063/1.2771548.
 14.
Gillespie DT. The chemical Langevin equation. J Chem Phys. 2000;113:297–306.
 15.
Gillespie DT. A rigorous derivation of the chemical master equation. Phys A. 1992;188:404–25.
 16.
Andrews SS, Arkin AR. Simulating cell biology. Curr Biol. 2006;16:R523–7.
 17.
Andrews SS, Dinh T, Arkin AP. Stochastic Models of Biological Processes. In: Meyers RA, editor. Encyclopedia of Complexity and Systems Science. New York: Springer; 2009. p. 8730–49.
 18.
Zheng Q, Ross J. Comparison of deterministic and stochastic kinetics for nonlinearsystems. J Chem Phys. 1991;94:3644–8.
 19.
Grima R. NoiseInduced Breakdown of the MichaelisMenten Equation in SteadyState Conditions. Phys Rev Lett. 2009;102(21):218103. http://www.ncbi.nlm.nih.gov/pubmed/19519139.
 20.
Barrio M, Burrage K, Leier A, Tian TH. Oscillatory regulation of hes1: discrete stochastic delay modelling and simulation. Plos Comput Biol. 2006;2:1017–30.
 21.
Khanin R, Higham DJ. Chemical Master Equation and Langevin regimes for a gene transcription model. Theor Comp Sci. 2008;408:31–40.
 22.
Tian TH, Burrage K, Burrage PM, Carletti M. Stochastic delay differential equations for genetic regulatory networks. J Comput Appl Math. 2007;205:696–707.
 23.
Krishnarajah I, Cook A, Marion G, Gibson G. Novel moment closure approximations in stochastic epidemics. Bull Math Biol. 2005;67:855–73.
 24.
Lee CH, Kim KH, Kim P. A moment closure method for stochastic reaction networks. J Chem Phys. 2009;130(13):134107. http://www.ncbi.nlm.nih.gov/pubmed/19355717.
 25.
Kampen NG. Stochastic processes in physics and chemistry/N.G. van Kampen. Amsterdam; New York. New York: NorthHolland; sole distributors for the USA and CanadaElsevier NorthHolland; 1981.
 26.
Gray P, Scott SK. Autocatalytic reactions in the isothermal, continuous stirred tank reactorisolas and other forms of multistability. Chem Eng Sci. 1983;38:29–43.
 27.
Selkov EE. Selfoscillations in glycolysis .1. A simple kinetic model. Eur J Biochem. 1968;4:79.
 28.
Pearson JE. Complex Patterns in a Simple System. Science. 1993;261:189–92.
 29.
Elf J, Ehrenberg M. Spontaneous separation of bistable biochemical systems into spatial domains of opposite phases. Syst Biol. 2004;1:230–6.
 30.
Higham DJ. An algorithmic introduction to numerical simulation of stochastic differential equations. Siam Rev. 2001;43:525–46.
 31.
Gray P, Scott SK. Sustained oscillations and other exotic patterns of behavior in isothermal reactions. J Phys Chem. 1985;89:22–32.
 32.
Schnoerr D, Sanguinetti G, Grima R. The complex chemical Langevin equation. J Chem Phys. 2014;141:024103.
Acknowledgments
AG and TML were supported by subsidy from the Cabinet Office of Japan to the Integrative Systems Biology Unit (MarquezLago lab), Okinawa Institute of Science and Technology. AL was supported by subsidy from the Cabinet Office of Japan to the Okinawa Institute of Science and Technology.
Author information
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
TML and AL conceived and designed the study. AG performed derived mathematical expression and executed numerical simulations. All authors took part in the writing of this manuscript. All authors read and approved the final manuscript.
Additional file
Additional file 1:
Derivations of correlation functions from the Reaction Diffusion Master Equation and the Spatial Chemical Langevin Equation.
Rights and permissions
About this article
Received
Accepted
Published
DOI
Keywords
 Reaction Diffusion Master Equation
 Spatial Chemical Langevin Equation
 Turing patterns
 Noiseinduced phenomena