The Linear Noise Approximation for Spatially Dependent Biochemical Networks

An algorithm for computing the linear noise approximation (LNA) of the reaction–diffusion master equation (RDME) is developed and tested. The RDME is often used as a model for biochemical reaction networks. The LNA is derived for a general discretization of the spatial domain of the problem. If M is the number of chemical species in the network and N is the number of nodes in the discretization in space, then the computational work to determine approximations of the mean and the covariances of the probability distributions is proportional to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M^2N^2$$\end{document}M2N2 in a straightforward implementation. In our LNA algorithm, the work is proportional to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M^2N$$\end{document}M2N. Since N usually is larger than M, this is a significant reduction. The accuracy of the approximation in the algorithm is estimated analytically and evaluated in numerical experiments.


Introduction
Many biochemical networks are modeled by ordinary or partial differential equations at a macroscopic level of fidelity. Such continuous models may not be sufficiently accurate when the number of molecules involved in the chemical reactions is small. This is often the case in molecular cell biology (Elowitz et al. 2002;McAdams and B Per Lötstedt perl@it.uu.se Arkin 1997; Raj and van Oudenaarden 2008;Tsimring 2014). Chemical reactions are then best described as random events, and the discrete number of molecules is important when the copy numbers are low at a mesoscopic level of modeling. The macroscopic equation for the mean values is often satisfactory when the number of molecules is large. Analytical solutions to the governing macroscopic or mesoscopic equations can be obtained only for special systems. Computational methods are needed for quantitative information about the behavior of the systems.
The master equation (ME) or Kolmogorov forward equation is an equation for the time evolution of the probability density function (PDF) for the copy numbers of different species in systems with an intrinsic noise (Gardiner 2004;Gillespie 1992;van Kampen 2004). The systems are modeled as Markov processes with discrete states defined by the copy numbers of the chemical species in continuous time. The particular ME for spatially homogeneous, well-stirred problems in chemistry is the chemical master equation (CME) where reactions between two molecules occur with a propensity that depends on the copy numbers of the species. The ME is generalized in the reaction-diffusion master equation (RDME) to spatially heterogeneous chemical systems by introducing a discretization of the reaction volume into compartments or voxels (Gillespie et al. 2014;Gillespie and Seitaridou 2013). The state is then given by the copy numbers in each one of the voxels.
The computational work and the storage requirements to solve the RDME grows exponentially in the number of species and the number of voxels making the simulation of biochemical systems with the ME prohibitive except for very small systems. Analytical solutions are known only for limited classes of problems such as those with linear propensities. Instead, sample trajectories of well-stirred systems are generated by Gillespie's stochastic simulation algorithm (SSA) (Gillespie 1976(Gillespie , 1977. The original algorithm has been improved in many ways, e.g., for efficiency (Cao et al. 2006;Gibson and Bruck 2000;Gillespie 2001) and for systems with slow and fast reactions (Cao et al. 2005;E et al. 2007). The Gillespie algorithm is generalized to problems with spatial variation due to diffusion in Elf and Ehrenberg (2004), Engblom et al. (2009), Isaacson and Peskin (2006), Lampoudi et al. (2009) and implemented in software (Drawert et al. 2012(Drawert et al. , 2016Hattne et al. 2005). The computational effort may be quite large to simulate a system with many chemical species, many molecules, and many voxels, since many realizations of the process are required in a Monte Carlo method like SSA due to slow convergence to the mean and other moments of the distribution. An introduction and an overview of Markov models of chemical reactions are found in Goutsias and Jenkinson (2013). Recent reviews of computational methods at different levels of modeling are , , Mahmutovic et al. (2012) and Sokolowski et al. (2017).
There are ways to approximate the solutions to the CME with deterministic equations. The linear noise approximation (LNA) is obtained from the CME by deriving the equations for the moments and then expanding the solution in a large parameter , representing the size of the chemical system (van Kampen 1976Kampen , 2004. The means and covariances in LNA are exact for chemical systems with at most first-order reactions where the propensities are constants or linear in the copy numbers. The first and second moments are exact also for other systems with a special structure (Grima 2015). Different modifications have been proposed to improve the accuracy of LNA, see, e.g., Ferm et al. (2008), Grima (2012). Some of the improvements are compared experimentally in examples in Schnoerr et al. (2015). The LNA and similar approximations are used to quickly study biochemical networks in, e.g., Elf and Ehrenberg (2003), Thomas et al. (2013), Ullah and Wolkenhauer (2009) and more recently as a surrogate model to infer parameters in biochemical models from data in Fearnhead et al. (2014), Fröhlich et al. (2016, Ruttor andOpper (2009), Stathopoulos andGirolami (2013). A review of LNA and related methods and their use for inference are found in Schnoerr et al. (2017).
An alternative to the LNA is the EMRE approximation in Grima (2010) extended to spatial problems in Smith et al. (2017). The covariances satisfy the same Lyapunov equation as we derive here. The spatial EMRE algorithm is applied to gene regulation in a cell and to reactions in a aggregation of cells in two space dimensions in Smith et al. (2017). More compartments than one are also found in Challenger et al. (2012). The equations of LNA with spatial variation are derived in Scott et al. (2010) and applied to the modeling of spatial patterns. The equation for the covariances is replaced by an equation for the factorial cumulant. Equations similar to the LNA for spatial problems are used in Butler andGoldenfeld (2009), Anna et al. (2010) to investigate oscillatory systems. Turing patterns are studied in Asllani et al. (2013) Diffusive effects are important for the fidelity of models when the chemical reactions are localized in space in a cell and when the molecular transport is slow compared to the reactions. Some examples where the spatial effects are crucial are found in Fange and Elf (2006), Sturrock et al. (2013), Takahashi et al. (2010). The LNA is a level of modeling suitable for such systems, e.g., to infer parameters for the diffusion and the reactions from measurements, at least in the beginning of the iterative search process for the parameters, thanks to the relative simplicity of LNA.
In this paper, we develop a fast algorithm for computing approximations of the mean and the covariance of the PDF solving the RDME based on the LNA for spatial problems with reactions and diffusion. The equation for the expected values is a system of reaction-diffusion equations, and the equation for the covariances is a time-dependent Lyapunov equation with a source term localized in space. Let M be the number of chemical species and N the number of voxels. The structure of the covariance equations is utilized to compute an approximation of the covariance and to reduce the computational work and the memory requirements from being proportional to M 2 N 2 in a straightforward implementation to M 2 N in our algorithm. Since N usually is larger than M, this is a substantial reduction. A bound on the deviation of the true covariance from our approximation is proved in a theorem. The accuracy of the covariance approximation is demonstrated in numerical examples in one, two, and three dimensions (1D, 2D, and 3D).
In the next section, the RDME is given and a splitting of the operator is introduced. The equations of the LNA for spatially heterogeneous chemical systems are derived for general shapes of the voxels in Sect. 3. A continuous approximation of the equation for the covariances is analyzed in Sect. 4. The algorithm is presented in Sect. 5 for computation of the mean and the covariance. Numerical results are found in Sect. 6. Finally, some conclusions are drawn.
The notation in the paper is as follows. The ith element of a vector v is denoted by v i . The jth column of an array x with elements x i j is written as x · j , and x i· is the ith row. The derivative of v i (x) with respect to x j is denoted by v i, j . The time derivative ∂ p/∂t of p(x, t) is written as ∂ t p, andq is a shorter notation for dq/dt. The Euclidean vector norm is denoted by v = i v 2 i and the subordinate spectral norm for a matrix A is A . The set of integer numbers is written as Z, and Z + denotes the nonnegative integer numbers. In the same manner, R denotes the real numbers and R + is the nonnegative real numbers.

The Master Equation
Consider a biochemical system with M chemically active species. The system evolves on a line (1D), in an area (2D), or a volume (3D) V which is partitioned into N voxels (or compartments) V j such that they cover V, V = N j=1 V j , and are non-overlapping, V j V k = ∅, j = k. The voxels are defined by a computational mesh constructed for numerical discretization of partial differential equations, see Fig. 1. The size of a voxel is V i = |V i |, and the diagonal matrix V has the elements V i in the diagonal. Each voxel has a node in the center with coordinates x ∈ R d , d = 1, 2, 3 and the nodes are connected by edges. The molecular copy number of species i in voxel j is a random integer Y i j . The state of the system is time dependent and is given by y(t) which is an array of nonnegative integers, y ∈ Z M×N + . The state changes randomly with reactions between the species in a voxel and with diffusive jumps of the molecules between the voxels.
The CME is a Kolmogorov forward equation for the PDF p(y, t) for a system to be in the state y at time t (Gardiner 2004;van Kampen 2004). The state changes at discrete time points after a chemical reaction in a voxel. Ifỹ is the state before the Meshes with edges (solid), nodes j, k, and l, and voxel V j (dashed). The nodes j and k are connected by edge e jk . a An unstructured mesh, b a structured Cartesian mesh reaction and y is the state immediately after the reaction r , then the change in state in V j can be written asỹ (1) The reaction occurs with the propensity w r , i.e., with probability w r t in a short time interval t. The state change vector n r ∈ Z M tells how the state is updated after a reaction. Most entries of n r are zero and n ri = 0 only for those species involved in the reaction. In a system with R different reactions, the CME for p(y, t) is defining the master operator M. Diffusion of the molecules is modeled as jumps between voxels with a common boundary. Suppose that V j and V k share a point in 1D, an edge in 2D, or a facet in 3D. Then, a molecule of species i in V j jumps to V k with propensity q jkỹi j The probability for a molecule to jump is given by the jump coefficient q jk . The state change vector has two nonzero components: n jk, j = −1, n jk,k = 1. The diffusion master equation (DME) in a chemical system without reactions is defining the diffusion operator D. The RDME for the PDF of the state of a system with reactions and diffusion is then The jump coefficients q jk are determined by the geometry of V and the voxels V j and the diffusion coefficient γ . If V is a rectangle in 2D or a rectangular hexahedron in 3D and the voxels are squares or cubes, then the mesh partitioning V is Cartesian as in Fig. 1b. When V j and V k share a boundary and the size of an edge in the mesh is x, we have q jk = γ / x 2 . If V j V k = ∅, then q jk = 0. For a general shape of V, the voxels are defined by an unstructured mesh consisting of triangles (2D) as in Fig. 1a or tetrahedra (3D) in Engblom et al. (2009). If there is an edge e i j between node j in V j and node k in V k , then q jk > 0.
The diffusion equation for u(x, t) with Neumann boundary conditions in d dimensions is where ∂V is the boundary of V and n is the normal of ∂V. The equation is discretized in space by a finite element method in Engblom et al. (2009) to derive q k j . Let u i j be the concentration of y i j in V j such that u i j = y i j /V j . Then, u i· for one species i satisfies after discretizatioṅ where J ( j) is the set of nodes connected to node j by an edge. It is shown in Engblom et al. (2009) that With a mesh of good quality, S jk ≥ 0 for j = k, see Meinecke et al. (2016). The relation between the jump coefficients in (3) and the diffusion coefficients in (7) is To simplify the notation, the assumption here is that the diffusion speed is equal for all species. Otherwise, q k j , D jk , and S jk in (9) would be scaled by γ i /γ to account for the different diffusion coefficients γ i of the different species i.
A numerical solution of (5) is seldom computationally feasible due to the high dimension of y ∈ Z M N + . Suppose that the lattice for y has L points in each coordinate direction, i.e., y i j ∈ [0, 1, . . . , L]. Then, the lattice size for y is L M N . A simplification is possible by first splitting the operator M + D into two parts (Hellander et al. 2014;MacNamara and Strang 2016;Strang 1968). Suppose that the solution is known at t n . A timestep t is chosen, and then, the reaction part (2) is integrated from t n to t n+1 = t n + t followed by integration of the diffusion part (4) in The splitting error in p at t n+1 is of O( t) in (10). Second-order accuracy is obtained by evaluating the reaction equation at half the step at t n+1/2 = t n + 0.5 t, then solving the diffusion Eq. (4) for a full step, and finally taking half a step with (2) (Strang 1968) The solution has been advanced from t n to t n+1 . The error in p( In (2), the reactions occur independently in every voxel without being influenced by the species in the other voxels. Introduce the ansatz into (2) to arrive at Hence, for t ≥ t n and p(y · j , t n ) given, N separate solutions of the CME can be computed and then combined in (12) in the first step of (10) (or the first and third steps in (11)). The solution is computed N times on a lattice of size L M . This is smaller than the lattice for the full problem but may still be prohibitively large for numerical solution.
In the same manner, the species diffuse independently of each other in the second step in (10) and (11). Insert into (4) and rearrange the terms as in (13) to arrive at M separate equations for the diffusion of the species when t ≥ t n and p(y i· , t n ) is known Since the propensity is linear in y in (16), there exist analytical solutions to the subproblems with multinomial and Poisson distributions for p, see Jahnke and Huisinga (2007), but in practice they are not so useful due to the size L N of the lattice.

Linear Noise Approximation
The biochemical systems are assumed to have a scaling with a size parameter as in van Kampen (2004), Kurtz (1970), Kurtz (1971) where 1. In chemical applications, can denote a volume or a typical copy number of the species. Then, the copy numbers are rescaled by , z = −1 y, and the propensities can be written as Equations for approximation of the PDF of a system are derived below. The computational complexity of their solution is polynomial in M and N instead of exponential as in (5), (14), and (16).

The Mean Value Equation
Let m i j be the expected value E[Y i j ] of the copy number Y i j of species i in voxel j with a PDF satisfying the master Eq. (2). Multiply (2) by y i j and sum over Z N + . Then, m i j satisfies the equation (Ferm et al. 2008;van Kampen 2004) Suppose that every w r is linear in y. Then, and Eq. (18) are exact for the mean values. If w r is nonlinear in y · j , then an approximation is With this approximation, we obtain the reaction rate equationṡ The mean value equations scaled by the size parameter μ i j = m i j / arė The mean concentration

The Linear Noise Approximation
The scaled state variable z ·k in voxel k is split into a deterministic part μ ·k and a random part η by van Kampen in van Kampen (1976Kampen ( , 2004 for the chemical reactions. The random term is assumed to be proportional to −1/2 . The relation between the copy numbers y ·k , the scaled copy numbers z ·k , the fluctuations η, and the fluctuations in This expansion is inserted into master Eq. (2) with the propensities v r in (17) Replace p in (23) by in and expand the right-hand side of (23) in a Taylor series around μ ·k . Terms proportional to 1/2 vanish since μ ·k satisfies (20). If is the solution to then terms of O(1) cancel out. Terms of O( −1/2 ) and smaller are ignored in the expansion. This is the linear noise approximation (LNA) for the scaled copy numbers subject to chemical reactions in V k . The solution to (25) is the PDF of a normal distribution see Ferm et al. (2008), van Kampen (2004) and (Risken 1996, p. 156). The dimension of η is M in the CMEs (14) for all N voxels. The matrix for the covariance between the species i and j in V k is the solution oḟ Since η is normally distributed with 0 mean and covariance , η ∼ N (0, ), it follows from (22) that Y ·k , Z ·k , and the concentration U ·k = V −1 k Z ·k also have normal distributions The covariance of U ·k in (28) . Then, the differential equation satisfied by follows from (27) There are M nonlinear ODEs to solve in (20) for μ ·k in every voxel V k . The covariance matrix is symmetric, and we have to solve (M + 1)M/2 linear ODEs in (27) and (29). The structure of this equation is the same also for other approximations of the CME, e.g., EMRE in Smith et al. (2017). The accuracy of mean value Eq. (20) is improved in Ferm et al. (2008) by adding a term which is linear in the covariance.

The Diffusion Equation
The notation is simplified if we assume here that there is only one species, M = 1, but many voxels, N > 1. If M > 1, then the diffusion of the other species is treated separately in the same manner, see (16). In diffusion master Eq. (16), the propensity to jump from voxel k to j is linear in y ∈ Z N + with w r (y) = q k j y k and v r (z) = q k j z k . The linearity implies that there are explicit expressions for the mean value equations, ν i,k in (27), and W i j in (25).
The equations for the scaled mean values are obtained from (8), (9), and (20) The diffusion equation for the mean concentration is derived from (21) cf. (7). The equation for the covariance (27) between voxels i and j depends on ν j,k and W i j in (25). The derivative ν i,k in (27) and (29) is by (8) and (9) since n ri = 1 for the jump from k to i and n ri = −1 for all jumps from i to every j connected to node i by an edge e i j . Let E be the set of all edges in the mesh. The state change vector on edge e i j for a jump from i to j is n i j with the nonzero components n i j;i = −1 and n i j; j = 1. The contribution to W in (27) from e i j jumps in two directions: i → j and j → i. Hence, for all edges The nonzero elements of N i j = n i j n T i j are N i j;ii = N i j; j j = 1 and N i j;i j = N i j; ji = −1. Therefore, the elements of the symmetric W are The random component of the concentrations ψ = V −1 η is normally distributed The coefficients in (35) multiplying are In (35), W is scaled by V Then, the scaled W-term in (37) can be written in a symmetric form The covariance equation corresponding to (29) for diffusion is by (35), (36), and The copy numbers Y k· and Z k· are normally distributed as in (28). The covariance of the concentrations in space of a species k, U k· , is −1 and U k· ∼ N (u k· , −1 ).
In the stationary equation,˙ i j = 0 in (39) and it is a Lyapunov equation. A stationary solution of (31) is u i = const. and of (39) is where δ i j is the Kronecker delta. If the initial data i j (0) are symmetric, then the solution to (39) is symmetric for all t > 0. At the stationary solution of (40) where diag(x) is a diagonal matrix with x i in the diagonal. Thus, the stationary distribution of the copy numbers in the voxels Y k· for species k follows from (28) The stationary copy numbers in different voxels are uncorrelated, have a multivariate normal distribution and are therefore independent, and are approximately Poisson distributed since Y ki ∼ N ( V i u i , V i u i ) with equal mean and variance. If i j (0) = 0, then the time-dependent solutions to (39) will be proportional to u i /V i and the mean and the covariance of Y i j are both proportional to V i u i . The distributions of the solution to the DME in (4) are discussed in  based on the theory for linear propensities in Jahnke and Huisinga (2007). Their distributions are either multinomial, Poisson, or a combination. The stationary distribution is multinomial according to Anderson et al. (2010) and approximately Poissonian .
The components of the solution u of (31) are the node values of the finite element approximation of u(x, t) solving diffusion Eq. (6) for one species. Let (x 1 , x 2 , t) be the covariance between the solutions at the coordinates x 1 , x 2 ∈ R d . Then, i j in (39) can be interpreted as the value of (x 1 , x 2 , t) at the nodes at x 1i and x 2 j . The coefficient D i j in f i j in (39) is negative when x 1i = x 2 j , and positive when On a regular mesh with a typical length of an edge equal to x, the positive weight D i j depends approximately only on r i j = ξ i j and V i varies smoothly with a typical size V . On such a mesh, and the solution u(x, t) to (6). The scalings σ and ω are chosen such that σ ∝ x and ω = x −1 . When r = 0 then D ii V i ≈ ϕ(0) = −γ / x 2 V and when r = x we have Then, the continuous equation corresponding to discrete Eq. (39) is If exact initial conditions of the distribution of molecules are known, then (x 1 , x 2 , 0) = 0. Let u ∞ (x) be the stationary solution to the diffusion equation. Then, one can show that an approximate stationary solution to (43) is As x → 0, this solution approaches u ∞ (x 1 )δ( x 1 − x 2 ) where δ is the Dirac measure. The solution in (40) to discrete Eq. (39) is similar to (44).

Analysis of the Covariance Equation
A property of the continuous approximation (x 1 , x 2 , t) of the covariance in (43) is derived in this section. We show that decays exponentially when x 1 − x 2 grows, indicating that the discrete variance i j in (39) Consider (43) in free space x 1 , x 2 ∈ R d and for t ≥ 0 with initial data (x 1 , x 2 , 0). The concentration u(x 1 , t) is nonnegative and is assumed to be bounded by C u for all x 1 and t ≥ 0. Introduce a change of variables The diffusion equation in (43) is in the new variables Here f (ξ 1 , ξ 2 , t) = u(x 1 , t) + u(x 2 , t) is nonnegative and bounded by 2C u . The factor ϕ in the source term vanishes quickly when ξ 1 increases.
With the fundamental solution of the diffusion equation in 2d dimensions (Evans 2010;Stakgold 2000), the solution of (46) can be written as a sum of two integrals depending on the initial data and the source term where and The integral with the source term is bounded by The spatial integral I d2 over ζ 2 ∈ R d in (50) is The integral I d1 of the product over ζ 1 ∈ R d in (50) is With τ = 4γ (t − s) and α = τ −1 + σ −2 , we have Using (53), (52), and (51), a bound on I src in (50) is Assume that the initial data are localized close to ξ 1 = 0 such that | (x 1 , x 2 , 0)| ≤ 0 exp(− ξ 1 2 /χ 2 ) for some χ > 0. A bound on the integral in (47) due to the initial data is then Hence, a bound on the covariance solution in (47) is obtained by (54) and (55). The assumptions and conclusions are summarized in a theorem: Theorem 1 Assume that |u(x, t)| ≤ C u and that the initial data satisfy The relations between the x and ξ coordinates are Then, the solution of (46) with ϕ defined by (42) for t > 0 is bounded by The function f d depends on the dimension d and is 0 at t = 0. The first term in (56) is proportional to C u /V = sup u/V since σ ∝ x. The solution in (56) decays exponentially in ξ 1 = 1 √ 2 x 1 − x 2 for a fixed t in all dimensions and is small when ξ 1 > σ ∝ x. For a given ξ 1 , the first term in (56) increases slowly with t in 1 and 2 dimensions and is bounded by C u σ 2 / x 2 V (d − 2) when d ≥ 3. The second term in (56) decreases when t ≥ 0.5 ξ 1 2 − 0.25χ 2 for d = 1 and for all t ≥ 0 when d ≥ 2.
Our bounded domain V for x 1 and x 2 has a boundary that is not taken into account in (56). The bound on is a good estimate when the main part of the solution is away from the boundary, e.g., when t is not too large and u(x, t) is nonzero only in the middle of V.
Since (43) is a continuous approximation of (39) we expect the discrete variance i j to behave in a similar way and be negligible when the nodes i and j are not neighbors and not directly connected by an edge in the mesh. This property will be exploited in the algorithm in the next section.

Algorithm
The algorithm to compute the solution to the LNA for both reactions and diffusion is based on the operator splitting in Sect. 2, the derivations in Sects. 3.2 and 3.3, and Theorem 1 in Sect. 4.
The mean value equation in (21) is added to the diffusion equation in (31) to obtain the reaction-diffusion equation for the concentration u ik of species i in voxel k with M N componentṡ The covariance between the concentrations of the species i and j in voxels k and l is written as i j;kl and has M 2 N 2 components. The equation satisfied by the covariance is obtained from (29) and (39) The reaction source term vanishes if the concentrations are from different voxels, since there is no reaction between molecules in separate voxels. The diffusion source term is zero if the species are different since a diffusion event occurs when the same species changes location between adjacent voxels by a jump.
The equations for the mean and the covariance (58) and (59) are solved by splitting the operator on the right-hand side and advancing the solution one timestep from t n to t n+1 = t n + t as in (10) with u ik (t n ) and i j;kl (t n ) as initial data, cf. (10). The algorithm is Algorithm 1 The discretization error in u ik (t n+1 ) and i j;kl (t n+1 ) will be of O( t). The ODEs in steps 1, 2 and 3 update u and in a voxel (steps 1, 2) and between two adjacent voxels (step 3) due to the reactions as in step 1 of (10) and (14). In the ODEs in steps 4, 5 and 6, u and change due to diffusion between voxels without any influence of the other species as in step 2 of (10) and (16). A more accurate splitting than in Algorithm 1 with an error of O( t 2 ) is possible in the same manner as in (11).
It follows from Theorem 1 that if i j;kl (t n ) decays rapidly when the nodes x k and x l are separated then this property is preserved in˜ i j;kl (t n+1 ) where C u > 0 in (56) in step 2 and C u = 0 without the source term in step 3. Using the same arguments in steps 5 and 6, we find that if i j;kl (t n ) decays rapidly when x k − x l increases, then after one timestep i j;kl (t n+1 ) also decays rapidly in x k − x l . Supported by the analysis in Sect. 4, we assume that i j;kl is negligible when node l and node k are not neighbors, l / ∈ J (k), and we letˆ i j;kl = 0 in a sparse approximation of i j;kl . Then, only elements of i j;kl when k = l and l ∈ J (k) need to be stored and updated inˆ i j;kl by Algorithm 1. The sparsity (or nonzero) pattern ofˆ i j;kl for each pair i, j is the same as that of S and D in (8) since S kl and D kl are nonzero only on the diagonal and if nodes the k and l are neighbors connected by an edge in the mesh and l ∈ J (k). Moreover, i j;kl is symmetric in both i and j and k and l. With M different species and N voxels, i j;kl in general has 1 2 M N (M N + 1) different elements butˆ i j;kl has only C d M 2 N nonzero elements that are necessary to store taking the symmetry into account. The coefficient C d depends on the dimension and the structure of the mesh. In a Cartesian mesh, C d = 2(1D), 3(2D), or 4(3D) and in an unstructured mesh C d = 2 in 1D but C d depends on the particular mesh in 2D and 3D. The mean value vectors u andũ have the dimension M N .
In order to estimate the computational work in the steps of the algorithm, we assume that ν i depends on a limited number of u jk independent of M. Then, there are also a limited number of derivatives ν i, j different from zero and independent of M. Thus, the work to compute the right-hand side (RHS) in step 1 in (60) is independent of M and N and it is computed once for every species i and voxel k, i.e., M N times. Since there are a limited number of nonzeros in ν i,α , the sums and W in step 2 in (61) are computed independently of M and N . Hence, the work is proportional to M 2 N for the covariances between the species in every voxel. In step 3 in (62), M 2 covariances are computed for every combination of voxels k and l where˜ is nonzero. This is the case when k and l are neighbors and each k has a limited number of neighbors. This number is independent of N . Therefore, the work to compute the full RHS in step 3 is of O (M 2 N ). The number of D kβ = 0 in the RHS in step 4 in (63) is independent of N according to the previous paragraph. The work to determine all derivatives of u ik is then proportional to M N . For ii;kl to be nonzero in step 5 in (64), voxels k and l are neighbors, l ∈ J (k). Furthermore, the products in the sums are nonzero only if β ∈ J (k) ∩ J (l). The work to calculate the sums is independent of N , and the RHS is computed O(M N ) times. In the same manner, the RHS in step 6 in (65) is computed O(M 2 N ) times. The conclusion is that the work to determine the RHS in the ODEs for u ik andˆ i j;kl in the algorithm has linear complexity in N and is proportional to Since there are additional administrative costs in Algorithm 1, the straightforward algorithm ignoring the sparsity of will be faster when N < N * for some small N * which is problem dependent. However, for N > N * Algorithm 1 will be the winner and its advantage is greater, the greater the N is.
If the diffusion coefficient is different for different species i, then D kβ and D lβ in steps 4, 5 and 6 would depend on i but the algorithm and its properties remain the same.
In summary, the algorithm in words is for one timestep t from t n to t n+1 = t n + t: 1. Solve the ODE in (60) numerically for the mean values with initial data u(t n ) to obtainũ(t) 2. Solve the ODE in (61) numerically for the covariances between the species in the same voxel k withũ from step 1 and initial data (t n ) to determine˜ ··;kk (t) 3. Solve the ODE in (62) numerically for the covariances between the species in different voxels k and l satisfying l ∈ J (k) withũ from step 1 and initial data (t n ) to determine˜ ··;kl (t) 4. Solve the ODE in (63) numerically for the mean values with initial dataũ(t n+1 ) from (60) to obtain u(t) and u(t n+1 ) 5. Solve the ODE in (64) numerically for the covariances between voxels k and l satisfying l ∈ J (k) for the same species i with u from step 4 and initial datã (t n+1 ) from steps 2 and 3 to determine ii;·· (t n+1 ) 6. Solve the ODE in (65) numerically for the covariances between voxels k and l satisfying l ∈ J (k) for different species i and j with u from step 4 and initial datã (t n+1 ) from steps 2 and 3 to determine i j;·· (t n+1 ) In the first three steps, the mean values and the covariances change due to the reactions and in the last three steps due to the diffusion.
Theorem 1 and numerical experiments in Sect. 6.1 indicate that the accuracy in increases when the dimension grows. By storing and updating only the sparse approximation in steps 3, 5 and 6 in Algorithm 1, considerable savings are possible in computing time and computer memory when N is large, e.g., in 3D.

Example
Consider the reversible reaction for association and dissociation of the species A, B, and C with copy numbers μ T · j = (a j , b j , c j ) in voxel j and propensities and state change vectors The macroscopic reaction coefficients are k a = V k k a and k d = k d . Then, Eq. (21) for the concentrations in step 1 of the above algorithm in V k iṡ Order the means and the covariances such that The Jacobian J of the propensities with J i j = ν i, j in (31) and the source term in step 2 in V k are since n 1 n T 1 = n 2 n T 2 . Introduce K and G using the identity matrix I N of size N , J, and W in (70) , where ⊗ denotes the Kronecker product. Then, the equation in steps 2 and 3 in matrix form is˙ The submatrix D is the approximation of the Laplacian in (31) and (39). If the diffusion varies between the species, then D in (73) would be replaced by γ i /γ D, i = 1, 2, 3, on the diagonal. The sparsity or nonzero pattern in H i;·· is the same as in D. The diffusion equation for the mean values in step 4 in Algorithm 1 is as in (31) The matrix form of steps 5 and 6 in Algorithm 1 is (cf. (38) and (39) In 1D, D is a tridiagonal matrix and if V = x is constant then

Numerical Results
The algorithm is tested for computing the mean and the approximation of the covariance in the LNA of systems with diffusion in 1D, 2D, and 3D and a system in 2D with the reversible reaction (66).

Diffusion
A Cartesian grid in d dimensions is generated with a constant step size x and a diffusion coefficient γ = 0.01. The number of dimensions is d = 1, 2, 3, and the domain is the unit cube [0, 1] d . The number of grid points is n = 1/ x + 1 in each dimension yielding N = n d components in u. A straightforward implementation of Algorithm 1 in steps 5 and 6 will generate N 2 elements in . By updating only those elements of that correspond to nonzeros in D and S, the number of nonzero elements in the approximationˆ will be of O(N ).
The initial data u(0) are sampled from a uniform distribution u k (0) ∼ U[0, 1] and = 0. The ODEs in (74) and (75) are solved numerically for t ≥ 0 by the forward Euler method for simplicity. Then, the RHS in each step of Algorithm 1 is evaluated once requiring a computational work proportional to M 2 N in every timestep from t n to t n+1 . Better numerical accuracy is achieved by splitting the computations according to Strang (1968) as in (11) and by using a higher-order method. Better numerical stability is obtained by an implicit method.

1D
The covariance (x 1 , x 2 , t) is computed in 1D on a grid with x = 0.025 and N = n = 41 using the full without zeros, the sparseˆ with the same nonzero pattern as S, i.e., the diagonal, the subdiagonal, and the superdiagonal are nonzero in a tridiagonal matrix as proposed in Sect. 6.1, and the sparseˇ where another two diagonals below and above the diagonal are nonzero in a pentadiagonal matrix. One row of the three matrices is shown in Fig. 2. The approximationsˆ (0.5, x 2 , t) anď (0.5, x 2 , t) agree fairly well with (0.5, x 2 , t) in particular for larger t in Fig. 2. The PDF of the multivariate normal distribution N (u, ) is The covariance matrix is factorized by = Q Q T where Q is orthogonal and has the positive eigenvalues of on the diagonal. The expression in (76) in the exponential defines surfaces of ellipsoids in R N with equal probability, and the eigenvalues of are the lengths of the principal axes of the ellipsoids. Another way of comparing and its approximations is then to compare the eigenvalues to see the difference between the lengths of these axes.
The eigenvalues of ,ˆ , andˇ are displayed in Fig. 3. They also agree well except for the one or two smallest ones in the figure. Using five diagonals and four neighbors inˇ improves the approximation somewhat compared toˆ . Including more than the nearest neighbors in 2D and 3D with an unstructured mesh is possible but makes the algorithm more complicated.

2D and 3D
In 2D, (x 1 , x 2 , t) is computed with the full matrix and with the approximationˆ that has the same sparsity pattern as S on a grid with x = 0.05 and N = n 2 = 441. One row of corresponds to one coordinate x 1k and its covariance with the 2D x 2 . The variance is high at (x 1 , x 1 , t), and the covariance (x 1 , x 2 , t) is very low when x 1 = x 2 . This is depicted in the left panel of Fig. 4 where x 1 and x 22 are fixed and x 21 varies in x T 2 = (x 21 , x 22 ). The differences in covariance between andˆ are very small and not visible in the figure. Since is symmetric, the result is similar in other directions in x 2 and for other x 1 . The steady-state solution (40) with u = 0.5 and V = 1/400 is here 200.
One section of the 3D covariances andˆ is shown in the right panel of Fig. 4. As in 2D, x 1 , x 22 , and x 23 in x T 2 = (x 21 , x 22 , x 23 ) are fixed and (x 1 , x 2 , t) is plotted as a function of x 21 . The step size in the grid is x = 0.0833 and N = n 3 = 2197. The scaled difference between the covariances andˆ is defined by The dominant elements inˆ are the variances on the diagonal. With small elements in compared to 1, the difference between the covariances in andˆ is small relative to the variances inˆ . In Fig. 5, for the 2D example is shown at two time points. The values of are low in blue color in most parts of the matrix. The peaks in the left panel are at 0.035 in isolated points. In the right panel, max i j < 0.02. The eigenvalues of andˆ in 2D and 3D are compared in Fig. 6. The sparse approximation captures all the eigenvalues except for one or two of the smallest ones.
The covariance of the fluctuations in concentration between different parts of the domain is well approximated by the sparseˆ , especially in 2D and 3D in Figs. 4 and 6. This is expected from Theorem 1 in Sect. 4 where the decay of (x 1 , x 2 , t) is slower in 1D than in 2D and 3D when ξ 1 = 1 √ 2 (x 1 − x 2 ) is growing.

Reactions and Diffusion in 2D
The time evolution of the chemical reaction (66) on the Cartesian mesh in 2D in Sect. 6.1.2 is computed with the LNA as in the example in Sect. 5.1. The parameters are k a = k d = 0.1, and the diffusion is low with γ = 0.01. The dimension of u is M N = Mn 2 = 1323, and the initial values in u(0) are uniformly distributed between 0 and 1 and (0) = 0. The eigenvalues of the covariances with the full matrix and with the sparse matrixˆ are compared at different t in Fig. 7. An approximate stationary solution u ∞ is determined at t = 20. The convergence of the three subvectors u 1· , u 2· , and u 3· in (69) corresponding to the concentrations of A, B, and C is displayed in the lower right panel in the figure. The difference between u and u ∞ is measured in · for the species. In the resolution of the figure, it is not possible to distinguish between the differences in convergence between the species. The balance equation k aāb = k dc is satisfied with a relative error less than 0.008 by the mean valuesā,b, andc of the components in u i· (20), i = 1, 2, 3. At ∞, u ∞1· , u ∞2· , and u ∞3· are constant in space.
The convergence plot in Fig. 7 shows that the variation in the solution is larger for small t and decreases with t. The covariances andˆ agree well for large eigenvalues for small t and they agree well for all eigenvalues when t grows. The off-diagonal submatrices i j;·· , i = j, in (69) are comparable in size to the diagonal submatrices ii;·· when t is small but as t grows ii;·· , i = 1, 2, 3, will dominate and be closer and closer to diagonal matrices. There is a jump in the spectrum for larger t, e.g., at t = 6. This is explained by the difference in the size of the stationary valuesā,b, and c whereā ≈b, andc/ā ≈ 0.6. The approximation in the covariance equation behaves in the same way with reactions as in Sect. 6.1.2 without the reactions.

Conclusions
The master equation is a model for biochemical reactions and diffusion but the numerical solution of it is impossible except for simple, well-stirred systems with special properties. An alternative for large systems with spatial variation is to use the linear noise approximation (LNA). We have derived the equations for the LNA for diffusion and chemical reactions on general meshes. The reactions involve M species, and the mesh consists of N voxels. The covariance of the concentrations is approximated by a sparse representation in an algorithm such that the computational complexity is reduced from O(M 2 N 2 ) in a straightforward implementation to O(M 2 N ) here. Also the memory to store the solution is reduced in the same way. The approximation is supported by analytical expressions showing that the higher the dimension is, the better the approximation is. Consequently, the quality of the approximation and the savings in work and storage are more prominent in 3D when N is large. The accuracy of the approximation is evaluated by comparing the elements and the eigenvalues of the full covariance matrix and its sparse approximation in numerical examples with only diffusion in 1D, 2D, and 3D and an example in 2D with a reversible reaction and slow diffusion.