## Abstract

In small-scale chemical reaction networks, the local density of molecules is changed by discrete jumps owing to reactive collisions, and through transport. A systematic perturbation scheme is developed to analytically characterize these non-equilibrium intrinsic fluctuations in a multispecies spatially varying system. The method is illustrated on a variety of model systems. In all cases, the continuous approximation method is corroborated with extensive stochastic simulation. As an example of our technique applied to a spatially varying steady state, we demonstrate that a model for embryonic patterning mediated by regulatory mRNA is surprisingly robust to intrinsic fluctuations.

## 1. Introduction

For intracellular biochemical reaction networks, elucidating the mechanisms and physical constraints underlying spatial patterning is essential for understanding important phenomena such as bacterial cell division (Fange & Elf 2006) and embryonic development (Gregor *et al.* 2007). At a macroscopic level, the modelling of patterns typically takes the form of a partial differential equation, with a local reaction term characterizing the interaction among molecules and a transport term governing movement in the spatial domain. Classical mechanisms of pattern formation (Turing instabilities for diffusive transport, Turing 1953; Murray 2007; differential-flow instabilities for advective transport, Malchow *et al.* 2008) are useful in principle, but suffer from serious shortcomings in application to physical processes. Both mechanisms often require unphysically large disparity among the transport rates, and instabilities are predicted to exhibit perfect ordering with patterns generated at a fixed wavenumber. In reality, imperfect patterning is invariably observed. Despite these obvious shortcomings, models of the macroscopic dynamics have had much success in describing the mean behaviour of these spatially varying systems; however, this approach implicitly assumes that the concentrations are sufficiently large for the continuum approximation to be valid. Indeed, this classical approach is significantly less accurate in regions of low concentrations.

Although chemical and biological networks are often modelled using continuous concentration fields, when viewed at the microscopic level, it is readily observed that the number density is in fact discrete, leading to fluctuations about the mean-field behaviour. In a stochastic formulation, it is useful to categorize the origin of the fluctuations into two broad classes (van Kampen 1992; Swain *et al.* 2002): intrinsic (or demographic) noise arising from the discrete stochastic interactions of the individual constituents, and extrinsic (or environmental) noise arising from un-modelled dynamics outside the system of interest. For a spatially varying stochastic model, fluctuations are included in the governing equations in one of two ways: either as multiplicative noise in a stochastic partial differential equation (Garcia-Ojalvo & Sancho 1999), or by representing the dynamics in terms of a probability distribution for the state using a reaction-transport master equation (van Kampen 1992; Baras & Mansour 1998; Gardiner 2004).

Considerable progress has been made in characterizing, approximating and simulating realizations of stochastic partial differential equations, typically in the context of modelling extrinsic, or environmental, noise (Garcia-Ojalvo & Sancho 1999; Shardlow 2005; Malchow *et al.* 2008). Less work has been directed towards the approximation of the reaction-transport master equation. Emerging experimental methods point to the significant role intrinsic fluctuations (or more specifically, the control of these fluctuations) play in embryonic development (Gregor *et al.* 2007) and front propagation (Kessler & Levine 1998; Hansen *et al.* 2006; Hallatschek *et al.* 2007). Furthermore, in spatially homogeneous models, it is well known that intrinsic fluctuations can give rise to a rich variety of novel behaviour not observed in corresponding deterministic models (Wio & Lindenberg 2003; Scott *et al.* 2007). The most dramatic and unexpected examples are those in which fluctuations generate an ordered response: regular oscillations out of deterministically stable systems (Steuer *et al.* 2003; McKane & Newman 2005; Kuske *et al.* 2007; Pineda-Krch *et al.* 2007), or stochastic switching among deterministically stable fixed points (Walczak *et al.* 2005). In spatially varying models, intrinsic noise coupled with the additional spatial degrees of freedom open up further possibilities. These include the appearance of new fixed points (Togashi & Kaneko 2004), phase separation in bistable systems (Elf & Ehrenberg 2004) and, importantly, pattern formation in reaction–diffusion models that are deterministically stable (Howard & Rutenberg 2003; McKane & Newman 2004; Fange & Elf 2006; Butler & Goldenfeld 2009; Biancalani *et al.* 2010).

Several powerful numerical methods have been developed to generate stochastic simulation of the reaction-transport master equation (Ander *et al.* 2004; Elf & Ehrenberg 2004; Lampoudi *et al.* 2009). Nonetheless, for a systematic mapping of parameter space, simulation methods rapidly become prohibitively time-consuming, even for reasonably simple models. Noise-induced spatial structure has been demonstrated in patch-dynamic models by McKane and coworker (Lugo & McKane 2008), where a maximum occupancy density is enforced for a given patch of the discretized domain. Under those conditions, anomalous nonlinear cross-diffusion terms appear, introducing coupling among species not observed in classical mean-field models. Although these Moran-like processes are certainly appropriate for many ecological and evolutionary applications (e.g. Nowak 2006), they appear to be less applicable to open chemical systems where maximum closed-packing density is never approached. Other work focusing on noise-induced spatial structure in open systems maps the dynamical equations to a Lagrangian formalism (Butler & Goldenfeld 2009), although diffusion is treated as a deterministic transport process, leading to an underestimation of the magnitude of the intrinsic fluctuations (§4*a*). More recent work by Biancalani *et al.* (2010) describes a remarkable feature of stochastic spatially varying models—that noise-induced patterns can appear in a subset of reactants (§4*b*). In all of this previous work, the equilibrium fluctuations are studied via Fourier transform, and are therefore limited to systems exhibiting a spatially homogeneous steady state.

Here, we develop a general approximation scheme for the solution of the multi-species reaction-transport master equation, built upon previous work by van Kampen for the single-species case (van Kampen 1976*b*, 1992), that extends and unifies past work. The governing equations are easily constructed once the details of the stochastic model are specified. We use this method to compute exactly the non-equilibrium standard-deviation envelope for a concentration field undergoing purely diffusive transport, illustrating an error in recent work (Butler & Goldenfeld 2009). Around a homogeneous steady state, the equation for the second moment is linear so that the equilibrium spatial correlation spectrum is determined by a set of algebraic equations. The spectra of the resulting solutions are used to identify regions of parameter space exhibiting noise-induced spatial structure. We apply this method to a two-species activator–inhibitor model to illustrate the differences that can arise in the fluctuating system in comparison to the mean-field approach. We show that the inherent stochasticity in the system serves as a destabilizing influence. Consequently, the parameter regime for diffusion-driven instabilities is enlarged, and the patterns are formed over a range of wavelengths. As in the work of Biancalani *et al.* (2010), we identify a parameter regime where noise-induced patterning is only apparent in one of the species, behaviour that cannot occur in the deterministic model. As an example of our technique applied to a spatially varying steady state, we examine the magnitude of intrinsic fluctuations in a model used to describe RNA-mediated spatial patterning in development (Levine *et al.* 2007). We demonstrate that, in addition to the sharp patterning exhibited by the deterministic model, the mechanism is surprisingly tolerant to intrinsic fluctuations. This is, to the best of our knowledge, the first time that the intrinsic fluctuations about a spatially varying steady state have been estimated analytically.

## 2. Non-equilibrium characterization of spatially varying intrinsic noise

Our goal is to compute the mean and variance owing to intrinsic fluctuations in a time- and space-varying concentration field. To that end, we artificially discretize the continuous spatial domain into subvolumes and apply the linear noise approximation in each subvolume (van Kampen 1976*a*). This discretization is a useful short cut that allows us to avoid functional integration, but at the expense of rather dense notation. We derive a large system of coupled differential equations that govern the evolution of the first two moments. In contrast to other methods (Lugo & McKane 2008), in the final stage, the discretization is reversed by taking the continuum limit, leading to a small set of partial differential equations governing the time-dependent mean and variance fields. As detailed below, the approximation relies on van Kampen’s system size expansion (van Kampen 1976*a*), which accommodates fluctuations owing to linear transport exactly, and therefore has a range of applicability that dependents only on errors incurred by the approximation of the reaction events (Grima 2009).

### (a) Approximate evolution of the first and second moments

If is the number of each species in cell *λ* with volume *Ω*, then the master equation governing the probability density function, *P*(*N*^{λ},*t*), is written in analogy with the spatially homogeneous case,2.1where transport is included as a linear reaction event moving the species between cells (figure 1). Here, the reaction rates (also called the reaction propensities) have units of concentration per time, **S** is the stoichiometry matrix specifying the network connectivity and is the step operator incrementing the *i*th species by an integer *k*: (van Kampen 1976*b*, 1992; Elf & Ehrenberg 2003). For nonlinear transition rates *ν*_{j}(**N**(*r*_{λ})/*Ω*), there is no known general solution for the master equation, necessitating methods of approximation. To that end, we make the linear noise approximation derived from the ansatz that inside each cell,2.2where, as above, *Ω* is the volume of each cell *λ*. The total number *N*^{λ} is written in terms of a deterministic concentration **x**^{λ} and a continuous stochastic variable *Ω*^{1/2}*α*^{λ} that scales as the square root of the number of particles (van Kampen 1961, 1976*a*). Within each subvolume, in the limit , the concentration **x**^{λ} is governed by the deterministic rate equation,2.3The fluctuations *α*^{λ} owing to the reaction kinetics are characterized by the probability density function *p*^{λ}(*α*^{λ},*t*) governed by the linear Fokker–Planck equation,2.4Summation over repeated Roman indices is implied (Greek indices denote location in space). The matrices **A** and **B** are defined by the rate vector ** ν** and stoichiometry matrix

**S**as follows:2.5Here, diag[

**x**] is a square matrix with the vector

**x**along the diagonal (see Elf & Ehrenberg (2003) for the derivation of the multivariate, spatially homogeneous Fokker–Planck equation). The matrix

**A**is simply the Jacobian of the deterministic model evaluated along the deterministic trajectory, and quantifies the local dissipation in the model, while the matrix

**B**is related to the magnitude of the fluctuations (Scott

*et al.*2006). In the absence of transport, the entire joint probability density function extending over all subvolumes

**={ … ,**

*α*

*α*^{λ−1},

*α*^{λ},

*α*^{λ+1}, … } is simply2.6

In practice, our focus is not on the joint distribution *p*(** α**,

*t*) but rather on the first two moments of the distribution, as these are the most accessible through experiments. Furthermore, taking the continuum limit of the Fokker–Planck equation is not a well-defined operation; in contrast, the continuum limit of the moment equations is perfectly well defined (van Kampen 1992).

To compute the evolution of the mean, multiply the time derivative of equation (2.6) by and using equation (2.4), integrate over all species and subvolumes2.7

Integrating by parts with respect to *α*^{λ}, the second term in the square brackets vanishes. The first term contributes only one non-zero element when *i*=*k* and *λ*=*μ*,2.8

The calculation of the second moment proceeds in a similar fashion— multiply equation (2.6) by and integrate by parts:2.9The first term in square brackets produces two terms, corresponding to *i*=*k* and *λ*=*μ*, and *i*=*l* and *λ*=*η*,2.10The second term in square brackets likewise produces two terms, corresponding to *i*=*k*,*j*=*l*, and *i*=*l*,*j*=*k* (and *λ*=*μ*=*η*),2.11where *δ*_{μη} is the Kronecker delta. Combining equations (2.10) and (2.11), in matrix form, the second moments **C**^{μη}=〈*α*^{μ}(*α*^{η})^{T}〉 obey the following system of ordinary differential equations:2.12

Using equations (2.8) and (2.12), it is straightforward to show that the second cumulants2.13obey the same equation (2.12) as the second moments. It is convenient in what follows to re-write the matrices **A**^{μ} and **B**^{μ} in terms of the propensity vector ** ν**(

**x**

^{μ}) and stoichiometry matrix

**S**introduced in the master equation, equation (2.1), so that the equation for the second cumulants is given by2.14Note that the cumulant matrix is no longer symmetric,

**K**

^{μη}≠

**K**

^{ημ}, as it is in the spatially uniform case when

**A**

^{μ}=

**A**

^{η}. Also note that the inhomogeneous term on the right-hand side appears locally (i.e. when

*μ*=

*η*), and is therefore multiplied by the Kronecker delta function.

Equation (2.14) characterizes the fluctuations owing to the jump-driven reaction kinetics. In a spatially varying system, local fluctuations can be transported to neighbouring subvolumes. In addition, diffusion itself acts as a source of fluctuations (§4*a*). If the transport is deterministic, then the fluctuations are transported by the same operator, and we simply append a transport term to the continuum limit of equation (2.14) above. Diffusion is, however, more appropriately modelled as a stochastic process, and as such, the transport of the variance becomes more complicated (van Kampen 1976*b*). Rather than derive an equation for the transport of the variance, it has been previously shown in van Kampen (1976*b*, 1992) that it is more compact to introduce active and diffusive transport to the equation for the factorial cumulant ** Ξ**,2.15A useful property of the factorial cumulant is that for the Poisson distribution, second and higher factorial cumulants vanish, making them particularly useful for characterizing the fluctuations owing to diffusive transport. By definition, the concentration in each cell

**x**

^{λ}is a sure variable (i.e. 〈〈

**x**

**x**

^{T}〉〉=0), and substituting equation (2.2) into the above equation yields an expression for the factorial cumulants in terms of the fluctuations,

*α*^{λ}, about the deterministic state:2.16To leading order in

*Ω*, from equation (2.3), the average number of each species in subvolume

*λ*is governed by the deterministic rate equation2.17Combining this with the evolution equation for the second cumulant of the fluctuations, equation (2.14), we obtain the evolution equation for the factorial cumulant:2.18

Up to this point, the system is described in terms of individual reaction subvolumes distributed on a discrete lattice. This large system of coupled equations is cumbersome to deal with, except at a homogeneous steady state where they can be converted straightforwardly to the Fourier domain for analysis (Lugo & McKane 2008). For a non-equilibrium analysis, it is more convenient to express the equations for the moments in a continuous representation, as an averaged concentration and fluctuation field. Upon taking the continuum limit, the indexing by the Greek letters gives way to a continuous domain denoted by the spatial coordinate **r**: **N**^{λ} → **N**(**r**). Notice that the continuum limit is not another level of approximation—it simply reverses the artificial discretization invoked to simplify the derivation of the governing equations. The applicability of the approximation still lies in the validity of the assumption that the local number density is sufficiently high that a perturbation expansion in inverse molecule numbers is meaningful.

The continuum limit is achieved by dividing the factorial cumulant *Ξ*^{μη} in equation (2.15) by *Ω*^{2}, and taking the limit of the lattice spacing going to zero to obtain the continuous limit of the factorial cumulant,where **n**(**r**) is the number density at location **r**, and subvolumes *μ* and *η* are centred on the points **r**_{1} and **r**_{2}, respectively. Using the fact that *δ*_{μη}/*Ω* → *δ*(**r**_{1}−**r**_{2}) as *Ω* → 0, equation (2.18) becomes2.19If, in addition to reaction events, the local population is changed through transport governed by a linear operator ** L**(

**r**), then the continuous limit of the mean-field equation (2.17) becomes2.20

Transport of the cumulant factorial ** Ξ** is somewhat more elaborate—we introduce the transport operator that is defined component-wise in terms of the mean-field transport

**(**

*L***r**) through the Hadamard product (Horn & Johnson 1985),2.21The advantage of using the factorial cumulants (instead of other expressions involving the second moment) is that the diffusive transport acts in an intuitive way: the total transport of the factorial cumulant,

*Ξ*

_{ij}(

**r**

_{1},

**r**

_{2})=[

*n*

_{i}(

**r**

_{1})

*n*

_{j}(

**r**

_{2})], is the sum of the transport of the first term,

*n*

_{i}(

**r**

_{1}), with respect to

**r**

_{1}and the transport of the second term,

*n*

_{i}(

**r**

_{2}), with respect to

**r**

_{2}. This macroscopic transport operator can be related to the microscopic transport rates, and indeed, it is that correspondence that motivates the use of the somewhat unfamiliar factorial cumulant; further details can be found in van Kampen (1976

*b*, 1992).

In terms of , the evolution of the factorial cumulants is governed by a system of linear partial differential equations that combine the effect of stochasticity in the reaction kinetics and in the diffusive transport, both to relative order *O*(*Ω*^{−1/2}),2.22This equation is the central result of this work. The leading order term in the perturbation expansion simply recovers the mean-field reaction–diffusion equations (equation (2.20)) that arise in classical deterministic models. The importance of equation (2.22) is that it determines the evolution of the magnitude of the fluctuations about the mean field, even under non-equilibrium conditions. Although the mean-field equations can be nonlinear in the state variables, the governing equation for the fluctuations is always a linear system of partial differential equations. Thus, solution methods are rapid and can often be obtained exactly.

In the next section, spectral transforms are used to obtain steady solutions to equation (2.22) in the case of purely diffusive transport, ** L**(

**r**)=

**D**∇

^{2}. We emphasize that the theory can be applied to more general linear transport operators, including anomalous diffusion and advection through a velocity field that itself may be governed by highly nonlinear equations.

### (b) Spectrum of the steady-state fluctuations

Pattern formation is ubiquitous in physical, chemical and biological systems. One mechanism through which they can arise is due to the evolution of a deterministically unstable steady state. For example, Turing showed that for a reaction–diffusion model tending to a homogeneous equilibrium state, diffusion can act to destabilize the steady solution (Turing 1953). Moreover, the system becomes destabilized to a narrow range of spatial modes, leading to the emergence of regular patterning.

The Turing space is the range of parameter space for which reaction kinetics are stable at equilibrium, but nonetheless, diffusion drives small perturbations to form spatial patterns. In Fourier space, a Turing instability in the deterministic case leads to a Dirac delta at the wavenumber of the pattern. Noise-induced pattern formation arises in a stochastic system for parameters outside of the Turing space. In analogy to the deterministic case, these patterns are revealed as peaks in the spatial correlation spectrum (Butler & Goldenfeld 2009; Biancalani *et al.* 2010). This particular mechanism of noise-induced spatial patterning is analogous to the ‘resonant amplification’ (McKane & Newman 2005) mechanism underlying temporal instabilities of spatially homogeneous stochastic models, and can be quantified through a local analysis about the equilibrium point. A peak in the noise spectrum corresponds to stabilized transients induced by the fluctuations, and maintained by the slow relaxation back to equilibrium (Kuske *et al.* 2007). Instabilities arising from excitable dynamics require a more global analysis (Kuske *et al.* 2007), and are not amenable to the methods described in this article. The spatial correlation spectrum is related to the Fourier transform of the factorial cumulants, ** Ξ**. Equation (2.22) governs the factorial cumulants for all time; however, the linearity of the moment equations allows particularly convenient evaluation of the steady-state fluctuations in the system. If the deterministic system approaches a spatially homogeneous steady state, then

**(**

*Ξ***r**

_{1},

**r**

_{2}) becomes a function of spatial separation

**(**

*Ξ***r**

_{1},

**r**

_{2}) →

**(**

*Ξ***r**), where

**r**=|

**r**

_{1}−

**r**

_{2}|. Furthermore, the factorial cumulant becomes a symmetric matrix, and the spatial dependence disappears from the coefficients, i.e.

**A**(

**r**

_{1})=

**A**(

**r**

_{2})=

**A**, reducing equation (2.22) to2.23

The presence of the delta function and the linearity of the constant-coefficient partial differential equation yield a simple expression for the Fourier transform of the factorial cumulant .

## 3. Applications

In the following three subsections, we use several examples to illustrate the applicability of the method derived in the previous section. The first application is to simple diffusive transport. By explicitly computing the time-dependent mean and standard deviation, we justify the use of the second factorial cumulant to quantify the magnitude of the fluctuations. The second example is a nonlinear reaction–diffusion model with a homogeneous steady state (Koch & Meinhardt 1994). We show that there are regions of parameter space for which the stochastic model exhibits noise-induced pattern formation. Finally, we consider the fluctuations about a spatially varying steady state in a recent model of regulatory RNA-mediated patterning in early development (Levine *et al.* 2007). We show that the proposed mechanism exhibits very robust and sharply delineated patterns, despite substantial intrinsic fluctuations.

### (a) Pure diffusive transport

Diffusive transport is well described as a stochastic process, exemplified by Smoluchowski’s model of Brownian motion (von Smoluchowski 1906). Beginning with a random walk on a lattice, the equation describing the average number density reduces to the diffusion equation in the continuum limit,3.1

From the same master equation for the random walk, it is straightforward to show that the factorial cumulant obeys the two-dimensional analogue of equation (3.1) (van Kampen 1976*b*, 1992),3.2

Using equation (2.15), the factorial cumulant is used to compute the variance of the density about the mean-field state by solving equation (3.2) with absorbing (Dirchlet) boundary conditions (figure 2). For a model with one spatial dimension, it is straightforward to discretize the spatial domain into subvolumes, and perform a stochastic simulation using the Gillespie direct method (Gillespie 1977), where the diffusive transport is encoded as an elementary reaction event (as depicted in figure 1; Elf & Ehrenberg 2004).

At the steady state, the second factorial cumulant vanishes (consistent with the one-step diffusion process converging to a Poisson distribution), with the relative magnitude scaling like the square root of the average density. It is important to note that, because the transport is linear, equation (3.2) is exact. In some recent work on intrinsic fluctuations in spatially varying systems (Butler & Goldenfeld 2009), diffusive transport is treated deterministically, with the result that the magnitude of the fluctuations is underestimated.

### (b) Activator–inhibitor model

In the case of a spatially homogeneous steady state, our general formulation recovers the analysis of Biancalani *et al.* (2010) in their dilute limit. As in that work, we observed noise-induced spatial patterning in a subset of reactants. A simple network that can be used to illustrate the effect of intrinsic fluctuations on pattern-forming instabilities is the activator–inhibitor model of Gierer & Meinhardt (1972). Rescaling time and space leads to a minimal toy model given by the following system governing the concentration of activator *A* and inhibitor *H* (Koch & Meinhardt 1994),3.3The reaction part of the dynamics admits a single steady state (denoted by the superscript ‘ss’),whose stability depends upon the two control parameters *σ*_{A} and *ρ*_{H}, with asymptotic stability guaranteed for3.4The model is temporally unstable for *ρ*_{H}<*ρ*_{c}, which is depicted as the hatched region in figure 3. The condition for the appearance of Turing instabilities is,3.5The Turing space for this model as a function of the inhibitor degradation rate *ρ*_{H} is depicted as the solid black region in figure 3. As the reaction kinetics become more stable (*ρ*_{H}≫*ρ*_{c}), a larger disparity in the diffusivity *D*_{H} is required to generate pattern-forming Turing instabilities. Notice that if the steady-state *A*^{ss}(*x*) exhibits stable pattern formation, then the synthesis term in the dynamics of *H* in equation (3.3) necessarily becomes spatially dependent. In that way, a Turing instability in only one reactant is impossible for the deterministic case.

The stochastic analogue of the deterministic system, equation (3.3), requires specification of the stoichiometry matrix, **S**, and the propensity vector, ** ν** (equation (2.1)). For simplicity, we assume the stoichiometry matrix consists of unit steps. Enforcing periodic boundary conditions in space, with subvolumes of unit length, figure 4 illustrates the behaviour of the model for the parameter choices indicated in figure 3. For parameters in the region of deterministic instability, the stochastic simulation shows regular temporal oscillations displayed as vertical bands in figure 4

*a*. In contrast to a simulation of the deterministic model (not shown), the peak heights are non-uniform across the spatial domain, and there is slight irregularity in the period of oscillation. Figure 4

*b*shows the result of a simulation in the Turing space of the deterministic model. Here, the stochastic simulation shows regular spatial patterning displayed as horizontal bands. As in figure 4

*a*, there is some irregularity in the wavelength of the pattern, and movement of the edges in time. In all, for parameters chosen from deterministically unstable regions of the phase plot, the stochastic simulation and the deterministic model are in qualitative agreement (so long as the fluctuations remain subdominant to the mean-field behaviour). Of interest are those parameter choices for which the behaviour of the stochastic and deterministic models are no longer in agreement.

Close to the Turing space of the deterministic model, the steady state is only weakly stable (figure 3, dark grey region). The action of the fluctuations brought about by the discrete reaction events and diffusive transport are enough to destabilize the system and allow spatial patterns to form (figure 4*c*). Compared with the spatial patterns in figure 4*b*, the boundaries are far more ragged and the wavelength of the pattern is distributed over a range of values. Nevertheless, there is obvious pattern formation in a parameter regime where the deterministic model is asymptotically stable. Figure 4*d* shows an example simulation from a region of the parameter space demonstrating a surprising difference between the stochastic and deterministic models (figure 3, light grey region). Here, the spatial pattern is less distinct than in figure 4*c*, though still observable. What is remarkable is that patterning in the activator is not accompanied by a patterning of the inhibitor, a feature that is more clearly seen in the spectra of the spatial correlation functions. Note, however, that the temporal stability of these spatial patterns requires analysis of the full spatio-temporal correlation.

At steady state, the Wiener–Khinchin theorem (van Kampen 1992; Gardiner 2004) can be invoked to relate the Fourier transform of the simulation data to the covariance spectrum of the fluctuations *S*_{ii}(*k*),3.6As a consequence, the analytical value of **S**(*k*), rather than , obtained from the linear noise approximation will be compared with the value computed in the numerical simulations. Spatial patterning is evident in the spectrum of the spatial correlation function as a peak at non-zero wavenumber *k* (Butler & Goldenfeld 2009; Biancalani *et al.* 2010). Figure 5*a* shows the spectra of the activator and inhibitor corresponding to the simulation data shown in figure 4*c*. There is a narrow peak in both the activator and inhibitor spectra close to *k*=0.6, corresponding to the spatial pattern of wavelength approximately 10. The activator spectrum in figure 5*b* likewise shows a peak close to *k*=0.6, but more broad than in figure 5*a*. Notice, however, there is no discernible peak in the inhibitor spectrum; i.e. there is no finite *k*>0 for which we find a local maximum in the spectrum characterized by d*S*_{HH}/d*k*=0. The asymmetry between the stability of the activator and the inhibitor arises from the positive feedback loop in activator synthesis—consequently, we would expect one-species patterning to be a generic consequence of autoactivator–inhibitor models. Although this particular example is a toy model, it clearly demonstrates the qualitative and quantitative differences that can arise owing to intrinsic fluctuations in spatially varying systems; effects that can arise in other more complex models. Furthermore, for sufficiently large positive feedback, it may be possible for a model to exhibit noise-induced spatial structure without requiring an associated disparity in the transport coefficients.

It is useful to highlight difficulties that may arise in the comparison of the continuous approximation to stochastic simulation data (figure 6). The first type of error (illustrated in figure 6*a*) arises from low spatial resolution in the local density estimation of the discrete simulation data. When converting individual-based data into a continuous density field, local number counts must be made in a finite-sized subvolume. It is important to choose a subvolume that is small enough on the scale of spatial structuring in the model, so that spurious spatial averaging is not incurred. This problem applies equally to lattice-based simulation methods, like a discretized Gillespie simulation as employed here, or off-lattice methods (van Zon & Rein ten Wolde 2005). The second difficulty occurs at low number densities, where ensembles must be large enough that information in the moments and spectra are not washed-out by sampling effects (figure 6*b*), and higher order corrections to the approximation equations may be required. As for the approximation method itself, since linear transport fluctuations are treated exactly, the limits of applicability are the same as in the linear noise approximation for spatially invariant systems. Depending upon the strength of the feedback in the reaction network, and other nonlinearities, the approximation may still apply to very small populations of the order of a few dozen molecules. Nevertheless, for systems with only a handful of reactants, for example in fluorescence correlation spectroscopy (Qian 1990) or enzymatic reactions in small volumes (Grima 2009), it may be necessary to include higher order corrections to the first two moment equations, leading to a renormalization of the mean and, significantly, a second relaxation time in the time-correlation function (van Kampen 1976*a*). These important extensions of the present methodology will be discussed in more detail elsewhere.

### (c) Small regulatory RNA in developmental patterning

Levine *et al.* (2007) have recently proposed a mechanism to sharpen spatial mRNA expression profiles using small regulatory RNA. Denoting the small regulatory RNA by *μ*, and the mRNA target by *m*, the model assumes that the small regulatory RNA alone diffuses, and that the interaction between the species results in mutual annihilation,3.7Here, *β*_{m} and *β*_{μ} are the linear degradation rates of mRNA and regulatory RNA, respectively, *κ* is the RNA interaction parameter and *D* is the diffusion coefficient for the regulatory RNA. The spatially varying synthesis rates, *α*_{m}(*x*) and *α*_{μ}(*x*), are assumed to be anti-correlated sigmoidal functions,3.8where the spatial domain length *L* has been scaled so that 0≤*x*≤1. Under conditions of strong interaction, *κα*_{m}(0)/*β*_{m}≫*D*/*L*^{2}, solving equation (3.7) subject to reflecting (Neumann) boundary conditions results in an mRNA profile that exhibits a very sharp transition from high-to-low expression states (figure 7, solid black).

Levine *et al*. demonstrated that the mean behaviour is robust to intrinsic fluctuations, but it remains of interest to ask what effect fluctuations have on the variance of the interface position over an ensemble of realizations (i.e. how do fluctuations affect the accuracy of the interface position). Using equation (2.22), the system of partial differential equations governing the factorial cumulants ** Ξ** is simply3.9where3.10and ° is the component-wise Hadamard product defined in equation (2.21). Substituting the mean-field solution, equation (3.7), into the coupled two-dimensional Poisson equations for

**, the system is solved with absorbing (Dirichlet) boundary conditions. From the definition of the factorial cumulant, equation (2.15), the standard deviation about the averaged state is then computed (figure 7, dashed lines).**

*Ξ*There are several features of interest regarding the fluctuation estimates in figure 7. First, despite the spatially varying steady state, the standard-deviation estimates correspond well to the stochastic simulation traces, even with very small number densities (figure 7*b,c*). Second, from the standard-deviation envelope, it is clear that the mechanism is able to provide a sharp threshold with a location that is robust to intrinsic fluctuations—with less than 4 per cent relative error in the threshold location even down to densities of a few molecules. This remarkable accuracy of the threshold location is a feature of the mechanism that is not apparent from the deterministic model, and would require a large number of stochastic simulation runs to establish numerically.

## 4. Discussion

Stochastic models of physical systems often capture realistic features that are not possible in a deterministic framework (McKane & Newman 2004; Fange & Elf 2006). A stochastic formalism allows for more accurate representation of spatially varying models, particularly in the regime of low concentrations where the continuum assumption breaks down. This higher fidelity comes at a cost: stochastic models tend to be very difficult to analyse analytically, and stochastic simulation becomes the only means to gain insight into the model dynamics. For spatially varying systems, stochastic simulations are computationally expensive and can be prohibitively time-consuming in higher dimensions, so that a systematic mapping of parameter space is no longer possible. Here, we have derived an analytic perturbation method to approximate the evolution of the first and second moments of a stochastic spatially varying system governed by the reaction-transport master equation. In contrast to past work, the method can accommodate non-equilibrium (time-dependent) spatial structure and spatially inhomogeneous steady states. One advantage of our method is that the equations for the second moments are necessarily linear. As a consequence, for models with a homogeneous steady state, the equilibrium spatial correlation spectrum is simply determined by a linear algebraic system of equations (equation (2.23)). In that case, peaks in the spatial correlation spectrum indicate noise-induced spatial patterning, and facilitate the identification of regions of interest in the parameter space of the model, expediting stochastic simulation and providing a broad overview of the available model dynamics. In all examples discussed, the results of the approximation compare favourably with stochastic simulation. We note that the accuracy of the approximation depends upon the local density of the reactants being sufficiently large. In cases when that assumption breaks down, for example in fluorescence correlation spectroscopy (Qian 1990) or small-vesicle dynamics (Grima 2009), higher order corrections must be included in the reaction contribution to the first- and second-moment equations. Likewise, nonlinear transport requires specialized treatment (Lugo & McKane 2008; Biancalani *et al.* 2010).

The method is derived in full generality, and can be applied to arbitrary numbers of species, in arbitrary spatial dimension, transported through the action of a linear operator that can be mapped to a master equation (for example spatially dependent advection or anomalous diffusion). Of great interest is the application of these methods to constrain parameter estimates from spatially varying experimental data, as has been done in the spatially invariant case (Munsky *et al.* 2009). Finally, while exact stochastic simulation is prohibitively resource-intensive for spatially varying systems in three dimensions, hybrid accelerated methods could potentially be developed using the methods derived herein to separate fast-local reactions from slow-disperse reactions.

## Acknowledgements

The authors are grateful to Peter McHale and Marek Stastna for comments on the manuscript, and NSERC for financial support (to F.J.P.).

- Received May 28, 2010.
- Accepted August 3, 2010.

- © 2010 The Royal Society