# Reaction rates for mesoscopic reaction-diffusion kinetics

###### Abstract

The mesoscopic reaction-diffusion master equation (RDME) is a popular modeling framework, frequently applied to stochastic reaction-diffusion kinetics in systems biology. The RDME is derived from assumptions about the underlying physical properties of the system, and it may produce unphysical results for models where those assumptions fail. In that case, other more comprehensive models are better suited, such as hard-sphere Brownian dynamics (BD). Although the RDME is a model in its own right, and not inferred from any specific microscale model, it proves useful to attempt to approximate a microscale model by a specific choice of mesoscopic reaction rates. In this paper we derive mesoscopic scale-dependent reaction rates by matching certain statistics of the RDME solution to statistics of the solution of a widely used microscopic BD model: the Smoluchowski model with a Robin boundary condition at the reaction radius of two molecules. We also establish fundamental limits on the range of mesh resolutions for which this approach yields accurate results, and show both theoretically and in numerical examples that as we approach the lower fundamental limit, the mesoscopic dynamics approach the microscopic dynamics. We show that for mesh sizes below the fundamental lower limit, results are less accurate. Thus, the lower limit determines the mesh size for which we obtain the most accurate results.

## I Introduction

The reaction-diffusion master equation (RDME) is a commonly used mesoscopic model in the field of computational systems biology and it is a natural extension of the classical well-mixed Markov-process formalism for reaction kinetics gillespie (). Having a long history in the study of fluctuations in chemical reaction systems gardiner1976 (); Baras (), recently it has been successfully applied to study diverse biological phenomena such as yeast polarization Lawson:2013 (), pattern formation in E. Coli Howard (); FaEl (), and noisy oscillations of Hes1 in embryonic stem cells Sturrock1:2013 (); Sturrock2:2013 ().

In the RDME framework, spatial heterogeneity is modeled by dividing space into voxels in a computational mesh, where molecules are assumed to be well mixed inside each voxel. In individual voxels, reactions are simulated using the Stochastic Simulation Algorithm (SSA) gillespie (), while diffusion is accounted for through discrete jumps of molecules between voxels. Discrete diffusion and well-mixed SSA are combined in the Next Subvolume Method (NSM) ElEh04 (), an efficient kinetic Monte Carlo method. For moderate mesh resolutions, simulations of the RDME are typically orders of magnitude faster than microscopic particle-tracking models, such as the popular Green’s Function Reaction Dynamics (GFRD) algorithm ZoWo5a (); ZoWo5b (); SHeLo11 (). This contributes to the popularity of the method for applications where the system of interest needs to be studied on the minute to hour timescales that are typical for cellular events like gene expression, signaling and cell division. The RDME underlies software for spatial stochastic simulation such as MesoRD mesord (), URDME URDME_BMC (), pyURDME (www.pyurdme.org), and STEPS steps ().

Despite the proven usefulness of the RDME—and the extensive work put into speeding up simulations using approximate RosBayKou (); marquez-lago:104101 (); LaGiPe () and hybrid FeHeLo2009 () methods, as well as extending it to include additional transport phenomena HeLo:2010 (); Bayati:2013 () and to simulate it on complex geometries IsP (); EnFeHeLo ()—its fundamental numerical properties and its ability to approximate microscopic particle tracking models at high mesh resolution remains poorly understood. In order to discuss the accuracy of the RDME on small length- and timescales, we need to specify a fine-scaled alternative as our gold standard. A microscale model often utilized for that purpose is the Smoluchowski model, in which particles are modeled by hard spheres that diffuse according to Brownian motion, and reactions are modeled by a partially absorbing boundary condition at the surface of the spheres. This model has a long history in chemical physics, going back to ideas by Smoluchowski Smol (). In systems biology, the Smoluchowski model is being popularized through software packages such as E-Cell Tomita01011999 () and Smoldyn AnAdBrAr10 (). This paper is concerned with the accuracy of the RDME when viewed as an approximation to the Smoluchowski model.

The principal way in which the mesoscale and the microscale are connected is through the mesoscopic bimolecular reaction constant. A classical result by Collins and Kimball Kimball () provides effective rates in terms of the microscopic, intrinsic, reaction parameters of the Smoluchowski model. When scaled by the volume of the voxels, they can be used as approximate mesoscale reaction rates in the RDME for simulations in 3D. The same constant was derived more recently from first principle physics by Gillespie gillespie:164109 (). The constant is valid when voxels are large in comparison to the molecules. Since, in the conventional implementation, reactions occur only between molecules occupying the same voxel, the average time until molecules react diverges with vanishing size of the voxels Isaacson2 (); ErCha09 (); HHP (). A consequence is a lower bound on the size of the voxels, below which no mesoscopic reaction rate can be chosen so that the average reaction times match in the RDME and Smoluchowski models HHP ().

In previous work HHP () we analyzed the scenario of a single, irreversible, bimolecular reaction on a Cartesian mesh. For the case of perfect absorption, we obtained analytical lower limits on the mesh size in both 2D and 3D. Above those limits, it is theoretically possible to construct a mesoscopic rate such that the mean binding time in the two models match. Below that critical mesh size, no such rate can be constructed. Hence, in the presence of bimolecular reactions, a fundamental limitation of the RDME results from the inherent bound on the accuracy to which we can represent diffusion. As a direct consequence of our previous analysis, we also obtained mesoscopic reaction rates which ensure that the mean binding time in the RDME matches that of the Smoluchowski model for an irreversible bimolecular reaction.

In this paper we extend our previous analysis. First, we study the case of reversible reactions, and ask whether there are additional constraints on the admissible mesh sizes compared to the irreversible case. We derive a critical voxel size under the following three assumptions: the average time until a reaction fires should match between mesoscopic and microscopic models, the steady state levels should match, and the mesoscopic dissociation rate should be smaller than or equal to the intrinsic dissociation rate. The last condition is necessary for the dissociation to make physical sense, and if we are to match not only equilibrium distributions but also transient solutions. The result establishes the previously obtained critical voxel sizes—in the perfectly absorbing, irreversible, case—as a fundamental lower limit for the more general reversible case. Importantly, this means that there will be a non-trivial lower limit for the mesh size, independent of the intrinsic microscopic reaction parameters. We also study the accuracy of our reaction rates as compared to the Smoluchowski model, and show good agreement, both at steady state and during the transient phase, between mesoscopic and microscopic simulations as the mesh size approaches the critical lower limit. In particular, we show how the multiscale propensities can provide a better mesoscopic approximation to a diffusion limited model of a MAPK cascade, compared to the widely used propensities of Collins and Kimball Kimball (). It has not previously been possible to accurately simulate this model with a fully local RDME, although it has been simulated successfully with a non-local extension of the RDME FBSE10 (), and using a hybrid microscale-mesoscale method HeHeLo ().

For simplicity we have chosen not to write out units explicitly. We are using SI units throughout.

## Ii Background

A system of molecules of different chemical species, diffusing and reacting inside a finite reaction volume, can be modeled at several different scales, and the accuracy of the different models depends on the properties of the system. In this section we briefly review the microscopic Smoluchowski model and the mesoscopic RDME model, and discuss previous work on connecting the two models.

### ii.1 Microscopic scale: the Smoluchowski model

Several microscopic models have been proposed and studied in some detail ErCha09 (); SamDoiSDLR (); AnB (). We have chosen to focus on the Smoluchowski model, given the extensive attention it has received for instance in ZoWo5a (); ZoWo5b (); SHeLo11 (). In the Smoluchowski framework, molecules are modeled as hard spheres. The radii of the spheres are referred to as the reaction radius of the molecules. The Smoluchowski equation, extended with a Robin boundary condition at the sum of the reaction radii, determines the probability of a reaction occurring between colliding molecules.

Let and be the positions of two molecules. Consider their relative position . The governing equation is given by

(1) |

where is the probability distribution function (PDF) of the relative position at time , given the relative position at time , and where is the sum of the diffusion constants of the molecules. The boundary condition is given by

(2) |

where

(3) |

Here is the sum of the reaction radii of the molecules, and is the microscopic, intrinsic, reaction rate. The initial condition is .

To update a pair of molecules we solve for , and then sample the new relative position at time . Single molecules are updated by sampling from a normal distribution in all directions. The time until a dissociation is sampled from an exponential distribution, and the resulting products are placed at a distance equal to . A system of molecules becomes an -body problem, and a direct solution is generally unattainable. Instead we can simulate the system with the Green’s function reaction dynamics (GFRD) algorithm ZoWo5a (); ZoWo5b (). The core of the algorithm is a reduction of the problem to a collection of single- and two-body systems, accomplished through an appropriate restriction of the time step of the algorithm. In HeHeLo () the GFRD algorithm was extended to include complex boundaries, and in SHeLo11 () improvements were suggested with the aim of making the algorithm more flexible and efficient.

The GFRD algorithm is efficient for dilute systems where the free space between molecules allows for an efficient grouping in pairs while using a relatively large time step, and the computational benefit over brute-force BD methods can be orders of magnitude. If the system contains some species that are present in higher copy numbers, or if a very high spatial resolution is not required, the mesoscopic RDME can in turn be orders of magnitude faster than GFRD.

### ii.2 Mesoscopic scale: The reaction-diffusion master equation

The RDME extends the classical well-mixed Markov-process model gillespie (); GillespieRig () to the spatial case by introducing a discretization of the domain into non-overlapping voxels Gardiner (). Molecules are point particles and the state of the system is the discrete number of molecules of each of the species in each of the voxels. A common choice for the discretization is a uniform Cartesian lattice, where each voxel is a square with area (2D), or a cube with volume (3D). Simulations can also be conducted on unstructured triangular and tetrahedral meshes for better geometric flexibility EnFeHeLo (). The RDME is the forward Kolmogorov equation, governing the time evolution of the probability density of the system.

For brevity of notation, we write for the probability that the system can be found in state at time , conditioned on the initial condition at time . For a general reaction-diffusion system, the RDME can be written as

where denotes the -th row and denotes the -th column of the state matrix where is the number of chemical species. The functions define the propensities of the chemical reactions, are stoichiometry vectors associated with the reactions. and have the same meaning as for a well mixed system but are now defined for each of the voxels. are propensities for the diffusion jump events, and are stoichiometry vectors for diffusion events. has only two non-zero entries, corresponding to the removal of one molecule of species in voxel and the addition of a molecule in voxel . The RDME is too high-dimensional to permit a direct solution. Instead realizations of the stochastic process are sampled, using algorithms similar to the SSA but optimized for reaction-diffusion systems ElEh04 ().

The propensity functions for the diffusion jumps, , are selected to provide a consistent and local discretization of the diffusion equation, or equivalently the Fokker-Planck equation for Brownian motion GillHellPetz (). For a uniform Cartesian grid, a finite difference discretization results in diffusion jumps with propensities

(5) |

where is the diffusion coefficient of and is non-zero only for adjacent voxels. For a triangular or tetrahedral unstructured mesh, finite element or finite volume discretizations result in propensities that account for the shape and size of each voxel EnFeHeLo (). With this model of diffusion, setting reactions aside, the solution of the RDME will, with vanishing voxel sizes, converge in probability to Brownian motion.

In the case of mass action kinetics, the propensity functions for bimolecular reactions take the form

(6) |

We will refer to as the mesoscopic association rate. In practical modeling work, is often obtained by fitting the model to some phenotypic experimental observation, or provided directly as a macroscopic or mesoscopic model parameter obtained from experiments. If the association rate is instead given in terms of the microscopic reaction rate, , a result of Collins and Kimball Kimball () provides the effective rate, , of the system, which can be used for mesoscale simulations in 3D through , where

(7) |

Gillespie derives the same relation by applying classical results from gas kinetics, and calls it the diffusional propensity function gillespie:164109 (). From Gillespie’s physically rigorous derivation it is possible to relate the intrinsic reaction rate to fundamental physical constants.

A natural question, given model parameters and a mesh size, is how well a mesoscopic simulation can capture the microscale dynamics. It has been shown that the solution of the RDME diverges with respect to the Smoluchowski model Isaacson (); HHP (). Due to the point particle assumption, bimolecular reactions occur with successively lower probability, to eventually vanish in the limit of small voxels.

This means that there is a lower bound on the voxel size, below which bimolecular reactions cannot be accurately simulated in the RDME. There is also an upper bound on the voxel size, above which the reaction-diffusion dynamics will be insufficiently resolved. The question of how to choose the voxel size for sufficient accuracy has not yet been satisfactorily answered, but some attempts at establishing lower bounds on the voxel size have been made. A trivial bound on the voxel size follows from physical arguments GPS:2014 (); the voxels must be large enough for the reacting molecules to remain dilute and well mixed inside the voxel. This condition translates to

(8) |

Apart from physical common sense, the above condition is an explicit assumption in the derivation of (7). Only if this assumption is valid can the rate constant (7) be expected to provide an accurate simulation with respect to the Smoluchowski model. Unfortunately, condition (8) offers little guidance on how to choose the mesh size in practice.

The reaction constant (7) depends on the microscopic parameters but not on the spatial discretization. By allowing the rate to depend on the mesh, and by matching certain properties of the microscopic model, it is possible to derive alternative forms for that perform better than (7) for diffusion limited systems and for fine meshes FBSE10 (); HHP (). Those approaches also yield constants in 2D, where no expression based on a physical derivation is available. In FBSE10 () the rates are derived by matching the mean equilibration time of a reversible reaction on a spherical discretization. These rates are then used on Cartesian lattices.

Let be the mean binding time for two molecules in the Smoluchowski model in dimension . In the case of a single irreversible association reaction on a uniform Cartesian discretization of a square or cube of side length , Hellander et al. HHP () showed that if the mesoscopic propensity is chosen as , where

(9) |

then the mean binding time on the mesoscopic scale will match the mean binding time in the Smoluchowski model. For , this is possible down to in 3D, and in 2D. Below , the matching of the mean binding time will not be possible in the perfectly absorbing case. For simple geometries, such as a disk or a sphere, can be obtained analytically. For other geometries, provided that , the analytical result for a disk or a sphere with matching volume provides an excellent approximation. The analytical lower bounds on the voxel sizes are obtained by considering the extreme case of and using the above mentioned approximation for . In ErCha09 (), based on another assumption not involving the microscopic Smoluchowski model, Erban and Chapman arrived at an expression in 3D that is similar to (9) and that establishes the same critical mesh size.

In what follows, we set out to expand on this theory with the aim of an improved understanding of the range of the mesh sizes for which the RDME will accurately approximate the Smoluchowski model. In particular, we derive critical mesh sizes for more realistic kinetics like reversible bimolecular reactions, and we show that under mild assumptions, the rates (9) are effectively independent of . We also obtain error estimates that provide a way to estimate the needed mesh resolution to achieve a certain accuracy in the rebinding time distributions.

## Iii Results

The case of a single irreversible reaction was studied in HHP (). The analysis provided analytical lower bounds on the mesh size in 2D and 3D only for the case of , and reaction rates for matching mean association times. The reaction rates depended on the size of the domain. As reactions occur locally in space, that dependence is not intuitive.

In this section we expand on that theory. First we show that the reaction rates are independent of the size of the domain, under the assumption that the domain is much larger than the molecules. From this follows analytical lower bounds on the mesh size for the case of an irreversible reaction with .

Second, we study a reversible reaction on a square or cubical domain and proceed to derive mesoscopic reaction rates under the following three assumptions: the average reaction time should match between mesoscopic and microscopic models, the steady state levels should match, and the intrinsic dissociation rate should be larger than or equal to the mesoscopic dissociation rate. Under these assumptions we derive a lower bound on the mesh size independent of the reaction rates.

We also show how the multiscale mesoscopic reaction rates behave in the limit of small voxels as well as in the limit of large voxels, effectively connecting the microscopic and mesoscopic scales. Finally, we provide error estimates that relate the mesh size to the error in rebinding-time distributions, and show that the mean mesoscopic rebinding time approaches the mean microscopic rebinding time, as the mesoscopic mesh size approaches the lower bound.

### iii.1 Irreversible reactions

The analytical expressions for the lower bounds on the voxel size derived in HHP () are valid for the case of irreversible reactions with perfect absorption. In this section we assume , to obtain analytical expressions for the lower bounds in the general case of . We show that the reaction rates are independent of the size of the domain, thus depending on local parameters only.

###### Theorem 1.

Let be the mesoscopic reaction rate in dimension , and assume that . Then

(10) |

where , and

(11) |

where

(12) |

Let

(13) |

Then is a well-defined reaction rate only for , and

(14) |

###### Proof.

In HHP (), the reaction rate was derived under the assumption that . This effectively implies that , since we have adopted the basic assumption (8). Let and be the mean association times for uniformly distributed particles in dimension . While individual voxels may be too small for the Collins and Kimball approximation to be valid, the system as a whole may be well-mixed and dilute. Thus, if we assume that we can approximate the global mean reaction time by

which, when inserted into (9), yields

With defined by (7), we obtain as in (10). Thus, for large enough domains, the effect of the outer boundary is small in 3D, and the reaction rates are defined by local parameters only.

The situation in 2D is more complicated. In FBSE10 () they derive for a disk. Assuming , this will provide an excellent estimate of in the case of a square. Thus

where

yields

### iii.2 Reversible reactions

In this section we extend the analysis to the reversible case. We limit our considerations to the conventional RDME, thus allowing reactions only between molecules occupying the same voxel. The system consists of one -molecule and one -molecule, diffusing on a Cartesian lattice consisting of voxels of width , with periodic boundary conditions. The molecules react reversibly through the reaction

(17) |

Henceforth we assume that , the microscopic association rate, and , the microscopic dissociation rate, are given model parameters, and that is defined by (10).

It remains to define . A plausible approach would be to match the steady-state levels of the species at the mesoscopic and microscopic scales. For a system of one - and one -molecule, the steady state is governed by the ratio of the average time the molecules are unbound to the total time. Thus, to match the steady state at the different scales, we choose such that

(18) |

where and are the respective times until an - and a -molecule rebind, given that they have just dissociated. The quantities and are the average times until a -molecule dissociates (thus and ).

We can compute analytically by expressing it in terms of the mean binding time . We note that the rebinding time is given by

(19) |

When the molecules have reached the same voxel, the system is in the same state as immediately following a dissociation. Thus, by subtracting from the total mean binding time the time it takes to reach the same voxel, we are left with the average rebinding time.

Unfortunately is not easily computed by analytical means (for general geometries), but by noting that similarly as in the mesoscopic case, we obtain the estimate

(21) |

for . The argument is the same as for the mesoscopic case; represents the time until two molecules are in contact for the first time. By subtracting that time from the total time for the two molecules to react given a uniform initial distribution, we are left with the rebinding time. By using (20) and (21) in (18) it follows that

(22) |

which, after some simplifications, results in

(23) |

When rearranged, (23) can be recognized as the detailed balance condition. Maintaining this relation between the rate constants is sufficient to ensure that the equilibrium values of and are the same in the two models.

However, though we can match the mean association time as well as the equilibrium values for a given voxel size, we may have a mesoscopic dissociation rate that is faster than the microscopic dissociation rate. That would result in more dissociation events on the microscopic scale, and since the microscopic scale is assumed to be more fine-grained, this is unphysical. This can be seen by considering that the inverse of the mesoscopic dissociation rate is a combination of the expected time for a complex to break apart and the time for the products to become well-mixed by diffusion inside a voxel. The mesoscopic dissociation rate thus includes the possibility of fast microscopic rebinding events that occur before the molecules have become well-mixed. Thus, while matching the mean association time and the equilibrium values of the molecules, we may end up with incorrect transient dynamics due to too fast dissociations on the mesoscopic scale. The intrinsic microscopic dissociation rate is therefore an upper bound on the mesoscopic rate. As a consequence we want to determine for what size of we can match the mean binding times while satisfying

(24) |

###### Theorem 2.

Let

Then is the smallest for which we can satisfy as well as , and we have

(26) |

###### Proof.

We obtain by using (23) in (24):

(27) |

which holds if and only if

(28) |

However, from (10) we have

Thus we find that (28) is equivalent to

(29) |

since is excluded by the condition as shown in (16). Since and (we do not consider the case here), (25) follows.

We now observe that , and that, immediately from (10) and (13), it follows that . The last equality is an immediate consequence of (25), and we obtain (26).

∎

In Figure 1 we show how relates to .

### iii.3 Some limit cases

In this section we investigate the behavior of the reaction rate in some limit cases. We would expect to behave similarly to for large voxels. In the limit of small voxels, it is of interest to see how the mesoscopic reaction rates relate to the microscopic reaction rates.

###### Corollary 1.

For we have

(30) |

Also, as , with constant and , we have

(31) |

In the limit of we obtain

(32) |

In 3D, (32) implies that

(33) |

as .

###### Proof.

In words, as we approach the critical mesh size in the mesoscopic model, the mesoscopic rates approach the microscopic, intrinsic rates. For large voxels in 3D, we note that , as expected, approaches .

As an immediate consequence we have

###### Corollary 2.

For and , the following holds:

(34) |

###### Proof.

For , we match the mean association time, dissociation time, and the mean rebinding time. For any below , we cannot simultaneously match both the association and dissociation times. Consequently the rebinding dynamics will be less accurately captured on the mesoscopic scale for than for .

In Figure 2 we illustrate these limits in 3D for different values of the intrinsic reaction rate . The more diffusion limited the reaction is, the larger the discrepancy between and becomes.

### iii.4 Error estimates

Given a system with known intrinsic reaction rates, we would ideally want to know how to choose the mesh size for sufficiently accurate mesoscopic simulations. While solving this is hard for a general system, we can choose the mesh size such that we limit the relative error in the average rebinding times. We do not guarantee an accurate mesoscopic solution by doing so, but if the error in the average rebinding time is large, and we have fine-grained dynamics that we wish to capture, we may have to reduce the voxel size to decrease the error.

Now consider a reversible reaction, with intrinsic reaction rates given by and .

###### Corollary 3.

For a given error tolerance , the following holds:

(35) |

if and only if

(36) |

where

(37) |

Furthermore, for , we have

(38) |

In 2D we have

(39) |

In 3D, we obtain

(40) |

as

(41) |

For all , it holds that

(42) |

in 3D.

###### Proof.

We see that as the reactions become very diffusion limited, the difference between the mesoscopic and microscopic reaction rates can grow large, since as . The less diffusion limited a reaction is, the closer the reaction rates will be (and the difference is bounded). This makes intuitive sense, since less diffusion-limited reactions mean that the system is more well-mixed; thus the system can be accurately simulated on a coarser mesh.

In Section IV.2 we demonstrate how this theory can be applied to increase the understanding of the behavior of a relevant biological system.

## Iv Numerical experiments

In this section we present two numerical examples that will demonstrate the scope of validity of the mesoscopic reaction rates derived above. In the first example we consider the rebinding time of a pair of molecules. We compute the distributions and compare mesoscopic results for varying to microscopic simulations.

In the second example we study a model proposed by Takahashi et al. in TaTNWo10 (). It was shown to have fine-grained dynamics, captured at the microscale but not at the deterministic level; we will simulate it at the mesoscopic scale, and show that as the mesh size approaches , the mesoscopic and microscopic scales agree.

### iv.1 Rebinding-time distribution

Consider one molecule of species and one molecule of species , subject to a reversible reaction

(49) |

in a cubic domain of width .

We showed in Section III.3 that for the mean rebinding time at the mesoscopic scale will agree with the mean rebinding time at the microscopic scale. This does not, however, automatically guarantee that this particular choice of mesh size will yield the best agreement between the distributions of the rebinding times at the different scales.

In Figure 3 we plot the distribution of the rebinding times for different mesh sizes, and compare to the distribution obtained by simulations at the microscopic scale. For we get a distribution that matches the microscopic results well for . For other values of , we get a distribution shifted relative to the microscopic distribution. For , the microscopic simulations behave differently than the mesoscopic ones, regardless of mesh size. This is due to the fact that at this timescale the spatial resolution is coarser than the temporal resolution, thus being the limiting factor for the accuracy. During a time , a molecule diffuses an average distance proportional to in each direction. Thus, for , the distance the molecule diffuses is less than the size of a voxel.

On the mesoscopic scale, the distribution of the rebinding time is approximately exponential at short time scales (time scales smaller than the time it takes for a molecule to diffuse on average the distance of a voxel.) That gives the first plateau in Figure 3. At longer time scales the diffusion causes the rebinding time not to be exponential. In this region the mesoscopic simulations behave similar to the microscopic. At longer time scales, the distribution is again approximately exponential; once the molecules have become approximately well-mixed in the enclosing volume the assumption of an exponential reaction rate is quite accurate. This is seen as the second plateau in Figure 3.

In Figure 4 we plot the mean rebinding time as a function of , and show that for the mean rebinding times match between the different scales.

### iv.2 MAPK cascade

An example of when small errors in the transient dynamics of bimolecular equilibration processes can have a large impact on the system’s dynamics was given in TaTNWo10 () to illustrate the need of the microscale resolution provided by the GFRD algorithm. The model considered is two steps of the omnipresent mitogen activated phospatase kinase (MAPK) cascade. Here, a transcription factor is phosphorylated in two steps by a kinase and dephosphorylated by a phosphatase :

(50) | |||

(51) | |||

(52) | |||

(53) | |||

(54) | |||

(55) | |||

(56) | |||

(57) |

During the phosphorylation and dephosphorylation steps, the kinase and phosphatase turns into an inactivated form and . This can model e.g. a conformation change due to conversion of ATP to ADP, resulting in the need to reactivate the enzymes before proceeding with the next reaction. If the timescale for this reactivation step is short and the system very diffusion limited, rapid rebinding of to the newly phosphorylated molecules can have a big impact on the overall system dynamics, as illustrated and discussed from a biological perspective in TaTNWo10 (). Numerically, this means that the system is very challenging to simulate with lattice based methods, due to the need for very fine spatial resolution in order to resolve the rebindings on the fast timescale. This was noted by Fange et al. in FBSE10 (), where length scale dependent rates were derived based on the ansatz that the equilibration time should match on the two scales for a spherical discretization. They managed to resolve the microscale dynamics of the model by using these propensities in combination with extending the RDME to allow for reaction events between molecules occupying neighboring voxels as well as molecules occupying the same voxel. Without that extension they were not able to resolve the microscale dynamics; below we show that with reaction rates as defined in (10) and (23) we are able to resolve the microscale dynamics without considering a non-local extension of the RDME.

In Figure 5 we show the results obtained when simulating this model using our local rates for different mesh resolutions. As can be seen, when is close to , we obtain a good approximation of the results of the GFRD algorithm from TaTNWo10 (). In the figure, we show the time until half-activation (i.e. the time to reach half the steady state level of ) for varying diffusion constants, making the system range from reaction limited to diffusion limited. As can be seen, for large values of , as expected, the rates proposed here and the rates of Collins and Kimball give similar results, but for the strongly diffusion limited cases, our rates result in a much better agreement with the microscale model. Notably, for small we obtain comparable accuracy to using with for , resulting in simulations that run approximately sixteen times faster, due to the scaling of the computational time.

In Figure 6a we have computed the upper bound on obtained by requiring that the relative error in the average rebinding time is bounded by . As expected, we can see that for the less diffusion limited cases, when is larger, we can choose the voxel size larger and still obtain accurate results. As becomes even larger, the well-mixed assumption will be satisfied, and the system can be simulated at a much coarser scale. For smaller values of the restriction on is quite severe, and we are required to approach in order to accurately simulate the system. This has an implication for the computational complexity of the simulations; a small makes the system less stiff, meaning it becomes less computationally challenging, but also forces us to choose a smaller voxel size, while a large makes the system more stiff, but at the same time allows for a larger voxel size. In Figure 6b we have computed the maximum error in the average rebinding time, as a function of but independent of the size of the voxels.

## V Discussion

As we have seen, by taking a multiscale approach and deriving reaction constants by matching certain statistics of the Smoluchowski model, more accurate simulations can be obtained compared to the classical approach using rates from Collins and Kimball. An important reason why this is possible is that we start out with a given discretization of space, and then derive scale-dependent multiscale propensities (discretization-first), while the CK-rate was not derived with a mesh in mind, and hence the propensity depends only the volume of the voxels and not their shapes. As a consequence, we can better approximate binding times for bimolecular reactions that need high spatiotemporal resolution.

While it is important to be able to accurately resolve the reaction kinetics of a given diffusion limited system, it may also be important to accurately resolve complex geometries and external and internal boundaries, modeling for example cell membranes. In this case, uniform Cartesian grids have distinct drawbacks compared to unstructured triangular and tetrahedral discretizations in that they require more voxels in order to resolve the boundaries EnFeHeLo (); URDME_BMC () resulting in unnecessarily long computational time. While we expect fundamental limits for a triangular and tetrahedral discretization to be close to the ones obtained herein for Cartesian grids, whether it will be possible to extend the approaches taken here to obtain sharp estimates and reaction rates also for unstructured meshes has yet to be seen. Another approach to the problem of making an RDME-type model approximate a microscopic model is taken in Isaacson (), where a mesoscale model is constructed by discretizing the Doi model SamConvergent (). This approach seems more directly amenable to be used on general grids, but it has not yet been applied to the Smoluchowski model, and except for the case of irreversible, perfect absorption SamDoiSDLR (), the relationship between the Doi and Smoluchowski models is not well understood.

We have illustrated in numerical examples that by using the propensities proposed here it is possible to accurately simulate the MAPK system discussed in TaTNWo10 () on the mesoscale using a purely local RDME implementation, and we have shown theoretically why is the optimal mesh size. The same system has previously been successfully simulated with an RDME-type model by Fange et al. FBSE10 (), using another set of multiscale propensities and by relying on a non-local implementation of the RDME in which molecules occupying adjacent voxels are allowed to react. Relying on neighbor interactions leads to an increased computational cost due to the increased number of updates in each step of the kinetic Monte Carlo algorithm and hence we should expect the propensities derived here to provide a computational speed advantage for comparable accuracy.

Although more efficient than a non-local implementation, the simulations herein require a uniformly fine mesh and hence expensive simulations for systems that require very high spatial resolution. Unless there are species in the model present in high copy numbers which would cause microscopic simulation to become very time consuming, it is not unlikely that an efficient implementation of e.g. GFRD is more efficient than the purely mesoscopic simulation. With this in mind, for very diffusion limited systems with multiscale properties, a compelling approach is the use of hybrid methods. Such a method, blending the RDME and GFRD algorithms, has previously been proposed by the present authors HeHeLo (), in which it was demonstrated that an accurate hybrid simulation of the MAPK model TaTNWo10 () can yield accurate results with only a small part of the system simulated on the microscopic scale. However, there are outstanding challenges in making such hybrid methods easy to use by practitioners. Dynamic and adaptive partitioning of the system into microscopic and mesoscopic parts according to accuracy requirements is needed both for computational efficiency and for robustness of the simulator. In this paper we advance the fundamental theoretical understanding of the RDME when viewed as an approximation to the Smoluchowski model, something that is a prerequisite to develop adaptivity criteria for hybrid methods.

## Vi Acknowledgment

The authors thank Per Lötstedt for constructive comments on the manuscript. This work was funded by NSF award DMS-1001012, NIGMS of the NIH under award R01-GM113241, Institute of Collaborative Biotechnologies award W911NF-09-0001 from the U.S. Army Research Office, NIBIB of the NIH under award R01-EB014877, and U.S. DOE award DE- SC0008975. The content of this paper is solely the responsibility of the authors and does not necessarily represent the official views of these agencies.

## References

- (1) D. T. Gillespie. A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. J. Comput. Phys., 22(4):403–434, 1976.
- (2) C. W. Gardiner, K. J. McNeil, D. F. Walls and I. S. Matheson Correlations in stochastic theories of chemical reactions. J. Stat. Phys., 14(4):307–331, 1976.
- (3) F. Baras and M. Malek Mansour Reaction-diffusion master equation: A comparison with microscopic simulations. Phys. Rev. E., 54(6):6139–6147, 1996.
- (4) M. J. Lawson, B. Drawert, M. Khammash, and L. Petzold. Spatial stochastic dynamics enable robust cell polarization. PLoS Comput. Biol., 9(7):e1003139, 2013.
- (5) M. Howard and A. D. Rutenberg. Pattern formation inside bacteria: fluctuations due to the low copy number of proteins. Phys. Rev. Lett., 90(12):128102, 2003.
- (6) D. Fange and J. Elf. Noise induced Min phenotypes in E. coli. PLoS Comput. Biol., 2(6):e80, 2006.
- (7) M. Sturrock, A. Hellander, A. Marsavinos, and M. Chaplain. Spatial stochastic modeling of the hes1 pathway: Intrinsic noise can explain heterogeneity in embryonic stem cell differentiation. J. Roy. Soc. Interface, 10(80):20120988, 2013.
- (8) M. Sturrock, A. Hellander, S. Aldakheel, L. Petzold, and M. Chaplain. The role of dimerisation and nuclear transport in the hes1 gene regulatory network. Bull. Math. Biol., 76(4):766–798, 2013.
- (9) J. Elf and M. Ehrenberg. Spontaneous separation of bi-stable biochemical systems into spatial domains of opposite phases. Syst. Biol., 1:230–236, 2004.
- (10) J. S. van Zon and P. R. ten Wolde. Green’s-function reaction dynamics: A particle-based approach for simulating biochemical networks in time and space. J. Chem. Phys., 123(23):234910, 2005.
- (11) J. S. van Zon and P. R. ten Wolde. Simulating biochemical networks at the particle level and in time and space: Green’s-function reaction dynamics. Phys. Rev. Lett., 94(12):128103, 2005.
- (12) S. Hellander and P. Lötstedt. Flexible single molecule simulation of reaction-diffusion processes. J. Comput. Phys., 230(10):3948–3965, 2011.
- (13) J. Hattne, D. Fange, and J. Elf. Stochastic reaction-diffusion simulation with MesoRD. Bioinformatics, 21(12):2923–2924, 2005.
- (14) B. Drawert, S. Engblom, and A. Hellander. URDME: a modular framework for stochastic simulation of reaction-transport processes in complex geometries. BMC Syst. Biol., 6(1):76, 2012.
- (15) I. Hepburn, W. Chen, S. Wils, and E. De Schutter. STEPS: efficient simulation of stochastic reaction-diffusion models in realistic morphologies. BMC Syst. Biol., 6(35):425–437, 2012.
- (16) D. Rossinelli, B. Bayati, and P. Koumoutsakos. Accelerated stochastic and hybrid method for spatial simulations of reaction-diffusion systems. Chem. Phys. Lett., 451:136–140, 2008.
- (17) T. T. Marquez-Lago and K. Burrage. Binomial tau-leap spatial stochastic simulation algorithm for applications in chemical kinetics. J. Chem. Phys., 127(10):104101, 2007.
- (18) S. Lampoudi, D. Gillespie, and L. R. Petzold. The multinomial simulation algorithm for discrete stochastic simulation of reaction-diffusion systems. J. Chem. Phys., 130(9):094104, 2009.
- (19) L. Ferm, A. Hellander, and P. Lötstedt. An adaptive algorithm for simulation of stochastic reaction-diffusion processes. J. Comput. Phys., 229(2):343–360, 2010.
- (20) A. Hellander and P. Lötstedt. Incorporating active transport in mesoscopic reaction-transport models inside living cells. Multiscale Model. Simul., 8(5):1691–1714, 2010.
- (21) B. Bayati. Fractional diffusion-reaction stochastic simulations. J. Chem. Phys., 138(10):104117, 2013.
- (22) S. A. Isaacson and C. S. Peskin. Incorporating diffusion in complex geometries into stochastic chemical kinetics simulations. SIAM J. Sci. Comput., 28(1):47–74, 2006.
- (23) S. Engblom, L. Ferm, A. Hellander, and P. Lötstedt. Simulation of stochastic reaction-diffusion processes on unstructured meshes. SIAM J. Sci. Comput., 31(3):1774–1797, 2009.
- (24) M. v. Smoluchowski. Versuch einer mathematischen Theorie der Koagulationskinetik kolloider Lösungen. Z. phys. Chemie, 92:129–168, 1917.
- (25) M. Tomita, K. Hashimoto, K. Takahashi, T. S. Shimizu, Y. Matsuzaki, F. Miyoshi, K. Saito, S. Tanida, K. Yugi, J. C. Venter, and C. A. Hutchison. E-cell: software environment for whole-cell simulation. Bioinformatics, 15(1):72–84, 1999.
- (26) S. S. Andrews, N. J. Addy, R. Brent, and A. P. Arkin. Detailed simulations of cell biology with Smoldyn 2.1. PLoS Comput. Biol., 6(3):e1000705, 2010.
- (27) F. C. Collins and G. E. Kimball. Diffusion-controlled reaction rates. J. Colloid Sci., 4(4), 1949.
- (28) D. T. Gillespie. A diffusional bimolecular propensity function. J. Chem. Phys., 131(16):164109, 2009.
- (29) S. A. Isaacson. The reaction-diffusion master equation as an asymptotic approximation of diffusion to a small target. SIAM J. Appl. Math., 70(1):77–111, 2009.
- (30) R. Erban and J. Chapman. Stochastic modelling of reaction-diffusion processes: algorithms for bimolecular reactions. Phys. Biol., 6(4):046001, 2009.
- (31) S. Hellander, A. Hellander, and L. R. Petzold. Reaction-diffusion master equation in the microscopic limit. Phys. Rev. E, 85(4):042901, 2012.
- (32) D. Fange, O. G. Berg, P. Sjöberg, and J. Elf. Stochastic reaction-diffusion kinetics in the microscopic limit. Proc. Natl. Acad. Sci. USA, 107(46):19820–19825, 2010.
- (33) A. Hellander, S. Hellander, and P. Lötstedt. Coupled mesoscopic and microscopic simulation of stochastic reaction-diffusion processes in mixed dimensions. Multiscale Model. Simul., 10(2):585–611, 2012.
- (34) I. C. Agbanusi and S. A. Isaacson. A comparison of bimolecular reaction models for stochastic reaction–diffusion systems. Bull. Math. Biol., 76(4):922–946, 2014.
- (35) S. S. Andrews and D. Bray. Stochastic simulation of chemical reactions with spatial resolution and single molecule detail. Phys. Biol., 1:137–151, 2004.
- (36) D. T. Gillespie. A rigorous derivation of the chemical master equation. Physica. A., 188:404–425, 1992.
- (37) C. W. Gardiner. Handbook of Stochastic Methods. Springer Series in Synergetics. Springer-Verlag, Berlin, 3rd edition, 2004.
- (38) D. T. Gillespie, A. Hellander, and L. R. Petzold. Perspective: Stochastic algorithms for chemical kinetics. J. Chem. Phys., 138(17):170901, 2013.
- (39) S. A. Isaacson. Relationship between the reaction-diffusion master equation and particle tracking models. J. Phys. A: Math. Theor., 41(6), 2008.
- (40) D. T. Gillespie, L. R. Petzold, and E. Seitaridou. Validity conditions for stochastic chemical kinetics in diffusion-limited systems. The Journal of Chemical Physics, 140(5):054111, 2014.
- (41) K. Takahashi, S. Tănase-Nicola, and P. R. ten Wolde. Spatio-temporal correlations can drastically change the response of a MAPK pathway. Proc. Natl. Acad. Sci. USA., 107(6):2473–2478, 2010.
- (42) S. A. Isaacson. A convergent reaction-diffusion master equation. J. Chem. Phys., 139(5):054101, 2013.