Correction factors for boundary diffusion in reaction-diffusion master equations

One of the fundamental goals of computational cell biology is to deepen our understanding in the complicated dynamical processes taking place within cells, thus increasing overall systems predictability and control. In this context, the correct simulation of biological processes may imply simultaneously using different temporal and spatial modelling techniques. For instance, molecules are discrete entities, while their diffusion is continuous, and the chemical reaction events between them are stochastic. The impact that these features have on the dynamics of the biological process under investigation determines the suitability of a modelling approach. Namely, one has to assess if the biological system of interest can be best described as discrete or continuous, homogeneous or heterogeneous, deterministic or stochastic.

Temporal models are employed whenever there are no apparent spatial inhomogeneities, noticeable dependencies or – as is generally the case – whenever spatial data are simply unavailable. In such models, the participating molecules are readily assumed to be well mixed within the system. The standard deterministic model is then given by the reaction rate equations representing the evolution of chemical species as continuously varying concentrations.

On the other hand, the stochastic, mesoscopic description of the system's dynamics is specified by the corresponding chemical master equation (CME) which, given an initial condition, describes the time evolution of the probability P(X, t) for having X = (x₁, …, x_N) molecules in a system at time t, where N is the number of interacting chemical species and x_i is the number of molecules of species i, i ∈ [1, N].

In principle, one can “add” spatial resolution to the CME framework, yielding a reaction-diffusion master equation^1,2 (RDME). Here, the spatial domain of the reaction system is discretized into subvolumes (SV) or “voxels,” usually equally sized, and under the assumption that particles within each subvolume are well mixed. In fact, in 2D domains one is effectively dealing with subareas; however, it has become customary to use the term “subvolume” in 2D scenarios as well. In the RDME, the system's state is then described by the number of molecules per species per subvolume, and first-order reactions are used to model particle diffusion, approximated by jumps between (usually neighboring) subvolumes.

Quite naturally, obtaining the probability function described by a chemical or reaction-diffusion master equation can be computationally challenging, if not unfeasible, due to the combinatorial explosion of possible states of the system. This is one of the main reasons why trajectory-based approaches, commonly known as stochastic simulation algorithms^3–7 (SSAs), are so popular. Their outputs are trajectories that are exact realizations of the CME or RDME, the first moments of which can be approximated by averaging over a large number of individual trajectories. Additionally, the stochastic mean behavior of the reaction-diffusion system can sometimes be approximated by the solution of associated systems of differential equations. To be more precise, for reaction systems solely composed of zero- and first-order reactions, exact analytical solutions of the CME can be derived.^8,9 Since diffusion jumps are modelled as unary reactions, this also applies to RDMEs when the highest reaction order is at most one. Moreover, it has been shown for a linear gene transcription model that the first and second moments of the associated solution match those of the corresponding chemical Langevin equation (CLE).¹⁰ Additionally, the solutions of CLE, CME, and the corresponding system of rate equations (REs) match only for such linear systems, while, in general, this cannot be expected.¹¹ Recently, modifications of the conventional RE formalism have been suggested, many times yielding reasonable approximations of the evolution of the mean of species concentrations (mass-fluctuation kinetics,¹² effective mesoscopic rate equations¹³). However, to the best of our knowledge, none of these results has been explicitly extended to the spatial setting.

The RDME framework entails additional issues that have to be considered when modelling spatial aspects of reaction systems. For instance, molecules are described without a physical dimension. Hence, such models are many times insufficient to describe volume (hard-sphere) exclusion and other effects induced by molecular crowding. In this respect, it is known that molecular crowding prevents reacting molecules from reaching regions of the domain, due to the high concentration of macromolecules impeding their passage.¹⁴ More intuitively, if we consider an arbitrary RDME simulation, a voxel may be populated by a number of molecules that is actually too large to physically fit. Also, the artificial nature of the lattice may not only limit the spatial resolution of the simulation, but also introduce lattice anisotropy.¹⁵ However, one should note that a modified CME has been recently proposed¹⁶ that employs “corrected” propensity functions according to the scaled particle theory of hard-particle fluids. This approach has been shown to predict mean concentrations and intrinsic noise statistics in environments with low to intermediate levels of crowding.

On top of these shortcomings, a critical step in posing a correct RDME model is to choose an appropriate discretization level and, by consequence, upper/lower limits for subvolume sizes in RDMEs have long been discussed. We will review such limits in Sec. II. Moreover, in order to guarantee accurate results, the reaction rates commonly used in RDMEs have to be carefully reassessed when assuming reactive boundary conditions and when second-order (or higher order) reactions are considered, as they are entirely dependent of diffusion. Dissociations, unlike other first-order reactions, are driven by diffusion-dependent separation and, by consequence, also depend on the spatial discretization.¹⁷ 

The main focus of this paper is the correct implementation of RDME models for systems with (semi-)permeable or reactive boundaries in 2D/3D. For this, we will provide a simple approach for deriving discretization and diffusion dependent correction factors of boundary jump rates (Sec. III). In addition, we briefly discuss the correction of bimolecular reactions. We conclude with a summary, while all mathematical derivations and justifications of formulas can be found in the Appendixes. Our methodology represents an alternative approach to current correction strategies and provides straightforward ways to correct readily available RDME simulation algorithms in significantly wider ranges of diffusion settings.

The inherent problem of choosing appropriate discretization limits was already noticed in the mid 1990s. Back then, Baras and Malek-Mansour¹⁸ suggested the reactive mean free path (RMFP) to function as an upper limit for a suitable discretization size, where the RMFP refers to the average distance that a particle travels before undergoing a reactive collision with other moving particles. Since then, it has been argued that the RMFP is a suitable upper limit to ensure homogeneity within each subvolume. Particularly, in hard-sphere collision models of gas molecules, the RMFP is typically calculated as

$1/(2πd2nv)$ where d is the particle diameter, and n_v is the number of molecules per unit volume. More recently, Fange et al.¹⁷ have argued that microscopic re-binding events happen on length scales smaller than the RMFP and, hence, suggest to use subvolumes smaller than the RMFP together with modified, mesoscopic rate constants.

Separately, Bernstein¹⁹ proposed the condition τ_D ≪ τ_C to guarantee homogeneity, where τ_D is the time between diffusion events and τ_C is a typical time between reaction events. It should be noted that τ_D can be approximated by h²/(2dD), where D is the diffusion coefficient of a particle, d is the dimension of the domain (d = 2 or d = 3), and h is the length scale of a subvolume. Thus, in such a view, the discretization size h should be chosen to be small enough such that the above condition is met.

Yet another criteria for choosing a proper spatial discretization size has been given by the Kuramoto length^20–22 defined as

$L Kuramoto S=DSτS$⁠, where D_S is the molecular diffusion coefficient and τ_S is the average lifetime of a reactant molecule of species S (for an example on how to estimate τ_S, see Ref. 20). Within such a perspective, if one considers a cubic domain with side length L, and if

$\min _{S} {L^{S}_{\rm {Kuramoto}}} \gg L$

$minSL Kuramoto S≫L$⁠, then it is argued the particles of each species will diffuse fast enough to prevent local inhomogeneities.

Now, as a lower limit, Baras and Malek-Mansour also suggested the RMFP “to not compromise the separability of reaction and diffusion viewed as independent elementary processes.”¹⁸ In particular, the side length of the subvolume should be much larger than any reaction radius such that molecules can be considered uncorrelated (fully dissociated) within single subvolumes.⁵ 

A more rigorous lower limit for a discretization in three dimensional (3D) has been derived in Ref. 23, where it has been shown that the critical value of a subvolume's side length h is

where k is the macroscopic rate constant of the bimolecular reaction A + B → …, D_A (D_B) is the diffusion constant of A (B), and β_∞ ≈ 0.25272. It is important to note that h_crit stems from the derivation of modified bimolecular reaction rates for chemical kinetic systems in the limit of fast diffusion. In fact, it defines the lower limit when the modified bimolecular reaction rates are no longer applicable. As we will show below, such a lower limit also exists for the discretization length scale in 2D.

Lastly, an overview of discretization criteria (not including those stemming from correction methods) is shown in Table I.

The consideration of semi-permeable or reactive boundaries allows for particles to be absorbed by diffusion through the boundary with a user-specified probability. In such scenarios, choosing the rate for boundary jumps in RDME models is not straightforward.

In Ref. 24, the authors study a 1D stochastic position jump process with fixed time stepping and reactive boundary conditions, among other stochastic processes. By using an experimentally measurable reactivity constant ρ, they show how the probability for boundary absorption has to be chosen in order to match a corresponding partial differential equation (PDE) model with one Robin and one Neumann boundary condition (here, we use the symbol ρ instead of the K_r in Ref. 24, in order to avoid confusion with the symbol K that is defined below as a measure of the domain's discretization). In the stochastic process, the probability for a particle to jump across the boundary (R) as opposed to being reflected (1 − R) is R = hρ/D, where D is the particle's diffusion constant and h is the distance between jump positions. Thus, for given h, D, and ρ, one can calculate the value of the probability R. However, as the probability R is limited to values below or equal to 1, this imposes an upper limit on h, namely,

In this case, the value of h can be varied alongside R, while the values of ρ and D should be left untouched. Variations of the latter would necessarily imply changing the problem specification (i.e., the experimentally measured parameters).

On the other hand, a very large value of ρ (close to full absorption) in Eq. (1) might introduce an upper limit for h that is simply too low for the RDME approach to still be considered useful or computationally feasible. For instance, if we consider a diffusion coefficient D = 0.1 and boundary reactivity ρ = 10, it follows from Eq. (1) that h ⩽ 0.01. In other words, a minimum of 100^d subvolumes is required for a domain with length scale L = 1 in d dimensions. In fact, if one considers a fully absorbing boundary (ρ = ∞), an unfeasible value of h = 0 follows. Moreover, one cannot compensate for h > 0 by choosing an appropriate R. The latter shows how the problem of discretization-dependent solutions has not been fully addressed, despite the clear understanding of the relationship between PDE boundary conditions and probability of absorption. In this paper, we will show how one can keep the physical description of the reactive boundary with a fixed ρ and instead modify the jump rate for boundary diffusion by using a permeability and discretization-dependent correction factor, essentially removing the restriction on h.

In a regular mesh, the rate k_d of a diffusion event

$Ai→kdAj$ (i.e., a molecule “jumping” from subvolume i to subvolume j) is calculated^1,25 as k_d = D/h², where D is the diffusion constant. A more general rate formula has been recently proposed in Ref. 26 but assumes an underlying Cartesian mesh, while Engblom et al.²⁷ derived jump rates from finite element discretizations of the corresponding reaction-diffusion equation composed of triangles/tetrahedra in two/three dimensions. Nonetheless, both approaches recover D/h² for a regular mesh.

In what follows, we propose a methodology for deriving discretization-dependent factors correcting boundary jump rates.

We introduce a discretization- and permeability-dependent correction factor β for diffusion jumps through the boundary, that is, for molecules leaving the domain from a subvolume adjacent to the boundary. The new rate formula for such diffusion jumps is β × D/h². Using this rate in the associated RDME, we can derive equations for the stochastic mean (first moment) in each subvolume at steady state. However, in addition to the first moments, the correction factor β introduces one extra degree of freedom to the equation system. By imposing one more condition, namely, that the steady-state stochastic mean summed over all subvolumes equals a pre-defined value (σ), we are able to solve the equation system for β. In order to determine the total mean σ, we require a reference solution such as that of the corresponding PDE or from off-lattice simulations.

To illustrate our methodology, let us consider a test scenario where we assume a quadratic domain [0, 1] × [0, 1] with a uniform mesh with side length h = 1/K (K ⩾ 1), consisting of K² SV with area ν_K = h². Furthermore, let A_ij be the number of molecules of species A in SV (i, j) and D_A their uniform diffusion coefficient, while the domain is considered to be empty as an initial condition. For the sake of simplicity, let us first consider only one chemical reaction, namely, the constitutive creation of A_ij in each subvolume:

$∅→ksAij$⁠. The corresponding propensity is given as a_ij(t) = k_sν_K. Aside from this chemical reaction, we allow diffusion reactions between neighboring subvolumes and jumps through the boundary, which is assumed to be fully absorbing. The full derivation of the corresponding RDME can be found in Eq. (A1).

In a deterministic context, such a reaction-diffusion scenario can be described by the following PDE,

where u(x, y, t) is the concentration of molecules A at point (x, y), 0 < = x, y < = 1, at time t > 0, while u(x, y, 0) = 0. We can describe the fully absorbing boundaries of the domain with zero-Dirichlet boundary conditions u(0, y, t) = u(1, y, t) = u(x, 0, t) = u(x, 1, t) = 0. For semi-permeable boundaries (i.e., boundaries other than fully absorbing), Robin boundary conditions D_A∇u = −ρu are used instead. It should be noted that the latter is equivalent to the zero-Dirichlet boundary condition when ρ = ∞. The solution of the PDE using parameter values

$k¯s=100$ and D_A = 0.1 at steady state is shown in Fig. 1 and, quite logically, the steady-state solution of the PDE corresponds to the solution of Poisson's equation ∇²u = −k/D_A.

In a first step, we compare the (uncorrected) RDME for fully absorbing boundary jumps with the solution of the PDE. We set k_d = D/h² for boundary jumps as well, as one could argue the jump rate through any face of a subvolume should be equal everywhere. In Fig. 2, one can see a large discretization length scale h (small K) can produce considerable errors in the RDME solution, even for such a simple reaction-diffusion system on a uniformly discretized domain. As could be naturally expected, the RDME solution approaches that of the PDE solution for increasing values of K.

In order to provide a more accurate result, we can derive a set of K² + 1 equations with K² + 1 unknown variables, out of which 4K − 4 equations are nonlinear (as described in detail in Appendix  A). The unknown variables are the average number of molecules A at steady state in each SV

$v$ (M_v), and the correction factor β. One can obtain values for β and solutions for the mean of the RDME at steady state for each subvolume, by numerically solving this system of equations for different values of K and starting from an initial guess (e.g., the steady state of a reference solution). The dependency of β on K and ρ for this sample scenario is shown in Fig. 3. One can directly compare RDME and PDE solutions for each subvolume

$\mathbf {v}$

$v$ by subdomain integration (Fig. 4). As can be readily observed, the error map shows good agreement between the RDME with corrected boundary diffusion rate and the PDE domain integration results (Fig. 4(b)).

The position jump process described in Ref. 24 can be equivalently described by a RDME. This becomes evident by noting that the probability for jumping through a boundary face is R × D/h². Hence, the corresponding RDME equals Eq. (A1), where β is exchanged for the boundary jump probability R. In other words, the probability R functions as a correction factor and is directly comparable to β. Figure 5 shows the differences between R and β depending on the discretization value K for different values of ρ (D = .1 and k_s = 100). As was mentioned before, R is bounded by 1, hence imposing a lower limit for K given ρ and D. In contrast, the correction factor β yields accurate results using values of K beyond this limit. Nonetheless, one should note that the difference between β and R vanishes in the limit of K → ∞ or ρ → 0. The observation with ρ → 0 can be intuitively expected since the boundary condition approaches the no-flux condition, where no correction is needed.

Our proposed methodology readily generalizes to solving problems with any degree of adsorption in three dimensions, and scenarios where subdomains have distinct boundary conditions (Fig. 6).

Additionally, our simulations suggest that the correction factor is basically independent of any additional zero- and/or first-order reactions (results not shown). In other words, re-deriving the value of β while considering the appropriate forcing in the reference solution does not yield a significantly different match. In fact, earlier results support this finding by stating that, in the limits of h → 0 and δt → 0, additional first-order and/or higher order reactions will not influence the Robin boundary condition.²⁴ In such cases, it is much easier to obtain the value R, since β and R agree in the limit h → 0. However, away from this limit, one should obtain the value of β instead. In Sec. III F, issues related to the addition of higher order reactions, such as closure of moments and choosing an appropriate reference solution, will be discussed.

As already mentioned earlier, the reference solution is required to define an additional constraint and, hence, reduce the extra degree of freedom that the correction factor β introduced to the set of first moment equations. Here, the reference solution is used to determine the total stochastic mean and β is calculated such that the RDME yields the same total mean. The key question at this point is discerning which reference solution is most appropriate.

Whenever molecular concentrations are high (i.e., when the intermolecular distance is relatively short) and in the absence of noise-induced phenomena (such as signal amplification/damping, oscillations, bursts, or switches, to name a few), the deterministic and stochastic approaches for higher order reaction systems may yield sufficiently close results in the first moment. In the case of systems solely composed of zero- and/or first-order reactions, these approaches yield identical solutions. In such cases, the CME and stochastic differential equation (SDE) systems are identical in their first two moments and, since the drift term of the SDE system is linear, their first moment is also identical to the ordinary differential equations (ODE) solution.

However, in cases where the intermolecular distance is relatively large and stochastic effects are not negligible, solutions from particle simulators become essential. The consideration of particle simulator solutions as reference dynamics does not significantly affect the accuracy or performance of our approach. Nonetheless, and as could be naturally expected, good approximations of such reference solution require large numbers of sample simulations. In order to illustrate the applicability of our method to such scenarios, we performed a series of tests where a reference solution from the particle tracking software Smoldyn (Ref. 28) was obtained. For instance, we derived appropriate correction factors and obtained accurate RDME dynamics for various absorption probabilities using the above described system of molecule synthesis along with boundary diffusion absorbtion (data not shown).

In addition, we derived correction factors for a reaction system combining molecular synthesis of species A with a unary reaction

$A→kaB$⁠, and degradation of B solely through the boundary. For such a system, Fig. 7 presents mean trajectories and analytic results using 0% and 10% boundary absorption probability for molecules A and B, respectively (for the interpretation of absorption probabilities we refer to the Smoldyn manual). We derived a boundary correction factor β from the corresponding first moment equations and compared the results with the dynamics obtained when using the boundary factor

$(\hat{\beta })$

$(β̂)$ from a simpler system (synthesis and diffusion of one species with absorption probability of 10%; data not shown). Our results show for both approaches dynamics that are very similar to Smoldyn's reference solution, while analytical solution and corresponding ensemble average seem to match slightly better for the tailored correction factor (β).

Lastly, computing the coefficient of variation (CV) of both simulation types (off-lattice vs. corrected spatial SSA) shows a slight scaled difference between both approaches, but such difference can only originate from the second moment, as the first moment has already been shown to match. Thus, the lower CV could naturally stem from using strictly different simulation approaches. However, it is worth mentioning that the CV can many times yield misleading results, especially so in the case of nonlinear dynamics.²⁹ In our case, such concern goes a step further. Namely, the CVs of the uncorrected and corrected spatial SSAs were observed to be similar for species A, while for species B the uncorrected spatial SSA would apparently be a slightly better match (results not shown). Nonetheless, it is worth remembering that the first moments of B in the uncorrected and corrected SSAs differ considerably and, hence, any better match obtained from the uncorrected version naturally implies large differences in the second moment. Thus, the CV in this case would misleadingly prefer an erroneous simulation regime where none of the first two moments for B actually match.

Binary reactions require additional discretization-dependent corrections in the RDME. Recently, such corrected binding/unbinding rates have been derived for different scenarios and with different aims. In Ref. 23, a bimolecular reaction rate has been inferred for a two-reaction system in 3D such that the stationary distribution of the corresponding CME matches the RDME for any uniform discretization. This is valid given that the system can be considered homogeneous in the case of fast diffusion and, by consequence, RDME results should match the CME solution at any discretization length scale. Even though, this can also be done for bimolecular reactions in 2D (see Appendix  B for a derivation) the approach is not valid for any scenario where the homogeneity assumption does not hold. In Ref. 17, the authors present new mesoscopic rate formulas for diffusion-limited scenarios in 2D and 3D. Exact expressions for association and dissociation rates have been derived by considering a spherical domain discretized as shells, based on the microscopic Smoluchowski framework. For computational reasons these expressions were simplified without losing significant accuracy, while the transition to regularly meshed domains involves an extension of the conventional RDME scheme, where bimolecular reactions can now occur between molecules in neighboring subvolumes as well. We argue that appropriate bimolecular rates could also be derived by adding a correction factor in the corresponding set of algebraic moment equations. However, this is not as straightforward as in previous sections.

For kinetic systems with zero- or first-order reactions only, closed expressions for the moments can be derived. However, this is no longer the case when one or more bimolecular reactions are involved. In this case, the first moment equations become dependent on cross-order and higher order moments that in turn depend on other cross-order and higher order moments and so on. In addition, moment equations for RDMEs of a kinetic system with a binary reaction A + B → C depend on cross-moments of two (or even all three) participating species in different subvolumes. In fact, assuming particles can freely diffuse between any two subvolumes, cross-moments of all possible combinations of two (or three) subvolumes appear in the formulas (Appendix  C gives an example for such moment equations). Evidently, the number of cross-moments with species in different subvolumes increases rapidly with the number of subvolumes. Hence, for computing reasons, one needs to employ moment closure techniques to approximate the moment equations,^30,31 not only with a closure on the overall moment degree (the sum of all individual moments) but also with respect to the distance between subvolumes of participating species (shortest path between subvolumes). Naturally, any type of closure introduces an error that, as a rule of thumb, decreases as the truncation degree is increased. Preliminary studies of the reaction system in Appendix  C showed that the approximated moments deviate from the ensemble averages of corresponding spatial SSA simulations as K is increased, when truncating the system at summed order three and limiting cross-moments to species in identical subvolumes. Interestingly, considering moments of species in adjacent subvolumes increases the error as well, as opposed to improving the approximation. A similar, counterintuitive observation has been described,³¹ where a truncation at an even order k leads to a larger error then a truncation at k − 1. A more detailed analysis, in particular with respect to such spatial truncations, is left for future work.

In this paper, we have presented a method for obtaining accurate jump rates for systems with (semi-)permeable or reactive boundaries when using RDME-based simulation algorithms. For such purpose, it is sufficient to choose a minimal reaction system consisting of diffusion reactions (including boundary jumps) alongside zero- and first-order reactions, since such additional reactions do not alter the correction factor.²⁴ Such (un-corrected) system has to exhibit a non-zero steady state that varies according to changes in the domain discretization. Additionally, we require a reference solution, e.g., from a PDE or a particle simulator, the selection of which is entirely problem dependent. With this, the corresponding algebraic first moment equations are closed and only weakly nonlinear. There is only one nonlinear term, namely, the product of the correction factor with the first moment for subvolumes facing the domain boundary, and we can expect a unique solution for the correction factor. Once the correction factors are derived, they can be used in more complex systems, provided certain parameters such as diffusion rates or boundary permeability do not change.

The alternative approach in Ref. 24 is handy as it provides a closed formula for a correction factor. However, as we illustrated for simple reaction-diffusion scenarios, this is only useful as the discretization lengths h → 0. Moreover, in such case, the maximal h decreases alongside values of the diffusion coefficient D. Thus, in scenarios of very slow diffusion, h could become too small for the RDME approach to be actually useful. Our approach presented here remedies this issue.

Our method works in a straightforward manner whenever second-order reactions are not considered, due to the fact that such reactions require additional discretization-dependent correction factors. We argue here that our method can also be used in such scenarios, to derive the additional factors, where moment closure techniques have to be employed beforehand. However, the algebraic moment equations are likely to have more than one solution, due to the system being multistationary or, simply, because of the approximation itself. In fact, the number of solutions can depend on the truncation order, while most of those solutions are actually negative or complex.³¹ Obviously, when calculating the correction factor for a binary reaction together with those correction factors associated with permeable boundary jumps, additional equations comparable to Eq. (A9) are required. Those can usually be obtained from the reference solution (e.g., summed first moments of participating species).

Lastly, it should be noted an alternative method to approximate the first moments of the CME, based on Van Kampen's power series expansion of the master equation, has been proposed in Ref. 13 and is also applicable to RDMEs. However, it should be noted the so-called linear noise approximation has many drawbacks, even for solely temporal problems, as discussed in Ref. 29.

The authors would like to thank the reviewers for their helpful comments. The authors declare no funding sources for research and preparation of this manuscript. The authors took part of the modelling, drafted, read, and approved the final manuscript. A.L. performed simulations and derived moment equations.

We assume a quadratic domain [0, 1] × [0, 1] regularly meshed into K² SV with K ⩾ 1. Each squared subvolume has side length h = 1/K and area ν_K = (1/K)² = h². Subvolumes are indexed according to their position (i, j), 1 ⩽ i, j ⩽ K, in the domain. We denote with I the set of all subvolumes. More specifically, we denote with

• I₀ all inner subvolumes that do not have any boundary face,

• I₁ all subvolumes with exactly one boundary face, and

• I₂ the four corner subvolumes (characterized by two boundary faces).

Let

$E={[1,0],[0,1],[−1,0],[[0,−1]}$ be the set of 2D diffusion directions and

$\mathbf {E}_{\mathbf {v}} = \mathbf {E}$

$Ev=E$ the set of directions from SV

$\mathbf {v} \in I$

$v∈I$⁠. We divide this set into two subsets, namely, the set of boundary (outgoing) diffusion directions

$\mathbf {E}_{\mathbf {v}}^{\rm {out}}$

$Ev out$⁠, representing diffusion through the boundary, and the remaining set

$\mathbf {E}_{\mathbf {v}}^{\rm {in}}$

$Ev in$⁠. Diffusion from one SV

$\mathbf {v}$

$v$ to a neighbor SV

$\mathbf {v}+\mathbf {e}$

$v+e$ for

$\mathbf {e} \in \mathbf {E}_{\mathbf {v}}$

$e∈Ev$ is introduced in form of reactions

for all

$v∈I$ and

$\mathbf {e} \in \mathbf {E}_{\mathbf {v}}$

$e∈Ev$⁠.

Let A_v be the number of molecules of species A in SV

$v$⁠. Their diffusion constant is D_A. As already mentioned in the main text, we assume that A_v is constitutively created in each subvolume:

$\emptyset \stackrel{k_s}{\rightarrow } A_{\mathbf {v}}.$

$∅→ksAv.$ The corresponding propensity is given as a_v(t) = k_sν_K. We denote with

$p(\mathbf {n},t)$

$p(n,t)$ the probability that A_v = n_v for all SVs

$\mathbf {v} \in I$

$v∈I$⁠, where

$\mathbf {n}=\lbrace n_{\mathbf {v}}&#x007c;\mathbf {v} \in I\rbrace$

$n={nv|v∈I}$⁠. Let

$J^{\mathbf {e}}_{\mathbf {v}}(\mathbf {n})$

$Jve(n)$ be the jump operator for non-boundary jumps that effectively changes the entries of

$\mathbf {n}$

$n$ for SV

$\mathbf {v}$

$v$ from n_v to n_v + 1 and for SV

$\mathbf {v}^{\prime }=\mathbf {v}+\mathbf {e}$

$v′=v+e$ from

$n_{\mathbf {v}^{\prime }}$

$nv′$ to

$n_{\mathbf {v}^{\prime }}-1$

$nv′−1$⁠. Using this notation, we can define the operator for boundary jumps

$\tilde{J}^{\mathbf {e}}_{\mathbf {v}}(\mathbf {n}) = \mathbf {n} + \mathbf {\delta }_{\mathbf {v}}$

$J̃ve(n)=n+δv$⁠, where

$\mathbf {\delta }_{\mathbf {v}}$

$δv$ is set to be 1 for SV

$\mathbf {v}$

$v$⁠, and 0 elsewhere. For the moment, we assume that diffusive jumps through the boundary have a rate k_d = βD_A/h² and all others have a rate k_d = D_A/h².

Then the corresponding RDME takes the form

Let M_v be the average number (or first moment) of A molecules in SV

$v$⁠. Then, M_v corresponds to the sum of probabilities

$p(\mathbf {n})$

$p(n)$ weighted by n_v over all possible

$\mathbf {n}$

$n$⁠,

Multiplying Eq. (A1) by n_v and summing over all possible

$n$ yields K² equations for M_v,

where e₁, e₂, e₃, e₄ represent the four possible jump directions for a v ∈ I₀,

$e1′$⁠,

$e^{\prime }_2$

$e2′$⁠,

$e^{\prime }_3$

$e3′$ the three jump directions associated with a v ∈ I₁, and

$e^{\prime \prime }_1$

$e1′′$⁠,

$e^{\prime \prime }_2$

$e2′′$ the two corresponding jump directions for a v ∈ I₂ (4 equations).

Assuming steady state, Eqs. (A3)–(A5) simplify to the following set of equations:

Together with

where σ is the mean number of As in the entire domain at steady state (as obtained from the reference solution, PDE or particle simulations), we have K² + 1 equations (out of which 4K − 4 are nonlinear) for K² unknown M_v and β.

Bimolecular reactions can lead to discretization dependent inconsistencies in the solution of the corresponding RDME. In Ref. 24, a study of a two-reaction system consisting of a catalytic bimolecular degradation and a constitutive synthesis in a three-dimensional domain is presented. The authors point out that with increasing domain discretization (i.e., decreasing subvolume size) the probability distribution for the total number of molecules at steady state shifts to larger molecular numbers. As observed for a two-reaction system in 3D, the distributions of numbers of molecules increase alongside K in a 2D domain as well. However, such a shift is not as pronounced in 2D as it is in 3D.

Following rigorously the steps in Ref. 24, we can derive a modified kinetic rate also for bimolecular reactions on surfaces. Here, we report only the basic idea and main results for the 2D heterodimer reaction example (catalytic degradation of species A by species B, constitutive synthesis of A), considering a geometry uniformly discretized into K² compartments (with volume h², each). The interested reader may refer to Ref. 24 for further details.

We start with the system's RDME (for arbitrary K), where the standard rate for the bimolecular reaction k/h² is replaced by a variable λ. From there, we derive a formula for λ such that the stationary distributions of molecules A in the entire domain are identical for all K ⩾ 1 (hence, equal to the first moment of the corresponding CME at steady state).

For λ we obtain

where

In Ref. 24, the authors derive a critical discretization length scale h_crit from the correction formula for bimolecular reactions in 3D. In our case, we can derive an equivalent – though less strict – result in 2D, from Eq. (B1). The two factors in the denominator lead to two limits, the trivial h > 0 and β < (D_A + D_B)/k₁. Since β depends on K = 1/h, this implies the existence of a critical discretization level K, which can be numerically calculated for given parameters D = D_A + D_B and k₁. Figure 8 shows the stationary distributions for K = 20, both with and without corrected reaction rate, and the analytic solution of the stationary distribution obtained from the CME. The distributions were obtained from long time simulations with 10 000 sample points, assuming a zero-flux boundary condition.

For a chemical reaction system with N different species, let

$p(x,t)$ be the probability of being in state

$\mathbf {x}=x_1,\ldots ,x_N$

$x=x1,...,xN$ at time t, given an initial condition

$\mathbf {x}(0)$

$x(0)$⁠. Then, the associated multivariate moment is defined as

where

$i=(i1,...,iN)$ and

$\mathbf {x^i}=x_1^{i_1} \times \ldots \times x_N^{i_N}$

$xi=x1i1×...×xNiN$ (see Ref. 30). In case of a reaction-diffusion master equation, where identical chemical species are distinguished by the subvolume they reside, we use the more compact notation

$\mu ^{\mathbf {\chi }}_{\mathbf {j}}$

$μjχ$ where

$\chi \subseteq \lbrace X_{i_v}&#x007c; i\in \lbrace 1\ldots N\rbrace ,v\in \lbrace 1 \ldots M\rbrace \rbrace$

$χ⊆{Xiv|i∈{1...N},v∈{1...M}}$ (ordered first by species and then by subvolume) with M being the total number of subvolumes, and

$\mathbf {j}$

$j$ is the vector of corresponding moment orders. For instance,

$\mu ^{A_{v1},B_{v2}}_{2,1}$

$μ2,1Av1,Bv2$ is a cross-moment, involving the second moment of species A in subvolume v1 and the first moment of B in subvolume v2.

Assume the following system of reactions¹⁷:

$A→kdA∅$⁠,

$B \stackrel{k_{dB}}{\rightarrow } \emptyset$

$B→kdB∅$⁠,

$\emptyset \stackrel{k_s}{\rightarrow } C$

$∅→ksC$⁠,

$A + B \stackrel{k_{a}}{\rightarrow } C$

$A+B→kaC$⁠, and

$C \stackrel{k_d}{\rightarrow } A + B$

$C→kdA+B$⁠. Using the general formula for multivariate moments derived in Ref. 30, we obtain the first moment equations

with discretization length scale h and diffusion constants D_A, D_B, and D_C.

For the cross-moment

$μ1,1Av,Bv$⁠, one obtains

This equation depends on the first-order cross-moments

$μ1,1Av,Cv$ and

$\mu ^{B_v,C_v}_{1,1}$

$μ1,1Bv,Cv$⁠, cross-moments including second moments such as

$\mu ^{A_v,B_v}_{2,1}$

$μ2,1Av,Bv$ and

$\mu ^{A_v,B_v}_{1,2}$

$μ1,2Av,Bv$ as well as cross-moments of species in neighboring subvolumes, namely,

$\mu ^{A_{v},B_{v+e}}_{1,1}$

$μ1,1Av,Bv+e$ and

$\mu ^{A_{v+e},B_v}_{1,1}$

$μ1,1Av+e,Bv$⁠. The equations of the latter moments depend in turn on other moments of spatially “separated” species such as

$\mu ^{A_{v},C_{v+e+\hat{e}}}_{1,1}$

$μ1,1Av,Cv+e+ê$ and

$\mu ^{B_{v+\tilde{e}},C_{v+e+\hat{e}}}_{1,1}$

$μ1,1Bv+ẽ,Cv+e+ê$⁠, where

$\hat{e}\in \mathbf {E}_{\mathbf {v+e}}$

$ê∈Ev+e$ and

$\tilde{e}\in \mathbf {E}_{\mathbf {v}}$

$ẽ∈Ev$⁠, and leads to cross-moments between species in any two subvolumes.