Adaptive mesh refinement for stochastic reaction–diffusion processes

Abstract

We present an algorithm for adaptive mesh refinement applied to mesoscopic stochastic simulations of spatially evolving reaction–diffusion processes. The transition rates for the diffusion process are derived on adaptive, locally refined structured meshes. Convergence of the diffusion process is presented and the fluctuations of the stochastic process are verified. Furthermore, a refinement criterion is proposed for the evolution of the adaptive mesh. The method is validated in simulations of reaction–diffusion processes as described by the Fisher–Kolmogorov and Gray–Scott equations.

Introduction

Simulations of systems that exhibit multiple spatial scales often employ adaptive mesh refinement methods (AMR) [1], [2] for the discretization of the governing partial differential equations (PDEs). However, when the mesh spacing approaches the molecular scale and thermal fluctuations affect the dynamics of the system, the use of partial differential equations may become unjustifiable for variables that are inherently random and discrete [3], [4], [5]. At the same time simulations of microscopic phenomena may not be efficiently computed using atomistic level modeling as in molecular dynamics (MD). There is an increased interest [6], [7], [8], [9], [10], [11], [12], [13], [14] in formulating methods for mesoscopic simulations in which the computational elements contain neither so many particles as to be considered in a continuum, nor so few as to warrant MD simulations. Such mesoscopic simulations employ variables that are assumed to be governed by the laws of probability theory. The transition probabilities of any discrete-state, continuous-time Markov process must obey the Chapman–Kolmogorov equation, which in turn is equivalent to the so-called master equation (M-equation) [15]. Since the number of states is large – possibly infinite – for all but the simplest systems, analytical or direct numerical integration methods for the M-equation are generally impractical. Being at a loss for a direct solution, numerical realizations of the stochastic process can alternatively be generated. These numerical realizations, performed via stochastic simulation algorithms [16], [17], amount to generating random variates from the unknown discrete probability density function.

The M-equation frequently appears in chemical kinetics since the time-evolution of chemically-reacting species in a system is inherently stochastic owing to Brownian motion [18], [19]. The applicability of the M-equation to chemical kinetics rests on the assumption that the system is in thermal equilibrium, such that any molecule has enough time, on average, to move throughout the volume before it participates in a reaction. If Brownian motion does not suffice, then the system is inhomogeneous and the problem has intrinsically local phenomena. The variables in the system are recast to represent the number of molecules of a species in a subvolume – or volume element – of space. This description has been previously studied theoretically [6], [7], [8] and numerically [9] using a uniform discretization in 1-D. Moreover, 2- and 3-D uniform discretizations have been used for simulating nonlinear reaction–diffusion processes [12]. Since inhomogeneities may arise at various scales throughout the domain, the placement of smaller or larger volume elements in certain regions of the domain readily lends itself to finite volume methods [20]. Non-uniform 1-D discretizations, which have borrowed techniques from finite volume methods, have been employed for propagating wavefront systems [10], [11]. Recently, Engblom et al. [13] and Ferm et al. [21] have used adaptive and unstructured meshes in simulations of stochastic processes. Additionally, Drawert et al. [22] have used the Finite State Projection algorithm presented in [23] to simulate reaction–diffusion processes.

This work is concerned with formulating the transition rates, or propensities, for reaction–diffusion processes on adaptive locally refined structured meshes in the spirit of adaptive mesh refinement (AMR) [1], [2], [24], and simulating mesoscopic reaction–diffusion processes using these rates. The convergence of the diffusion process is presented as well as a refinement criterion for the adaptive mesh. The method is applied to relatively high spatial-resolution simulations of the Fisher–Kolmogorov [25] and Gray–Scott [26] equations. The results indicate that the method is especially suited for, but not limited to, wavefront propagation and pattern formation problems.

The paper is organised as follows: Sections 1.1 Master equation, 1.2 Stochastic simulation algorithms introduce the master equation and stochastic simulation algorithms, respectively. Section 2 presents a derivation for diffusion propensities on locally refined meshes, and Section 3 recapitulates reaction propensities and their validity with respect to the spatial discretization. Section 4 presents the mesh refinement criterion and stochastic interpolation method used in this study. Section 5 shows both the validity of the method and numerical examples, and Section 6 concludes this work.

Let $\mathcal{G}\subset {\mathbb{R}}^d$ denote the domain in which stochastic reaction and diffusion processes reside. Following the notation presented in [27], let $\mathcal{G}$ be subdivided into a collection of voxels, or volume elements, labelled by vectors i that belong to an index set, namely $\mathbf{i}\in \mathcal{I}$. Given S species in $\mathcal{G}$, let the time-dependent random variable $U_{\mathbf{i}}^{(s)}(t)$ denote the number of molecules (or particles) of species s in volume element i at time t. Furthermore, define ${\mathbf{U}}_{\mathbf{i}}\triangleq \left(U_{\mathbf{i}}^{(1)}\text{,}\dots \text{,}{U}_{\mathbf{i}}^{(S)}\right)\text{,}{\mathbf{U}}^{(s)}\triangleq {\left\{U_{\mathbf{i}}^{(s)}\right\}}_{\mathbf{i}\in \mathcal{I}}$, and $\mathbf{U}(t)\triangleq \{{\mathbf{U}}_{\mathbf{i}}{\}}_{\mathbf{i}\in \mathcal{I}}$.

The reaction–diffusion master equation (RDME) describes the time evolution of a probability density function:

$$\frac{\mathrm{d}p(\mathbf{U}\text{,}t)}{\mathrm{d}t}=\sum _{\mathbf{i}\in \mathcal{I}}\sum _{\mathbf{j}\in \mathcal{I}}\sum _{s=1}^S\left\{a_{\mathbf{i}\text{,}\mathbf{j}}^{(s)}\left({\mathbf{U}}^{(s)}+1_{\mathbf{i}}^{(s)}\right)p\left(\mathbf{U}+1_{\mathbf{j}}^{(s)}-1_{\mathbf{i}}^{(s)}\text{,}t\right)-a_{\mathbf{j}\text{,}\mathbf{i}}^{(s)}\left({\mathbf{U}}^{(s)}\right)p(\mathbf{U}\text{,}t)\right\}+\sum _{\mathbf{i}\in \mathcal{I}}\sum _{m=1}^M\left\{a_{\mathbf{i}}^{(m)}\left({\mathbf{U}}_{\mathbf{i}}-{\mathit{\nu }}^{(m)}\right)p\left(\mathbf{U}-{\mathit{\nu }}^{(m)}{1}_{\mathbf{i}}\text{,}t\right)-a_{\mathbf{i}}^{(m)}({\mathbf{U}}_{\mathbf{i}})p(\mathbf{U}\text{,}t)\right\}\text{,}$$

where the diffusion process has been modeled as the movement of a molecule to a neighboring volume element [6], [10]; $a_{\mathbf{i}\text{,}\mathbf{j}}^{(s)}(·)$, the so-called propensity function, represents the probability per unit time of the diffusion of a molecule from volume element j to i; a_i^(m)(·) denotes the propensity of reaction m, where m = 1 , … , M; ν^(m) represents the stoichiometric vector of reaction m; and $1_{\mathbf{i}}^{(s)}$ denotes one molecule of species s at index i (i.e. U(t) + 1_i^(s) = U_i^(s)(t) + 1).

Eq. (1) constitutes a set of ordinary differential equations. Note that the evolution of the probability density function is determined by the time-dependent propensity functions $a_{\mathbf{i}\text{,}\mathbf{j}}^{(s)}(·)$ and $a_{\mathbf{i}}^{(m)}(·)$, which represent the unscaled probabilities per unit time of their representative reactions. Since the diffusion events are modeled as unimolecular transitions to neighboring cells, Eq. (1) can be completely characterized by a set of propensities that vary over time. Indeed, stochastic simulations work by generating random variates from the probability density function p(U, t) at every time-step. Stochastic simulations therefore simulate the underlying Markov process, which is defined by the time-dependent propensities.

Given an initial vector of the state space U_sim, stochastic simulation algorithms proceed as follows:

$$\Theta (t)\triangleq ({\hat{a}}_1(t)\text{,}\dots \text{,}{\hat{a}}_W(t))\text{,}$$

$$\tau \sim \zeta (\Theta (t)\text{,}ϵ)\text{,}$$

$$k^{(w)}\sim \Psi (\Theta (t)\text{,}\tau )\text{,}$$

$${\mathbf{U}}_{\mathit{sim}}(t+\tau )={\mathbf{U}}_{\mathit{sim}}(t)+\sum _{w=1}^W{k}^{(w)}{\mathit{\nu }}^{(w)}\text{,}$$

where Θ(t) represents propensities ${\hat{a}}_w(t)$ only from the current state to all other of the W reachable states, where W is the sum of the diffusion and reaction events; the vector ν^(w) denotes the change induced by transformation w; and the distributions ζ and Ψ vary depending on the algorithm that is used (see [18] for a review). These algorithms proceed by iterating through Eqs. (2), (3), (4), (5) until the predefined final time t_f is reached. Ultimately, U_ssa(t_f) is an approximate random variate from p(U, t_f). Specifically, the relationship between the approximate and exact solution is:

$$p(\mathbf{U}\text{,}{t}_f)-p({\mathbf{U}}_{\mathit{sim}}(t_f))\approx C_1ϵ+C_2{N}^{-1/2}+C_3{h}^2\text{,}$$

where C₁ϵ is the error from the integration over time of the stochastic simulation (ϵ = 0 in the case of exact simulation algorithms [17], [28]), C₂N^−1/2 is the well-known sampling error of Monte Carlo methods, and C₃h² is the error from the spatial-discretization of the diffusion process.