{'fraction_non_alphanumeric': 0.03314917127071823, 'fraction_numerical': 0.008832740964915367, 'mean_word_length': 4.4129523809523805, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Local error estimates for adaptive simulation of the Reaction-Diffusion Master Equation via operator splitting

The efficiency of exact simulation methods for the reaction-diffusion master equation (RDME) is severely limited by the large number of diffusion events if the mesh is fine or if diffusion constants are large. Furthermore, inherent properties of exact kinetic-Monte Carlo simulation methods limit the efficiency of parallel implementations. Several approximate and hybrid methods have appeared that enable more efficient simulation of the RDME. A common feature to most of them is that they rely on splitting the system into its reaction and diffusion parts and updating them sequentially over a discrete timestep. This use of operator splitting enables more efficient simulation but it comes at the price of a temporal discretization error that depends on the size of the timestep. So far, existing methods have not attempted to estimate or control this error in a systematic manner. This makes the solvers hard to use for practitioners since they must guess an appropriate timestep. It also makes the solvers potentially less efficient than if the timesteps are adapted to control the error. Here, we derive estimates of the local error and propose a strategy to adaptively select the timestep when the RDME is simulated via a first order operator splitting. While the strategy is general and applicable to a wide range of approximate and hybrid methods, we exemplify it here by extending a previously published approximate method, the Diffusive Finite-State Projection (DFSP) method, to incorporate temporal adaptivity.


Introduction
To understand biological systems on the cellular level, it is often essential to account for the impact of noise due to small molecule count. This was beautifully demonstrated in the now classic results on the effect of stochasticity in gene regulatory systems [23,8]. Spatial distribution of molecules in a cell can result in locally small populations of key chemical species, such that noise drives essential behavior, as in the case of symmetry breaking across many eukaryotic cell types [30]. Spatial stochastic modeling has already begun to yield new insights in systems such as spatiotemporal oscillators [10,27,26], MAPK signaling [28], self-organization of proteins into clusters [7] and polarization of proteins on the cell membrane in yeast [1].
On a mesoscopic level, the reaction-diffusion master equation (RDME) is frequently used to model these systems. Space is subdivided into subvolumes that can individually be treated as well-mixed. Reactions within a subvolume are expressed in the form of the chemical master equation (CME) [13] and realizations of the process can be generated using Gillespie's stochastic simulation algorithm (SSA) [12]. Molecules can move freely between neighboring voxels via diffusive jumps, which are modeled as linear jump events in a Markov process. Optimized exact simulation methods such as the Next Subvolume Method (NSM) [7] can be used to generate statistically correct realizations of the RDME.
As with all exact methods applied to RDME models, the NSM suffers from a potentially high computational cost due to having to explicitly simulate each diffusion event. The number of diffusive transfers between voxels grows rapidly as the mesh resolution is made finer, and as a result the majority of computation time tends to be spent on sampling diffusion events. Additionally, these methods are also inherently serial, which has thwarted attempts to increase efficiency via parallelization.
To speed up simulation of the RDME many methods rely on operator splitting. By splitting the operators, most often with a Lie-Trotter scheme [29], the reaction and diffusion steps can be solved independently. While diffusion carries the bulk of the computational burden in exact solvers, it is often possible to take advantage of the structure and linear nature of the discretized diffusion equation to speed up this step in an operator-split solver. Examples of approximate methods that have been proposed to speed up the simulation of the RDME by reducing the cost of the diffusive step include methods based on tau-leaping [25,21], the multinomial simulation algorithm [19], spatially adaptive hybrid methods [11] and the diffusive finite state projection method (DFSP) [6]. While splitting in itself does not resolve the issue of the inherent stiffness of the diffusion operator, the continued introduction of methods for simulating the RDME via operator splitting highlights the potential of this approach. Another recent use of operator splitting in the RDME context is the use of Lie-Trotter splitting to simulate fractional diffusion [2]. Yet another advantage of splitting is that it converts a largely serial problem, which is known to be difficult to parallelize in an efficient manner, into a naturally parallelizable one. For the existing approximate and hybrid methods for the RDME, splitting the physics (reaction and diffusion) is necessary. For parallel implementations, another possibility is to split the computational grid into blocks, as proposed by [14], where the error introduced by operator splitting at block boundaries was analyzed. Our analysis here is different since it applies to the case of splitting the reaction and diffusion operators.
Splitting the operators introduces an error that depends on the size of the splitting time step, however previous algorithms that rely on operator splitting have not attempted an a priori error estimator. Without such an estimator these methods have no way to automatically control the splitting error. This limits their usefulness for practitioners. From an efficiency point of view, not knowing and controlling the error might lead to the use of unnecessary small timesteps at the price of slower simulations.
In this work we seek to address theses issues by presenting estimators of the local error in probability, mean and variance for a first order Lie-Trotter splitting of the RDME. The estimators allow control of the splitting error for spatial stochastic simulation and enables any method based on operator splitting to be implemented adaptively.
This paper is organized as follows: in Section 2 we introduce the RDME. In Section 3 we outline how to simulate the RDME using operator splitting. We derive our estimator for the local error in the PDF, mean and variance, and demonstrate the accuracy of the estimator for an example problem. Finally, in Section 4 we present how the local error estimates can be used to extend an approximate method for the RDME to incorporate temporal adaptivity.
2 Spatial Stochastic Simulation using the RDME Given a system with N s chemical species X s reacting in a physical domain Ω, discretize Ω into N v non-overlapping voxels V i , with volume |V i |, and let the state of the system be described by the (N v × N s ) state matrix x, where the element x is is the copy number of species X s in voxel V i . Let x i,· denote the i-th row of x and x ·,s the s-th column. The reaction network consists of N r chemical reactions r = 1 . . . N r . The (1 × N s ) stoichiometry vector n ir describes the change in the state,x i,· = x i,· +n ir , when reaction r occurs in voxel i and the propensity function for that reaction is a ir (x i,· ).
Diffusion of species X s along the edge (2D) or face (3D) connecting voxels V i and V j is modeled as a linear jump event, or diffusive transfer, with propensity function µ(x) = d ijs x is . The change in state due to the diffusive transfer is described by the (N v × 1) stoichiometry vector ν ijs , such that the new state isx ·,s = x ·,s + ν ijs . All entries of ν ijk are zero except ν ijs (i) = −1 and ν ijs (j) = 1. In the case of reactions only, the probability density function p(x, t) ≡ p(x, t|x 0 , 0) obeys the master equation For the case of one subvolume, N v = 1, (2) reduces to the CME for a well stirred system. For a system with only diffusion, the master equation takes the form For a system with both diffusion and reactions, the full RDME is then given by The values of the diffusion rate constants d ijs depend on the diffusion constant γ s of species X s and the shapes and sizes of voxels V i and V j . Let u s (ζ, t) = E[X s /|Ω|], i.e. the the expected value of the concentration process corresponding to (3). In the thermodynamic limit X s /|Ω| → ∞, u s is governed by the diffusion equation Here, ζ is a spatial coordinate in a Cartesian coordinate system. Similarly, the Fokker-Planck equation that describes the time evolution of the probability density of a single particle undergoing Brownian motion is given by Hence, by choosing d ijs according to a consistent spatial discretization of the operator ∆, we obtain mesoscopic jump coefficients that are motivated both from a microscale and macroscale perspective. To comply with the description of the RDME, the discretization used needs to couple only nearest neighbors in the mesh. For a uniform, Cartesian mesh such as shown in Figure 1 (a), the most natural choice is a centered finite difference discretization, giving d ijs = γ s /h 2 where h is the side length of the voxels and γ s is the diffusion constant of species s. For unstructured, tetrahedral meshes which will be used later in this paper, d ijs can be obtained from a finite element (FE) discretization using linear Lagrange elements. For a detailed description of that methodology, see [9]. Figure 1 (b) shows a triangular, unstructured mesh. The local volume where the molecules are assumed to be well mixed are given by the dual mesh. In the conventional RDME, molecules are allowed to jump between immediate neighboring voxels.
3 Operator splitting method for the RDME In this section we describe how an operator splitting method can be applied to generate realizations from the RDME with a controlled temporal discretization error.
In the remainder of this paper we will tacitly assume that we can impose a bound on the state space. For every species in the system we assume that P (X is > x max , t) = 0, for all t. The state space is then finite (but very large). Technically, this will not necessarily hold true for an open system, but from a biophysical perspective it is a reasonable assumption since a system would require infinite energy to blow up. In the finite case, the operators in the RDME can be represented by finite matrices. Hereafter we will use the notation M and D to mean representations of the bounded operators resulting from the above truncation of the state space. The solution of (4) can then simply be written We point out that even though this representation of the solution is simple, it is not feasible to solve (4) this way since the state state space is too large except for trivial models and discretizations. In the well mixed case, however, deterministic methods to solve the CME have been developed for small to medium-sized models, see for example [24,3,20,18,22]. Using a first order splitting method (Lie-Trotter splitting), p(x, t + ∆t) can be approximated by p s (x, t + ∆t) = e ∆tM e ∆tD p(x, t).
Simulation based on (8) proceeds in two half-steps, with the diffusion operator and the reaction operator acting sequentially Provided that the time step ∆t is sufficiently small, the local error in the scheme (9) is O(∆t 2 ) and from standard theory for numerical solution of differential equations, the global error in the PDF is then proportional to ∆t, i.e.
We note that it is not in general necessary for the operators to be bounded for the splitting method to result in a global error proportional to ∆t. Jahnke shows in [16] that the global error for Strang splitting applied to the CME is second order under certain assumptions on the chemical reactions. Unfortunately, as Jahnke points out, while those conditions can be expected to hold for many systems, they are not easily verified in practice. Error bounds for exponential operator splitting have also been studied in other contexts in e.g. [17].
In the next subsection we derive estimates of the local error. In a subsequent section we will then illustrate how these estimates can be used to enable temporal adaptivity in an approximate method for the RDME.

Local error
Following standard theory for the analysis of scheme (9), the local error in the PDF (i.e. in probability) in the n+1-th timestep is given by where [D, M] = (DM − MD) is the commutator of the operators. This follows directly from a series expansion of the left hand side in ∆t. Computing n+1 from (11) is obviously not tractable in general since knowledge of p n requires the solution of the full RDME. Even if an approximation of p n were available, the state space is still very large and the cost of evaluating the commutator would be prohibitive. Also, our goal is a strategy to estimate the error during the course of the generation of individual trajectories, whereas an estimation of p n would require very large ensembles of trajectories. By conditioning on the currently observed state x n in timestep n of a specific realization, we obtain the following sequence of conditional errors, where δ is the Dirac delta function, hence δ(x n ) corresponds to P (X(t n ) = x n ) = 1.
Note that E n is a random vector for each n with unconditional expectation E[E n+1 ] = n+1 .

6
While there are many possible states that can have a non-zero value after one application of the commutator, hereafter called reachable states, it is possible to obtain simple and computable expressions for E n+1 in (12), due to the sparsity of δ(x n ). The x n + n kr and x n + n kr + ν ijs , s = 1 . . . N s , r = 1 . . . N r , i, j, k = 1, . . . , N v and N e is the number of connections in the mesh. Below, the term to the left of the colon is a reachable state and the term to the right is the value of that state after applying the indicated operator to δ(x n ).
For example, the application of the diffusion operator D on δ(x n ), Dδ(x n ) results in where i, j are connected subvolumes and d 0 ( Nv j =i d ijs x is . Equation (13) enumerates the new states with non-zero value after one application of D to δ(x n ).
Applying M to Dδ(x n ) then gives the following reachable states and associated values.
x n + ν ijs + n kr : where we have defined a 0 (x N ) = Nv k=1 Nr r=1 a kr (x). By computing the analogous values for D(M(x n )δ(x n )) and taking the difference, we obtain the following expressions for the commutator error E n+1 for the reachable statesx l x n + n kr + ν ijs : x n + n ir : − x n + ν ijs : x n : 0.
As seen in the expression for the cross-terms in (19), large parts of the operators commute. Only terms pertaining to degrees of freedom that are sharing an edge in the mesh will be non-zero. Furthermore, (19) will be non-zero only if the reaction described by n kr depends on species s, i.e. n krs is non-zero. A graphical representation of the reachable states is given in Fig. 2.
Local error in mean and variance. Equations (19) - (22) give expressions for the local error in the PDF in one timestep (at t n+1 ), given that we observe state x n at time t n . The obvious advantage of computing these expressions directly is that they can be used to derive many estimates of the error such as l 1 , l ∞ in PDF, Kolmogorov distance for marginal distributions, or moments. However, they are rather expensive to compute directly, even with appropriate optimizations. In many cases it is sufficient to control a weak error such as the errors in mean and variance. This is particularly true if the error in the individual methods for propagating the reaction operator and the diffusion operator is controlled only in a weak sense. For example, this would be the case if the diffusion operator were updated with τ -leaping such as in [11].
Based on (19) -(22) we can calculate the local error in the expected value of an arbitrary bounded function g(X is ) where the subscriptf denotes expectation under the approximate PDF obtained from solving with the operator split method (9). For brevity, we have here dropped the superscript n + 1 on x and E. We introduce the notation ∆a − irs = a ir (x i,· + ν ijs ) − a ir (x i,· ) and ∆a + irs = a ir (x i,· + ν jis ) − a ir (x i,· ) and ∆g(y) = g(x n is + y) − g(x n is ). By summing up outflow and inflow of differences in probability to state X is and using the fact that ∆E[g( n jrs a jr (x j,· ) x n x n +n ir +ν ijs x n +ν ijs x n x n +n ir +ν ijs x n +n ir Figure 2: Illustration of the reachable state in Eq. (19) for k = i. On the left side, the diffusion operator is applied first, as in Eq. (14), followed by the reaction operator, resulting in the value derived in Eq. (18). On the right side the order of the operators is reversed, resulting in a change in the final value for the reachable state x n + n kr + ν ijs (note that the reachable state is, by definition, the same in both orders, but the value associated with that state is different). The terms in the box represent the difference between the two orders of applying the operators and show terms that do not commute. Note that if we had defined the reachable state as k = i, j in Eq. (19) then the reaction would occur in neither the originating subvolume nor the destination subvolume for the diffusion event and the terms would commute (and the value in the box would be zero).
To obtain the error in the mean, set g(x) = x. After simplification, we obtain where the operator R(x) is defined by and The (N v × N v ) matrix D s is defined to have diagonal elements d ii = − j =i d ijs and off-diagonal elements d ij = d ijs . Equation (25) is an exact (up to O(∆t 3 )) formula for the estimate of the local error in mean for any functional form of the propensity functions a r (x). For zeroth-order mass action reactions we have σ irs = 0. For a first order reaction with a linear propensity (such as creation or monomolecular conversion) of the form a ir (x i,· ) = k 1 x is we have in agreement with the analogous expression for the commutator error for the reactiondiffusion PDE.
To find the error in the second moment, set g(x) = x 2 in (24), to obtain While this expression is more complicated than (25), the amount of additional work required to compute it is not that great since the complexity of the extra terms are all O(N ) and most of the expensive propensity function evaluations can be overlapped with the computations involved in (25). A more commonly used second order statistic is the variance . While this quantity is not readily obtained from (24), an approximate formula for V [X] can be found from the observation where E f [X] is the expected value under the operator split PDF. Hence In the asymptotic regime, O(∆t 3 ) terms are small by definition and ∆V [X is ] is a good approximation to the true local error in variance. We point out that (11), and hence (25) and (29) are only good estimates of the local error if the O(∆t 3 ) terms are small. It is well known that an error estimate based on (11) will deteriorate in quality for large timesteps if D or M are stiff. As the norm of D increases rapidly with finer mesh resolution, this will lead to overly conservative estimates for large ∆t >> h 2 /γ s (in the non-assymptotic regime). For this reason, error estimates that perform better in the non-asymptotic regime have been devised in the PDE case [4]. However, it is not clear how those approaches would apply to the RDME. In the following section we will demonstrate that (11) is simple enough to lead to a strategy that is both computable in practice during the coarse of a single realization of RDME and local in space so that it has potential to be efficiently implemented in parallel.
To illustrate the accuracy and correctness of the estimate of the local splitting errors (25) and (29), we simulated a model of Min oscillations in E. Coli [7] in one spatial dimension for a single timestep. To isolate the splitting error, the reaction and diffusion steps were solved exactly by simulating them with SSA and NSM respectively. Fig. 3 shows the estimated conditional error in mean (25) and the true conditional error in mean as a function of the spatial coordinate (left) and for different timesteps (right). Fig. 4 shows the corresponding estimate of error in second moment based on (29). As can be seen, the estimated error accurately captures the true local error, and the quality of the estimate improves for smaller timesteps as expected from the O(∆t 3 ) error in the estimates. When computing the error in mean, the observed error is a combination of sampling error caused by a finite number of realizations, and the error caused by operator splitting. For small splitting errors, a large number of realizations are necessary to distinguish the sampling error from the splitting error, especially if we want to measure the error in e.g. the L 1 norm, in which case the variances in the different voxels add up. We note that to achieve tight confidence intervals on the true error for this simple system (see to the top row of Figs. 3 and 4) required 10 11 realizations.  To illustrate the use of the local error estimate, we extended a previously published approximate method for the RDME, the Diffusive Finite State Projection (DFSP) method [6] with temporal adaptivity. This results in increased robustness of the solver. DFSP relies on Lie-Trotter splitting to separate the reaction and diffusion updates over a discrete timestep ∆t. The reactions are then updated in each voxel using Gillespie's direct method [12]. Diffusion is updated by sampling from probability density functions that are precomputed by solving the diffusion master equation (3) locally up to the given timestep ∆t. For an unstructured mesh, there is one such PDF for each vertex in the mesh and for each ∆t. The (spatial) locality of the PDFs are enforced by applying an absorbing boundary condition at a certain distance away from the vertices. The timestep is assumed to be small enough that the majority of probability is located close to the vertex where the molecule is located at the beginning of the timestep so that the error due to the truncation of the statespace is small. It is shown in [6] that this strategy can speed up simulations by effectively aggregating the effects of many fast diffusive transfers over the splitting timestep, and in [5] it is discussed under what conditions one can expect DFSP to be more efficient than NSM for the MinCDE model [15].
DFSP has been implemented previously as an add-on solver in the URDME framework [5]. URDME is a modular framework that uses unstructured meshes for spatial stochastic simulations. URDME includes interfaces for handling of the geometry, mesh and computation of diffusion jump rates for the unstructured mesh. It has a modular design which facilitates the implementation of new algorithms as add-on solvers. In the current implementation of DFSP in URDME, ∆t has to be chosen by trial and error by manually picking a timestep that results in satisfactory error in the computed solutions, judged by a posteriori checks.
We implemented our error estimation strategy in URDME for arbitrary processes (i.e. it uses the same input files as all the core solvers) for structured and unstructured meshes, and supplemented the DFSP solver with an adaptive selection of the timestep. To compute the next proposed timestep, we control the per species error in mean (25) in the L 1 norm such that with s being a user specified relative error tolerance. The L 1 norm is a natural choice in the context of DFSP since the error in the diffusion lookup-tables has a natural bound in the l 1 norm [6]. Figure 5 shows an overlay of the pole-to-pole oscillation pattern of membrane bound MinD in E. coli along with time steps selected by the adaptive step size selection algorithm. The timesteps are themselves stochastic variables and fluctuate during the course of the simulation. In the figure they have been smoothed by a windowed average to more clearly visualize how they adapt to the dynamics of the solution. We simulated the MinCDE model in 3D in three different ways: using the adaptive timestep control, using a fixed manually determined timestep and using the exact NSM solver, and compared the execution times. In DFSP, the diffusion step is conducted by computing new positions for each molecule individually by sampling from a precomputed lookup-table. All the molecules' new positions can then be sampled in parallel. During the reaction update, all the voxels can be treated independently and in parallel. Note that the local error estimators used here have the same characteristics when it comes to the spatial access pattern as the DFSP diffusion operator. Thus the overall adaptive algorithm should parallelize quite well on a shared memory multicore architecture. To demonstrate this, we implemented a naïve parallel version using openMP. Table 1 shows simulation speeds for the adaptive DFSP method, with error estimation conducted every tenth timestep, on a machine with an Intel(R) Core(TM) i7-2600 CPU @ 3.40GHz (Nehalem) processor and 6GB of RAM using one core and using 4 cores (8 threads). For comparison, we also include simulation times of the (serial) NSM solver in URDME. As can be seen, the DFSP algorithm shows an almost ideal speedup on 4 cores, and for this error tolerance, the adaptive code is roughly two times faster than NSM on a single core.
Adaptive DFSP ( s = 0.05) NSM # Voxels 1 core 4 cores (serial)  588  124  35  337  1009  228  65  537  2818  757  190  1282   Table 1: Simulation time(s) for the adaptive version of the DFSP solver and for the exact NSM solver for the MinCDE problem with varying mesh resolution. The trajectories were simulated to a final time of 2000 s, and the state was sampled every second.

Discussion
In this work we have derived local error estimators for first order operator splitting in stochastic reaction-diffusion simulations based on the RDME. Operator splitting provides a way to decouple reactions and diffusion in spatial stochastic simulations, and the local error estimates enable this to be accomplished with a controlled error. The advantage of decoupling the reaction and diffusion operators for simulation are two fold. The first is that decoupling the operators allows for approximate methods to be applied to the diffusion operator to reduce the cost of frequent diffusive transfers. The method we considered as an example in this paper falls into this category. Enabling temporal adaptivity for such approximate and hybrid methods has several benefits. Possibly the most important is the robustness it adds to the solver. While it is possible to prescribe an error tolerance that is likely to give reasonable result across a range of models, the same is not possible by selecting a fixed timestep a priori. Adaptivity has long been available in state-of-the art software for e.g. numerical solution of ODEs and PDEs, but was previously lacking for this class of approximate methods for the RDME. The other major benefit is one of efficiency; using a globally conservative fixed timestep will lead to sub-optimal simulation speeds. We found that the adaptive version of the DFSP method presented here was roughly 1.5 times faster than the corresponding fixed timestep method using a conservative timestep for the oscillatory MinCDE model problem we considered. Even though this experiment is somewhat artificial since the conservative fixed timestep could not be known before running the adaptive code (which again illustrates the utility of adaptivity), it illustrates the efficiency benefits of an adaptive solver. The second advantage of operator splitting in the context of the RDME is that splitting the operators converts the largely serial exact simulation method to one with better potential for parallel implementations. We have here demonstrated that the local error estimates can be efficiently applied also in a parallel implementation of an approximate method for the RDME. The parallelelization aspects of the current work will be explored in greater detail in a forthcoming publication.
Finally, operator splitting allows the flexibility of coupling different types of models (possibly at different scales) and solvers. Thus this error estimator opens avenues of research in the development of hybrid algorithms with controlled errors.