Local error estimates for adaptive simulation of the reaction–diffusion master equation via operator splitting

Abstract

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.

Stochastic Simulation Service: Bridging the Gap between the Computational Expert and the Biologist

“…On the stochastic well-mixed level, where the distribution of molecules is assumed to be spatially homogeneous, StochSS offers the Stochastic Simulation Algorithm (SSA) [16, 17] (commonly known as the Gillespie method) as well as adaptive time-stepping tau-leaping [18], as implemented in the StochKit2 [13] stochastic simulation package for well-mixed systems. For spatial stochastic simulation in up to 3 dimensions, StochSS offers highly efficient implementations of the Next Subvolume Method (NSM) [19] and the Adaptive Diffusive Finite State Projection (ADFSP) algorithm [20, 21], as implemented in the PyURDME [22] spatial stochastic simulation package. The meshes are unstructured, and the spatial geometries can be complicated and may include multiple domains, for example membrane and cytoplasm.…”

Stochastic simulation of reaction-diffusion systems: A fluctuating-hydrodynamics approach

“…The FHD approach has been also applied to concentration fluctuations in a ternary liquid mixture in equilibrium [68] and the Model H equations for binary mixtures [69].The key difference between the FHD and RDME descriptions lies in the more efficient treatment of fast diffusion. A number of approximate numerical methods for the RDME [42][43][44][45][46][47] are based on operator splitting using first-order Lie or second-order Strang splitting [70]. In Appendix A we review and discuss in more detail a split scheme that uses multinomial diffusion sampling [71] for diffusion and SSA for reactions.…”

Accurate reaction-diffusion operator splitting on tetrahedral meshes for parallel stochastic molecular simulations

“…In our method we want to keep τ at the upper limit calculated for diffusion, allowing multiple reactions per update, but also to minimize loss of accuracy by considering how the reactions affect diffusion probabilities locally. This differs from previous methods where a half-step method is applied and tau adapted to a tolerance of errors calculated from the half-step application 16 . When allowing multiple reactions to occur between diffusion updates, the problem becomes how to calculate the number of molecules to diffuse to a neighbor at τ, given that the population may have changed during τ by reactions.…”

“…However, there appears to be a cost to error estimation and it is shown that systems must be stiff for there to be any performance gain, with the algorithm in fact performing slower than the ISSA implementation known as the Next-Subvolume Method (NSM) for non-stiff, low molecule systems. Since the 'SSA' operation also contained a diffusive term, making the initial algorithm inherently serial, the algorithm was further developed in 2014 16 for parallelization. This algorithm is based on a two half step Lie-Trotter splitting implementation with an adaptive time-step controlled by local error estimates based on the half-step method, and an initial parallel version applied up to 4 cores shows good scaling promise.…”

“…For arguments that we have presented, the best way to achieve this is to separate reaction from diffusion and deal with them separately, similarly to previous approaches 11,16,21 .…”