Reaction-diffusion kinetics on lattice at the microscopic scale

Lattice-based stochastic simulators are commonly used to study biological reaction-diffusion processes. Some of these schemes that are based on the reaction-diffusion master equation (RDME), can simulate for extended spatial and temporal scales but cannot directly account for the microscopic effects in the cell such as volume exclusion and diffusion-influenced reactions. Nonetheless, schemes based on the high-resolution microscopic lattice method (MLM) can directly simulate these effects by representing each finite-sized molecule explicitly as a random walker on fine lattice voxels. The theory and consistency of MLM in simulating diffusion-influenced reactions have not been clarified in detail. Here, we examine MLM in solving diffusion-influenced reactions in 3D space by employing the Spatiocyte simulation scheme. Applying the random walk theory, we construct the general theoretical framework underlying the method and obtain analytical expressions for the total rebinding probability and the effective reaction rate. By matching Collins-Kimball and lattice-based rate constants, we obtained the exact expressions to determine the reaction acceptance probability and voxel size. We found that the size of voxel should be about 2% larger than the molecule. MLM is validated by numerical simulations, showing good agreement with the off-lattice particle-based method, eGFRD. MLM run time is more than an order of magnitude faster than eGFRD when diffusing macromolecules with typical concentrations in the cell. MLM also showed good agreements with eGFRD and mean-field models in case studies of two basic motifs of intracellular signaling, the protein production-degradation process and the dual phosphorylation cycle. Moreover, when a reaction compartment is populated with volume-excluding obstacles, MLM captures the non-classical reaction kinetics caused by anomalous diffusion of reacting molecules.

On the other hand, in lattice approaches, the average diffusion behavior is adopted and the reactions follow either the simple first-order process, or the second-order process when two reactive molecules meet on the same lattice voxel. Such approaches reduce the computational cost even in crowded space and provide an efficient way to simulate large numbers of molecules and reactions. Within lattice approaches, variation exists depending on how each molecule is represented and reaction is modeled. In the well-established reaction-diffusion master equation (RDME) models [27][28][29][30][31][32][33], space is discretized into lattice voxels called subvolumes. In each subvolume, point-like molecules are assumed to be dilute and well-mixed. To obey the wellmixed condition, there is a limit to the size of the subvolume [28,34,35], which in turn imposes a limit to the spatial resolution. Diffusion of molecules across subvolumes is modeled as a firstorder reaction with a concentration dependent rate. Unimolecular and bimolecular reactions only occur within each subvolume with a rate defined by the propensity function [36]. Compared to continuum-based schemes, RDME models RD from the mesoscopic to the macroscopic scale but not at the microscopic scale. However, there have been several efforts to overcome the wellmixed limit in RDME models and to bridge mesoscopic and microscopic scales [34,35,37,38].
Apart from the RDME lattice models, there is another class of schemes, which we refer as microscopic lattice method (MLM) that represents molecules at single particle resolution [39][40][41][42][43][44][45][46][47][48][49][50]. In most of these schemes [39, 41-45, 47, 50], the size of the voxel follows the molecule size, whereas in the small-voxel tracking algorithm (SVTA) [48], a particle can occupy multiple voxels, providing greater spatial resolution at the cost of higher computational complexity. In MLM, a molecule hops into a neighbor voxel at a constant rate such that normal diffusion is satisfied. Excluded volume arises naturally since the size of molecule is directly reflected by the voxel size and occupancy in the lattice. Similar to RDME models, unimolecular reaction is modeled as a first-order process. Bimolecular reactions are coupled to molecular collisions in all of these schemes except GridCell [45]. In the collision-based reaction schemes, the steady-state reaction rate follows the macroscopic effective reaction rate when the reaction is activationlimited. However, the reaction accuracy of MLM has not been studied in detail when it is diffusion-influenced. In a recent work, Sturrock [51] also reported several shortcomings in MLM, notably in the accuracy of Spatiocyte [47] when estimating steady-state bimolecular reaction rates.
Our focus in this work is to examine in detail the accuracy and consistency of MLM in solving diffusion-influenced reactions using theoretical analysis and numerical simulations. The theoretical framework here is constructed based on the hexagonal close-packed (HCP) lattice but is also applicable for any regular lattice arrangements such as the simple cubic lattice. We employ the Spatiocyte scheme to construct and analyze the general theoretical framework of MLM in both activation-and diffusion-limited regimes. We then describe the first-passage behavior of the method according to the random walk theory and obtain the analytical formula for the total rebinding probability of a pair of reacting molecules and their effective reaction rate constant. Next, we perform numerical simulations to evaluate the accuracy of the theory and investigate the time-dependent kinetics. We found that MLM exhibits the expected steady-state and asymptotic time-dependent behaviors of the reaction as in the collision-based continuum model. Subsequently, we evaluate the performance of MLM in comparison to other well-known off-lattice methods. As application examples, we show that the method correctly recapitulates the time-dependent behavior of proteins in the production-degradation process and the dual phosphorylation-dephosphorylation cycle, two fundamental building blocks of intracellular signaling. Finally, we demonstrate the effects of crowding obstacles on the kinetics of a simple bimolecular reaction with MLM.

II METHODS
We begin by presenting the theoretical background of the Collins-Kimball [25] approach in modeling irreversible bimolecular diffusion-influenced reaction. We highlight the particle-pair formalism for the reaction rate, which will be used in the MLM theory. We then briefly describe the Spatiocyte RD scheme, an MLM implemented on the HCP lattice. Finally, we construct the theoretical framework of reaction rate coefficient on lattice using the Spatiocyte scheme.

A Irreversible bimolecular diffusion-influenced reaction in continuumbased framework
Consider an irreversible bimolecular reaction involving two distinct species: where A and B are hard-sphere molecules with radii r A and r B , respectively. The molecules diffuse in three-dimensional (3D) space with diffusion coefficients D A and D B . The time evolution of the species concentration is well-described by a time-dependent rate coefficient k irr (t): Smoluchowski [24] derived the rate coefficient by relating the diffusion coefficient of the molecules with molecular collisions, leading to the product formation. In his work, A is made up of an immobile molecule and is surrounded by multiple diffusing B molecules. Collins and Kimball extended the Smoluchowski theory by modeling the reaction using radiation boundary condition and obtained the rate coefficient in 3D space as a function of microscopic parameters, namely the microscopic or the intrinsic reaction rate constant, k a , the contact distance of the reacting molecules, R = r A + r B and the relative diffusion coefficient of the pair, D = D A + D B : Here, k D = 4πRD is the collision rate, Φ(x) = exp(x 2 )erfc(x) and τ = (1/D)(k a R/(k a + k D )) 2 . The rate coefficient (3) starts (t = 0) at k a but decays rapidly to [52] k irr (t) k ef f 1 + k a k a + k D at long-time. k ef f is the steady-state or the effective reaction rate constant given by [53]: According to Noyes theory [54][55][56], the rate coefficient can be expressed equivalently using the particle-pair approach: k irr (t) = k a S(t; R), where S(t; R) denotes the survival probability of an isolated reactant pair at time t given that they were initially in contact. Additionally, let p reb (R, t|R, 0) denote the rebinding-time probability distribution for a reactive particle-pair separated by distance R at time t, given that the pair were initially in contact. In the case of radiation boundary condition, the probability distribution is given by (see Appendix A) with τ = tD(1 + k a /k D ) 2 /R 2 . Note that the survival probability S(t; R) is the same as the probability that the first rebinding event between an initially in-contact pair has not yet occurred at time t. Hence we can rewrite Eq. (6) as At long-time, we then have where the integrated term gives the total rebinding probability: Therefore, the effective rate constant (5) can also be written in terms of the total rebinding probability: The above relation was also described previously, but in the context of irreversible and reversible rate constants [57]. In subsequent sections, we use the relations described by Eqs. (8) and (11) as the central concepts to derive the rate coefficient in MLM. Figure 1: A voxel on the HCP lattice has 12 nearest neighbor voxels. The distance between the centers of two adjacent voxels is the voxel size l.

B Spatiocyte reaction-diffusion scheme
In the Spatiocyte scheme (see Algorithm 1) [47], space is discretized into HCP lattice because the arrangement allows the highest density of regular sphere voxels in 3D space. The voxel has a diameter l and can be occupied by at most, a single molecule. At each diffusion time step t d , a molecule can hop to one of its 12 nearest neighbor voxels (see Figure 1) with the step acceptance probability P w = 1.
Given the irreversible bimolecular reaction in Eq. (1), a collision arises when B meets A at the destination voxel. The collision is reactive with an acceptance probability P a = ∆N C /Z.
is the microscopic change in the number of product molecules in step interval t d and Z is the expected number of collisions between A and B in the interval (see Appendix B). The acceptance probability can then be expressed as [47] The above relation is applicable when the reaction is activation-limited (k a k D ). For diffusion-influenced reactions (k a k D ), the collision rate Z is reduced relative to the production rate ∆N C . The acceptance probability P a = ∆N C /Z would then have the issue of exceeding unity when ∆N C > Z. The Spatiocyte scheme overcomes this issue by reducing the simulation interval by a factor of α to t = t d α. With the reduced interval, the effective number of collisions in t d , Z is increased. The step and reaction acceptance probabilities are then decreased accordingly to P w = α and P a = P a α, respectively. Algorithm 1 describes how α is set. In summary, the Spatiocyte scheme operates with α = 1 when P a ≤ 1 (activation-limited case) and with α < 1 when P a > 1 (diffusion-influenced case).
where xy denotes the pair of reactive species x and y; reaction acceptance probability P axy = P axy α; step acceptance probability P wx = α; end Main loop: while S = {} and t sim < t end do t sim ← τ x = next event in S; get species identity x; get current voxel location s 0 ; reschedule next event, τ x = τ x + t for each molecule of species x do choose a random target voxel s 1 ∈ {nearest neighbor of s 0 }; if s 1 is vacant then draw rand; if rand ≤ P wx then walk succeeded, s 0 ← s 1 ; else walk rejected, s 0 ← s 0 ; end else if s 1 contains reactant species y then draw rand; if rand ≤ P axy then reaction xy accepted, s 0 ← s 1 else reaction failed and walk rejected, s 0 ← s 0 end else walk rejected, s 0 ← s 0 ; end end end Algorithm 1: Basic outline of the Spatiocyte algorithm for bimolecular reactions. t sim is the current simulation time, t end -t sim is the simulation duration, P axy is the reaction acceptance probability for the reactive pair of species x and y, t d = l 2 /6D x is the diffusion (hopping) time step of the current species x, l is the voxel size, D x is the diffusion coefficient of x, and rand is a random number drawn from the uniform distribution with the interval [0, 1).

C Rebinding probability and reaction rate on HCP lattice
As an alternative to the time-dependent reaction rate coefficient in Eq. (8), we define a discretespace version with a step-dependent rate coefficient on lattice as [39,58] where m is the simulation step, which is related to the simulation time by 6D x t = ml 2 , k a is the initial reaction rate constant on lattice (see Appendix B) and H n is the lattice analogue of the rebinding-time probability function p reb (R, t|R, 0) in diffusion step n. At long-time, the effective rate on lattice follows similarly to Eq. (9): where the summation term (14) corresponds to the total rebinding probability on lattice.
To obtain the analytical expression for H n , we consider again a reactive pair A and B, which are initially in-contact by occupying adjacent voxels on lattice. We are interested in the rebinding-time probability distribution as a function of the diffusion step n. Without losing generality, we can fix one of the molecules and diffuse the other with the relative diffusion coefficient D. Then, the rebinding-time probability distribution of A and B is related to the arrival-time probability distribution of a random walker to the origin for the first time, given that the walk started at one of the neighbor voxels of the origin with diffusion coefficient D. In the following sections, we define H n explicitly and use it to derive the rate coefficient on HCP lattice. Since the approaches for activation-limited and diffusion-influenced cases are different in the Spatiocyte scheme, we perform their derivations separately.
1 Activation-limited case (k a k D , α = 1) We denote s 0 as the voxel at origin, s 1 as an element of the set of immediate neighbor voxels of s 0 . We define F n (s a |s b ) as the first-passage time distribution for a random walker to walk from voxel s b to s a , that is, the probability of arriving at voxel s a for the first time at the nth step, given that the walk started at voxel s b . We first consider the rebinding-time probability distribution for the case P a = 1. Let F n (s 0 |s 0 ) and F n (s 0 |s 1 ) denote the first-passage time distributions to origin from origin and s 1 , respectively. The two probabilities are related via where p(s 0 → s 1 ) = 1 is the transition probability from s 0 to s 1 in a single step. This implies that the trajectory we are interested in, which is from an in-contact situation (e.g., A at s 1 and B at s 0 ) to the rebinding situation (A hops to s 0 ) in a single step, is equivalent to the 2-step trajectory, s 0 → s 1 → s 0 .
Therefore, the rebinding-time probability distribution is fully described by F n (s 0 |s 1 ) and is related to F n (s 0 |s 0 ). The latter can be obtained analytically from its probability generating function F (s 0 |s 0 ; z) = ∞ n=0 F n (s 0 |s 0 )z n [59] (see Appendix D).
As for P a ≤ 1, the trajectories that have undergone failed reaction attempts before step n are included in the rebinding-time probability distribution: where F j n (s 0 |s 0 ) is the probability to reach the origin for the jth time at the nth step (I.1.9 in [60]): where F 1 n (s 0 |s 0 ) = F n (s 0 |s 0 ). The generating function of H n (s 0 |s 1 ) in terms of F (s 0 |s 0 ; z) is (see Appendix A): By taking the limit z → 1 on H(s 0 |s 1 ; z), we obtain the total rebinding probability on lattice as where F (1) = F (s 0 |s 0 ; z = 1). It was shown previously that the probability generating function of the HCP lattice is topologically equivalent to that of the face-centered cubic (FCC) lattice [61]. Therefore, we have F (1) ≈ 0.256318 ( [62], p. 153) for HCP lattice. Finally, if we set the initial rate k a = 3 √ 2lDP a (see Appendix B) and substitute the total rebinding probability H reb from Eq. (19) into Eq. (14), we obtain the effective rate constant on lattice as 2 Diffusion-influenced case (k a k D , α < 1) The rebinding-time probability distribution G n (s 0 |s 1 ) of the diffusion-influenced scheme is defined as G n+1 (s 0 |s 1 ) = S n (s 1 |s 1 ) p(s 1 → s 0 ), for n ∈ N, where is the probability for a successful reaction, P 1 (s 0 |s 1 ) is the probability to select s 1 given that the molecule is in s 0 (= 1 12 for HCP lattice), and S n (s 1 |s 1 ) is the probability that a particle is in contact after n-steps (see Appendix B for more details). The probability generating function of G n (s 0 |s 1 ) on HCP lattice is given by (see Appendix B) G(s 0 |s 1 ; z) = P a α/12 where S(s 1 |s 1 ; z) is the probability generating function of S n (s 1 |s 1 ) .
Taking the limit z → 1, we get the total rebinding probability as (Appendix B ) which is identical to Eq. (19) in the activation-limited case. Similarly, by substituting the summation term in Eq. (14) with Eq. (24), we get the effective rate constant for the diffusioninfluenced case as which also follows Eq. (20). Henceforth, we adopt the same notations of the effective reaction rate and total rebinding probability for both the activation-limited and diffusion-influenced cases.

D Comparison with continuum-based theory
Since the effective rate on lattice (25) has the same form of Eq. (11) in continuum, we can match them by equating the initial rate and total rebinding probability of the two: k a = k a and G reb = P reb . With the former relation, the reaction acceptance probability is connected to the initial rate constant, diffusion coefficient, and voxel size by Employing the G reb = P reb relation, the voxel size is found to be about 2% greater than the molecule size: The Spatiocyte scheme is thus guaranteed to have the same effective rate and total rebinding probability as the continuum framework provided that Eqs. (26) and (27) are satisfied. In addition, the expression for lattice effective rate constant follows the same form of the continuum-based framework: where k D = 3 √ 2lD(1/F (1) − 1). According to Eq. (27), accurate matching of both the effective rate and the total rebinding probability requires the voxel size to be larger than the molecule size. Nonetheless, during modeling we can fix the voxel size to be the same as the molecule size, l = R. In this case, it is still possible to match the lattice effective reaction rate to the continuum-based rate by setting the reaction acceptance probability to However, this is done at the expense of losing accuracy in the total rebinding probability, since G reb = P reb .
For the reversible reaction, A + B ka k d C, the local detailed balance on lattice is achieved by choosing a lattice dissociation rate constant k d from the following equilibrium constant relation: The MLM method can simulate the dissociation reaction as a first-order process with rate k d and place the dissociated pair of molecules at an in-contact condition.

E Numerical simulations
We verify the main theoretical results presented above with numerical simulations using Spatiocyte. Spatiocyte is included in E-Cell System version 4 [63], an open-source biochemical simulation environment that supports multiple algorithms, time scales and spatial representations.

III RESULTS AND DISCUSSION
We first validate the theory of total rebinding probability and its time-dependent behavior on lattice using numerical simulations. We examine the accuracy of the reaction rate coefficient and its time-dependent behavior on lattice. We then compare the diffusion and reaction performances of MLM and several other off-lattice particle methods. Finally, we evaluate MLM in protein production-degradation process, dual phosphorylation cycle and a simple bimolecular reaction in a crowded compartment.

A Numerical validation of MLM theory 1 Rebinding probability
We examine the rebinding probability distribution of a reactive pair, A and B that are initially in-contact. The theoretical rebinding-time probability distribution H n (s 0 |s 1 ) and G n (s 0 |s 1 ) are validated against numerical results. In the activation-limited case (k a /k D 1), the expected first rebinding probability at nth step is obtained using Eq. (16), whereas in the diffusion-influenced case (k a /k D ≥ 1), the probability is calculated from the generating function G(s 0 |s 1 ; z) Since the theoretical rebinding-time probability distribution on lattice is validated by simulations, the analytical formulas for the total rebinding probability derived from it, Eqs. (19) and (24) are therefore valid.
To illustrate the dependency of total rebinding probability on k a /k D , we obtained the probability at various k a /k D up to n = 10. Table (2) shows the simulated and the expected theoretical Table 1: Theoretical and simulated rebinding-time probabilities on lattice for activation-limited and diffusion-influenced cases. Simulation parameters: l = 0.01 µm, volume = (100 l) 3 with periodic boundary, runs = 1 × 10 9 .  Table 2: Theoretical and simulated rebinding probabilities up to n = 10 steps with k a /k D ratios ranging from the highly activation-limited case (k a /k D = 0.01) to the strongly diffusion-influenced case (k a /k D = 100). Simulation parameters: l = 0.01 µm, D = 1 µm 2 /s, volume = (10000 l) 3 with periodic boundary, runs = 1 × 10 9 . α = 1/P a for diffusion-influenced cases. values for various k a /k D ratios. Both simulated and theoretical values coincide well, with discrepancies never exceeding 0.03%. Qualitatively, the total rebinding probability increases with larger k a /k D values, consistent with the continuum theory (10). We then evaluated the rebinding-time probability distribution by recording the time taken for A and B to associate immediately after a dissociation event. We performed the simulations for a large number of steps and independent runs. Figure 2 shows the average number of rebinding events per unit time at k a /k D = 0.1, 1 and 100. Lines depicting the rebinding-time probability distribution of the continuum-based model according to Eq. (7) are also shown as reference. It is clear that at times larger than t d , the time-dependent behavior of lattice simulations is consistent with the continuum-based model. The scaling behavior at long-time, p reb (t) ∝ t −3/2 is a well-known characteristic of a 3D random walker returning to the origin [64]. We have corroborated this result with detailed asymptotic analysis that is provided in Appendix 1.
Note that in the diffusion-influenced case (k a /k D = 1 and 100), finer step intervals generate rebindings at times smaller than the diffusion time step t d , denoted by the vertical dashed line in Figure 2. In this temporal regime, MLM behaves differently from the continuum-based framework because the MLM reaction kinetics approximates the Poisson process (see Appendix 3). Despite the difference, the rebinding behavior correctly converges to the continuum-based formalism for times larger than t d .

Reaction rate
We evaluated the accuracy of the effective reaction rate constant for irreversible bimolecular reactions (1) over various k a /k D regimes on lattice. We considered an immobile species A and a diffusing species B that are uniformly distributed at initialization with concentrations Rebinding time probability [A] and [B], respectively. We recorded the surviving fraction of A molecules at each time step. Figure 3a displays the survival probability of A and the expected theoretical curve [65]). From the survival probability, we calculated the time-dependent reaction rate coefficient using (Eq. 2.1 in [65]) We adopted the following discretization scheme for the time derivative to get the discrete rate coefficient: where n is the index of the discretized S A and t. The boundary cases are computed as where N denotes the final time step. The reaction rate coefficient obtained for various k a /k D ratios are shown in Figure 3b along with their corresponding theoretical curves from Eq. (4).
Recall that the long-time asymptotic variant of the Collins-Kimball theory (4) has the form where C 1 and C 2 denote the steady-state rate constant and the time-dependent term, respectively. We fitted Eq. (35) to the numerical data, omitting early time points to avoid nonsteady-state effects. The resulting C 1 and C 2 parameters after fitting are listed in Table 3.   The theoretical values correspond favorably to the estimated steady-state reaction rate constants and are well within the standard error, thus validating the lattice theory for the effective rate. The time-dependent terms are also in good agreement with the theory, especially in the diffusion-limited case, with discrepancy less than 1%. This is consistent with the asymptotic analysis carried out in Appendix 2 and 4. In the activation-limited case (k a /k D = 0.1), the fitted C 2 had the largest deviation from theory because the standard error was also the highest. The low number of data points contributed to the higher standard error. Nonetheless, we did not increase the data points because C 2 has a weaker influence in activation-limited reactions than C 1 .

B Performance 1 Diffusion
We compared the 3D diffusion performance of MLM using Spatiocyte (git 9757fb3) and three other off-lattice particle-based simulation methods, Smoldyn [66] (version 2.55), eGFRD [14] (in E-Cell System version 4.1.4) and fast Brownian dynamics (BD) [67] (C++ program example in Spatiocyte git 9757fb3). When the molecules are represented as hard-spheres with volume exclusion, Spatiocyte required shorter run times than Smoldyn in all cases ( Figure 4a). Spatiocyte achieves comparable or better performance than eGFRD in the typical concentration range of cytoplasmic macromolecules (0.1 to 10 µM). For example at 6 µM in volume 30 µm 3 , Spatiocyte is about 4.5 and 16 times faster than Smoldyn and eGFRD, respectively. In contrast to eGFRD, Spatiocyte and Smoldyn execution times increase with the number of molecules but not the molecular crowdedness (V = 30 µm 3 vs. 3 µm 3 ). The simulation times of Spatiocyte scale almost linearly with the number of molecules (T ∝ N ), which is not apparent Table 3: The steady-state rate constant C 1 and the time-dependent term C 2 of reaction (1) at various k a /k D were obtained by fitting the simulated reaction rate coefficient with Eq. (35). Uncertainty in the simulated data was used as a weight in the fitting. Theoretical values from Eq. (4) are listed for comparison. Simulation parameters: l = 0.01 × 1.0209 µm, D A = 0, D B = 1 µm 2 s −1 , volume = (350 l) 3 with periodic boundary, N A = N B = 4000, duration = 0.05 s, runs = 3 × 10 4 , α = 1/P a for the diffusion-influenced case. with Smoldyn and eGFRD. The drastic slowdown of eGFRD at higher concentrations is caused by the shorter time steps required to resolve many molecular interactions that take place in the densely occupied system [14]. If molecules are represented as dimensionless point particles, higher diffusion performance is expected since inter-molecular collisions can be ignored. Figure 4b shows the run times of Spatiocyte, Smoldyn and fast BD when diffusing point particles with the same simulation interval. eGFRD was not considered here since it only supports molecules with physical volume. Spatiocyte and fast BD execution times showed an almost linear scaling with the number of molecules. Although Smoldyn did not scale as well, it had the fastest run times when the number of diffusing molecules was 30,000 or less. Spatiocyte outperformed fast BD in all tests and is on average 2.5 times faster. As expected, the simulation times of all three methods were not affected by the crowdedness in the volume since molecular collisions are disregarded. On average, Spatiocyte takes about 2 times longer to diffuse hard-sphere molecules than point particles.

C Application Examples
We applied MLM to model two fundamental RD systems of intracellular signaling, the productiondegradation process, previously studied using lattice-based methods [51,70,71], and the dual phosphorylation-dephosphorylation cycle of the mitogen-activated protein kinase (MAPK) cascade [72][73][74], a common motif found in signal transduction systems but with a response function that is highly sensitive to the binding kinetics. We also report the effects of excluded volume on the kinetics of a simple bimolecular reaction using MLM.

Production-degradation process
Consider the production and degradation processes of protein A represented by a zero-order production coupled with a second-order degradation: The concentration of A will go through an initial transient state before settling down at a steadystate equilibrium, [A] = k 1 /(k 2 [B]) that fluctuates according to the Poisson distribution [70].
To confirm if MLM can recapitulate the production-degradation process correctly in 3D space, Molecule or voxel radius (r), simulation or diffusion step interval (∆t) and run time (T ) are as indicated. All simulations were executed on the same workstation with Intel Xeon X5680 3.33 GHz CPU, 48 GB memory and Ubuntu 16.04 LTS operating system.
we have simulated the process with Spatiocyte and compared the outcomes with eGFRD and the well-mixed model. To generate the results of the well-mixed model, we solved the rate equation using an ordinary differential equation (ODE) solver. The time-series of A is shown in Figure 6a, while the equilibrium values are provided in Table 4. As evident from the figure and table, Spatiocyte results are all in good agreement with both the well-mixed model and eGFRD.
Recently, the Spatiocyte scheme was reported to not only fail to reproduce the expected equilibrium value of A but also generate different values depending on the voxel size [51]. In the report, the effective bimolecular rate k 2 was used in the calculation of reaction acceptance probability instead of the intrinsic reaction rate constant k a , which inevitably caused the deviation from the well-mixed model (see first row of Table 4). As shown in Figure 6a and Table 4, there was no discrepancy when the intrinsic rate k a was used to compute the reaction-acceptanceprobability (26). The well-known relation between k a and k 2 is given by Eq. (5), wherein k 2 is represented by k ef f . Furthermore, just as in the well-mixed and eGFRD models, the resulting equilibrium concentration from Spatiocyte is also independent of the molecule radius or spatial discretization. Conversely, the RDME method deviated substantially from the well-mixed result when the voxel size is small, which is expected [34,37,38].
The well-mixed model assumes the time scale of diffusion to be always shorter than that of the reactions. As a result, molecules are expected to be uniformly distributed at all times and reactions can take place independent of spatial localization. The well-mixed assumption is valid when describing activation-limited reactions but when they are diffusion-influenced, the position of molecules should be taken into account. We therefore expected some disparity between the well-mixed model and MLM when the production-degradation process is diffusioninfluenced. In Figure 6b, at smaller diffusion coefficients (D x = 0.01, 0.02), the equilibrium concentrations are indeed lower with Spatiocyte than with well-mixed model. Spatiocyte behavior is consistent with eGFRD, which also accounts for molecule positions. RDME however,  has the same outcomes as the well-mixed model. The reduction in equilibrium value when the diffusion coefficient is decreased was previously Table 4: Equilibrium concentration of A in Eq. (36) simulated with Spatiocyte and eGFRD at different spatial discretizations. k 2 is the effective rate, k a is the intrinsic rate, l is the voxel size, K = 2 1/6 L/l is the compartment length in number of voxels, while L denotes the actual length [51]. At l = 0.01, K = 521; at l = 0.04, K = 130; and at l = 0.1, K = 52. The well-mixed equilibrium concentration is 5 µm  (1) described by the microscopic theory of Agmon and Szabo [52]. In contrast to the Collins-Kimball theory, Agmon and Szabo have considered the non-negligible effect of B concentration on the effective reaction rate, especially when the reaction is diffusion-influenced. The slow diffusion of molecules increases the effective contact radius, resulting in higher effective annihilation rate (see Appendix E for a detailed argument). The output of the production-degradation process according to the microscopic theory is shown in Figure 6b as a dashed line that coincides with Spatiocyte and eGFRD, further verifying the MLM theory. Given the same diffusion and macroscopic reaction rates, the change in the Spatiocyte voxel size does not affect the equilibrium behavior (at r = 0.1 and r = 0.05 in Figure 6b) since the reaction acceptance probability, P a is adjusted according to the voxel size to obtain the correct macroscopic behavior.
On the other hand, RDME shows large deviation from the expected values at slow diffusion. The inability of conventional rate equation and RDME to correctly capture diffusion-influenced reactions has previously been noted and worked on before [34,37,38,70,75]. By incorporating the diffusion coefficient into the bimolecular reaction propensity formula (Eq. 26 in [70]), the equilibrium concentration of RDME shows a better agreement with the expected values (see RDMEm, r = 1.0 in Figure 6b). However, when the reaction is diffusion-limited (D x = 0.01, 0.02), unlike MLM, the subvolume size of RDMEm cannot reach the microsopic resolution, r = 0.05. This is because the size is constrained by a critical value (Eq. 25 in [70]) that preserves the well-mixed condition. At D x = 0.01 for example, the critical subvolume size is about 13 times the molecule diameter, any size smaller is invalid.
We have also examined the fluctuation of A at equilibrium, as depicted in Figure 6c. At D x = 0.1, the histogram of Spatiocyte matches the distribution curves of eGFRD and the well-mixed model (Gillespie method [36]). At much reduced diffusion coefficient (D x = 0.02) however, both Spatiocyte and eGFRD shared similar distributions, with the width becoming narrower and the mean value shifting to the left. With the modified propensity function, RDMEm also exhibited similar distribution. The narrow width and the shifted mean are consistent with the characteristics of the Poisson distribution.
It was reported that MLM would not be able to solve the first-order production-degradation A accurately because of its spatial discretization scheme [51]. When the number of total voxels in the compartment, N v is less than k 1 /k 2 , the equilibrium concentration deviates from the well-mixed model. This deviation however, is a direct consequence of the volume exclusion property of MLM. Since each voxel can only occupy a single molecule, there would be an insufficient number of vacant voxels to accommodate new molecules when the degradation rate is not sufficiently fast. The maximum occupancy on HCP lattice simply reflects the maximum physical occupancy of voxel-sized molecules in the compartment because the HCP arrangement packs the highest density of sphere voxels [76]. Just as in the cellular compartment, no more molecules can be added into the system when the number of generated molecules exceeds available free space. Moreover, since only about 34 % of the cell volume is occupied by macromolecules [77], it is also an unlikely scenario to fully occupy the voxels of HCP lattice with macromolecules. With the multi-algorithm implementation of Spatiocyte [78], we can use the Gillespie's Next-Reaction method [79] to simulate small molecules that are in large abundance and are homogeneously distributed. In this case, the equilibrium result is independent of spatial discretization since the method assumes the well-mixed condition.

Dual phosphorylation-dephosphorylation cycle
In mean-field models, the spatio-temporal correlation of microscopic rebinding events is not resolved explicitly because the correlation usually does not cause a significant impact on the dynamics at the macroscopic scale. One case where the correlation does influence the macroscopic response is the dual phosphorylation-dephosphorylation cycle of the MAPK cascade [72][73][74], shown in Figure 7a. The substrate MAPK (K in Figure 7a) is phosphorylated in a two-step process by the MAPK kinase (KK) and dephosphorylated by a phosphatase P. The phosphorylation and dephosphorylation processes proceed according to the Michaelis-Menten kinetics and exhibit distributive property [73], wherein the enzymes must unbind from the substrate before they can rebind and modify the second site. Upon phosphorylation or dephosphorylation, the respective enzymes are inactivated (denoted as KK* and P*), and reactivated (KK or P) after some time τ rel . When the reactivation time is short and the enzyme-substrate reaction is diffusion-limited, the newly dissociated enzyme and substrate are close enough to rebind instead of escaping into the bulk. These microscopic rebinding events alter the response sensitivity of the phosphorylation state as shown by Takahashi et al. [14] using eGFRD. Processive behavior caused by rebindings of the same enzyme results in higher overall phosphorylation rate than the distributive case where the dissociated molecules can escape rebinding [73,74]. Such microscopic spatio-temporal correlation has been shown to change the response sensitivity of the phosphorylation state, which can cause the subsequent removal of ultra-sensitivity or bi-stability in the system [14,80]. Rebinding events taking place within very short time scales are difficult to be captured by RDME because of the fine spatial resolution required. To test whether MLM can resolve such events faithfully, we use Spatiocyte to model the dual phosphorylation cycle with the same parameters from [14]. Distributive and processive models are represented by Eqs. 1-5 of [14], and were solved using ODE solver. Figure 7b displays the steady-state response curves of Spatiocyte and reference theoretical models. Note that since the reactivation time τ rel is equal to or less than the diffusion time step t d (given in Figure 7b), the molecules can rebind soon after dissociation. The simulation result coincides very well with the switch-like response curve of the distributive model at fast diffusion (D x = 4 µm 2 s −1 ), whereas at much slower diffusion (D x = 0.06 µm 2 s −1 ), it converges to the graded response curve of the processive model. The influence of diffusion on the response curve can be understood through the rebinding events. When diffusion is slow, reactions become more diffusion-limited and rebinding occurs at higher frequency. The ensuing processive-like mechanism then leads to the loss of the switch-like response curve. Conversely, in the limit of fast diffusion as assumed in the mean-field model, a sharper switch-like response curve is recovered because of fewer rebindings.  The parameter ranges examined so far have a stable steady-state as demonstrated by the response curves in Figure 7b. When the total concentration of the substrate is increased fivefold, the mean-field theory generates hysteresis, shown by the dotted and dash-dotted lines. The dotted line represents the response when initialized with [K pp ]/[K] total = 1, whereas the dash-dotted line has the initial condition [K pp ]/[K] total = 0. MLM produced similar responses when the diffusion is fast (D x = 4) (diamond markers in Figure 7b). However, as diffusion slowed down to D x = 0.06, the bistability is lost (triangle markers). Bistable states appear when the diffusion is fast and the substrate concentration relative to enzyme is high. For example, at the initial state when almost all substrates are in the unphosphorylated form, most kinase will be bound to the substrates rapidly. Hence, a substrate that has been phosphorylated once is more likely to be dephosphorylated by free phosphatase than to be phosphorylated the second time by scarce and fast diffusing kinase. The inverse situation where all substrates are in the phosphorylated form would also respond similarly to phosphatase. On the other hand, when diffusion is slow, the kinase activity becomes processive because of the high rebinding probability. As a result, molecules are more likely to be phosphorylated or dephosphorylated consecutively before they could be disrupted by antagonistic enzymes from the bulk. This example highlights how local spatio-temporal correlation can change the binding behavior and results in a different global response than the one predicted by the mean-field model.
As a side remark, in the original Spatiocyte scheme [47], the voxel adopts the size of the diffusing molecules. However, as we found in the Methods section, the voxel needs to be about 2% larger than the molecule size (27) for the total rebinding probability and the effective rate constant to be exactly the same as in the continuum-based theory. Despite the 2% difference in voxel sizes, both new and original schemes displayed very good fit with the expected dual phosphorylation cycle response curves in Figure 7b. To be fully consistent with the continuumbased theory however, the size should be set according to Eq. (27). The voxel size is not hard-coded to be the same as the molecule size and can be easily specified in the Spatiocyte model file [78].

Effects of excluded volume on bimolecular reaction
Excluded volume in the cell arising from crowded obstacles such as macromolecules, Golgi apparatus or cytoskeletal elements can cause anomalous diffusion of reacting molecules [41,81]. Anomalous diffusion has been shown to generate non-classical reaction kinetics on 2D [42,43] and 3D lattices [49]. Here, we use MLM on HCP lattice to examine the effects of volume exclusion on the bimolecular reaction E + S − → ∅ in the presence of uniformly distributed immobile obstacles. E and S have the radius 5 nm and diffusion coefficient, Hence, D 0 is the diffusion coefficient in non-crowded dilute condition. Bimolecular intrinsic reaction rate constant k a = 10k D is chosen such that the reaction is diffusion-limited. We first consider the effects of immobile obstacles on diffusing molecules. We calculate the time-dependent diffusion coefficient from the mean-squared displacement of simulated particle trajectories. The time-dependent diffusion coefficient in Figure 8a indicates that the diffusion is anomalous at short times and normal at long times. The crossover time from anomalous to normal diffusion depends on the volume occupancy. The reduced long-time diffusion coefficient is well-described by [81,82], where φ p ≈ 0.77 is the percolation threshold for HCP lattice. We confirmed that the long-time diffusion coefficients obtained for φ in Figure 8a (dashed lines) are consistent with D in Eq. (37). Figure 8b shows that the survival probability of E decays slower when the volume occupancy, φ is increased. From the survival probability, we can calculate the rate coefficient according to Eq. (33) to obtain the kinetics. We replaced the constant concentration term [B] (33) with the time varying term [E](t) in the equation. For the dilute case (φ = 0) in Figure 8c, there is a good agreement for the simulated k(t) with the Collins-Kimball rate coefficient (3). As φ increases to 0.3 and 0.5, the overall reaction rate decreases, and thus progressively diverges from the Collins-Kimball rate. Despite the discrepancy, the rates can still conform to the Collins-Kimball theory when the long-time diffusion coefficient (37) is used.
As the volume occupancy approaches the percolation threshold (Figure 8c, φ = 0.7), the kinetics begins to deviate from the Collins-Kimball theory. The deviation is strongest at φ = 0.8, which is beyond the percolation threshold. Note that at lower volume occupany (φ = 0.3, 0.5), the anomalous to normal diffusion crossover time in Figure 8a is faster than the observation time in Figure 8c. Here, the kinetics is well described by the long-time effective diffusion coefficient. However, when the crossover time is comparable to the observation time because of the increased volume occupancy (Figure 8a, φ = 0.7), the effects of anomalous diffusion is visible in the kinetics (Figure 8c, φ = 0.7). At above the percolation threshold (φ = 0.8), anomalous diffusion does not crossover to normal diffusion. As a result, the long-time diffusion coefficient eventually decays to zero. In these highly crowded cases, the Collins-Kimball theory fails to describe the kinetics.
Grima and Schnell [44] have shown that reaction kinetics, either classical or non-classical, is not determined by the heterogeneity of the accessible space but rather by the reaction probability and the initial condition. In the Smoluchowski and Collins-Kimball framework, reaction follows classical kinetics when it is activation-limited (k a /k D 1) but non-classical kinetics is observed when it is diffusion-influenced (k a /k D 1). The non-classical behavior in the latter is well-described by Eq. (3) using microscopic parameters. The corresponding long-time behavior up to the second order term scales according to Eq. (4), which has the same general form of the Zip-Mandelbrot equation proposed by Schnell and Turner [43,49]. The Zip-Mandelbrot equation is valid for long-time kinetics whereas the Collins-Kimball rate (3) describes the kinetics for all time ranges.
Here, we have studied the kinetics of bimolecular reaction in the presence of immobile obstacles with MLM. When the total volume occupied by obstacles is much smaller than the percolation threshold and the observation time scale is longer than the anomalous to normal diffusion crossover time, the kinetics is still reproducible with the Collins-Kimball theory and Eq. (37). However, it deviates from the theory when the volume occupancy nears or crosses the percolation threshold, wherein anomalous diffusion dominates and the diffusion coefficient approaches zero at the long-time limit. Therefore, to better describe the non-classical kinetics analytically, we should incorporate the anomalous diffusion induced by fractal medium into the theory either phenomenologically [43,49,83] or by extending the Smoluchowski and Collins-Kimball framework using a generalized diffusion equation [84,85].

IV CONCLUSIONS
In contrast to macroscopic and mesoscopic approaches, particle-based methods have the advantage to directly link microscopic parameters to the observed RD behavior, thus providing insights about the underlying mechanisms of the system. MLM shares this same advantage, but with reduced computational costs owing to its fixed step lengths and voxel-based collision detection algorithm. The reduction in computational costs allows MLM to not only simulate non-dilute and crowded intracellular conditions [86] but also track individual molecules on large eukaryotic cells [87] and simulate membrane protein clustering in whole red blood cells [88].
Recently, Grima and colleagues developed a method called vRDME that incorporates volume exclusion into RDME [71,75]. The method can approximate the continuum model very well by matching the steady-state rate constants of both models. We note that vRDME is a type of MLM since each voxel can occupy a molecule and bimolecular reactions occur by colliding reactants. In contrast to vRDME, our work here employs random walk theory and particle-pair formalism to describe bimolecular reactions. Notably, both the effective rate constant and the total rebinding probability on lattice are matched to the corresponding continuum expressions to determine the correct reaction acceptance probability and voxel size. Contrary to the original assumption of Spatiocyte [47], the voxel should be larger than the molecule size (by about 2% for HCP lattice) to be quantitatively accurate. Numerical simulations showed that both the effective rate constant and the asymptotic time-dependent behavior have good agreements with the Collins-Kimball theory in activation-and diffusion-limited cases. MLM also displayed very good consistencies with eGFRD when simulating actual biochemical systems such as protein production-degradation and the dual phosphorylation cycle. Although MLM is analyzed based on the HCP lattice in this work, the theoretical framework is also applicable for other lattice arrangements such as cubic lattice, by simply updating the lattice density and the return probability F (1) (see Appendix B and C).
Despite achieving the same total rebinding probability as the Collins-Kimball theory, the time-dependent behavior of MLM at time scales shorter than t d is different than that theory (Figure 2). One potential solution to obtaining the same behavior at such fine time scales is to make the voxel size smaller than the molecule, similar to the SVTA approach [48]. This would reduce t d but increase the cost of computation significantly because of the finer time steps and the higher number of collision checks required.
MLM captures the effects of excluded volume naturally but comparing on-lattice behavior with continuum is not straightforward since the influence of volume exclusion and the resulting reaction kinetics vary according to the lattice arrangement [44,89]. Moreover, since all diffusing species in this work have the same molecule size, it is not possible to replicate the effects of relative size of interacting molecules. To minimize such lattice artifacts and to better approximate off-lattice volume exclusion, we can improve the size representation of each molecule on lattice by occupying multiple voxels as in the SVTA approach or by employing a hybridized on-and off-lattice approach. Higher spatial resolution of molecules would generate more realistic diffusion behavior in a crowded environment. Alternatively, we can introduce a density-dependent hopping rate as adopted by two previous RDME methods [90,91].
Realistic simulation of intracellular reaction-diffusion processes should also incorporate the influence of inter-molecular potentials such as van der Waals and hydrodynamic forces. By employing contact interactions on lattice as proposed by Fernando et al. [92] or the SVTA approach with interaction potentials [48], it may be possible to incorporate the above forces in MLM. The theoretical framework presented in this work serves as a building block for further development and integration of MLM-based algorithms.

VI ACKNOWLEDGEMENTS
We thank Kozo Nishida for technical advice and support and Kylius Wilkins for critical reading of the manuscript. We thank Ramon Grima for providing the Matlab source of fast Brownian dynamics (point particle version) and Steven Andrews for his help with Smoldyn models used in the performance benchmarks. W.-X.C. acknowledges RIKEN for supporting his doctoral research as an International Program Associate.

APPENDIX B LATTICE INITIAL RATE
Here we provide the derivation of the lattice initial rate, which was done previously in [47]. Given two reacting species A and B, in which A are stationary and B are diffusing. The initial rate constant at time step t , can be estimated using the rate equation as where N i denotes the number of molecules of species i, ∆N c denotes the change in N c and V is the compartment volume. The number of successful reactions in a single step t can be crudely estimated as ∆N C = ZP a where Z = N B N A /N v is the average number of encounter, N v = √ 2V /l 3 is the total number of voxels in a compartment volume V , and P a = P a α is the actual reaction acceptance probability during the encounter.
For the activation-limited scheme, where t = t d and P a = P a , The initial reaction rate is then given by Note that D is the sum of diffusion coefficients of the reacting pair, D A + D B . Similary, for the diffusion-influenced scheme, where t = t d α and P a = P a α, we have

(B.3)
Also note that the physical dimension of k a satisfies cm 3 s −1 .
The above derivation for HCP lattice can be generalized to other lattice arrangements: where d is the packing density of the lattice (e.g. d = π/6 for the simple cubic lattice).

APPENDIX C VOXEL SIZE
As shown in main text, in order to match the MLM with the continuum-based model, the voxel size of HCP lattice has to be chosen such that l = 4πR where R is the molecule size and F (1) ≈ 0.256318 (p. 153 in [62]) is the total return probability on HCP lattice. More generally, the voxel length of any regular lattice arrangement follows that For example, for the simple cubic lattice we have the voxel length: about 8% larger than the molecule size (F (1) for the simple cubic lattice is given in p. 153 of [62]).

APPENDIX D FIRST-PASSAGE TIME DISTRIBUTION ON HCP LATTICE
For n ∈ N, we define P n (s a |s b ) as the voxel occupation probability from s b to s a , that is, the probability of being at voxel s a after n steps, given that the walk started at voxel s b ; F n (s a |s b ) as the first-passage time distribution from s b to s a , that is the probability of arriving at s b for the first time on the nth step, given that the walk started at site s a ; s 0 as the origin voxel, s 1 as the element of the set of immediate neighboring voxels of s 0 , and s 2 as the element of the set of the second nearest neighbor voxels of s 0 . The probability generating function of F n (s 0 |s 0 ) and P n (s 0 |s 0 ) is related through (Eq. I.18 in [60]) where P (s 0 |s 0 ; z) = ∞ n=0 P n (s 0 |s 0 )z n is the lattice Green's function for the face-centered cubic (FCC) lattice as defined in Eqs.(2.6)-(2.9) of [59]: wherein K is the complete elliptic integral of the first kind. For the convenience of calculation, the voxel occupation probability is given as [59] P n (s 0 |s 0 ) = 1 12 n n j=0 n j The first-passage time distribution is related to the voxel occupation probability recursively via: A Activation-limited case (k a k D , α = 1) For P a = 1, the rebinding-time probability distribution F n (s 0 |s 1 ) is equivalent to the firstpassage time distribution F n+1 (s 0 |s 0 ) as mentioned in the main text. Whereas for P a < 1, the rebinding-time probability distribution is given by wherein F j n (s 0 |s 0 ) is the probability of reaching the origin for the jth time at nth step (I.1.9 in [60]): with F 1 n (s 0 |s 0 ) = F n (s 0 |s 0 ). With Eq. (D.10) we can obtain H n (s 0 |s 1 ) recursively via H n (s 0 |s 1 ) = P a n j=1 F j n+j (s 0 |s 0 )(1 − P a ) j−1 , for j ∈ Z + , n ∈ N. (D.11) The generating function of H n (s 0 |s 1 ) is related to the generating function of F n (s 0 |s 0 ): where in the last step we have j−1 k=1 F j k (s 0 |s 0 ) z n = 0 since for all k such that k < j − 1, the return probability is zero. Using (Eq. I.20 in [60]): ∞ n=0 F j n (s 0 |s 0 ) z n = F (s 0 |s 0 ; z) j , (D. 13) in Eq. (D.12) we then have: (D.14) Finally the total rebinding probability of an in-contact pair on lattice is obtained by taking the limit z → 1: where F (1) = F (s 0 |s 0 ; z = 1) ≈ 0.256318 (p. 153 in [62]) is the return probability on HCP lattice.

Rebinding probability at long times
The asymptotic behavior of the rebinding-time probability distribution H n (s 0 |s 1 ) at large n can be estimated directly from the generating function. First we expand the generating function of the return probability P n (s|s) for the HCP lattice around z = 1 up to the O(1 − z) term (see Eq. D.8b in [60] and Eq. A.237 in [62]) where P (1) = P (s|s; z = 1) ≈ 1.344661 and c 1 = 3 3/2 /2π. The corresponding expansion of the generating function of F n (s|s) is then where we have ignored the term equal to or higher than O(1 − z).
Recall that the generating function of the rebinding-time probability distribution for the activation- (D.21) By means of singularity analysis of the generating function (see Eq. 2.3 of [94]), the corresponding asymptotic behavior of H n (s 0 |s 1 ) as n → ∞ is therefore 2 Rate coefficient at long times From the definition of rate coefficient on lattice using the particle-pair formalism, we have the m-step reaction rate coefficient: The first summation term is the total rebinding probability while the second term can be evaluated using the Euler-Maclaurin formula: where we have used the definition nl 2 = 6Dt in the last step. Now we have the asymptotic reaction rate as Using the definition k a (1 − H reb ) = k ef f , and applying the expressions for reaction acceptance probability in Eq. (B.2) and voxel size in Eq. (C.1), we obtain the long-time approximation as which has the exact same form as the continuum case.
The derivation of the effective rate coefficient in the diffusion-influenced case differs from the activation-limited case due to the difference in the simulation scheme (see Algorithm 1 in main text), namely in the presence of non-unity step acceptance probability P w = α. The diffusion step n is therefore no longer the same as the simulation step. Specifically, a successful arrival at a new target voxel (or a successful reaction attempt with a reactant) after n = 1 step could have had multiple k simulation steps in the past with hopping failures (or failed reaction attempts). As a result, the actual simulation time corresponding to n steps is not a single value nt = nt d α, but follows some distribution. The purpose of this section is to derive the long-time asymptotic behavior of the rate coefficient, which is independent of the transient time-dependent behavior. Hence, we parameterize the rebinding-time according to the eventful step n (which will be incremented after a physical movement or a reaction attempt), rather than the actual simulation step k. The time-dependent behavior of rate coefficient on the other hand, will be treated in C.
As shown in the main text, the rebinding-time probability distribution G n (s 0 |s 1 ) is defined as G n+1 (s 0 |s 1 ) = S n (s 1 |s 1 ) p(s 1 → s 0 ), for n ∈ N (D.29) where p(s 1 → s 0 ) is the reaction probability and S n (s 1 |s 1 ) is the in-contact probability of a reactive pair after n steps. The reaction probability is defined as where the nominator term accounts for the probability of hopping to s 0 from s 1 and successfully reacting with the reactant located at s 0 in one diffusion step, while the denominator term comes from the infinite sum representing the total probability of unsuccessful escape to s ∈ {adjacent voxel of s 1 } \ s 0 at the previous simulation step 1 . When P w = 1 and α = 1 as in the activation-limited case, then the reaction probability becomes p(s 1 → s 0 ) = P a P 1 (s 0 |s 1 ).
Next, we derive the generating functions of two first-passage time distributions F n (s 1 |s 1 ) and F n (s 1 |s 2 ) that correspond to the current scheme. We start from F n+1 (s 1 |s 1 ) = s P 1 (s|s 1 )F n (s 1 |s), for n ∈ N = P 1 (s 0 |s 1 )δ n,1 + P 1 (s 1 |s 1 )δ n,0 + P 1 (s 2 |s 1 )F n (s 1 |s 2 ), (D. 31) where the first term on the right-hand side relates to the failed reaction attempt s 1 → s 0 → s 1 2 , the second term describes the hop from s 1 → s 1 , and the last term accounts for the trajectory s 1 → s 2 , which is continued by a series of n steps that have ended up in s 1 again.

C Continuous time limit of the diffusion-influenced scheme
In the diffusion-influenced scenario, Spatiocyte uses a different approach for hopping and reaction. Simulation progresses with a smaller time step t = t d α to resolve fast reaction events. We show that as α becomes smaller, the reaction and hopping events occur in a probabilistic manner that follows exponential time distribution. This property provides us with an approximation to study the time-dependent behavior of the reaction kinetics.

Hopping time distribution
Consider a single particle hopping on a completely vacant lattice. Let P w be the step acceptance probability for a particle heading to a vacant voxel. Then the probability of successful hopping after m trials is The survival probability (no hopping) until mth trial is then (D.62) If we perform the trial every δ sec such that P w = β 1 δ, where β 1 = t −1 d is the average hopping rate per second. The survival probability becomes where t = mδ. Similarly, we have Since P w = α, when α is small enough, the hopping time distribution of a particle approximates the exponential distribution ψ h (t) = exp(−β 1 t), (D. 66) with β 1 = t −1 d .

Reaction time distribution
Consider a reaction pair at an in-contact situation. The survival probability that they are still at the in-contact situation after n steps is S n = (1 − P r − P e ) n , for n ∈ N, (D. 67) where P r = P a α/12 is the reaction probability and P e = 11P w /12 = 11α/12 is the escape probability. Let the simulation trial performed at infinitesimal time δ, such that t = nδ = t d α. where β = (P a + 11)/12t d . Note that the survival probability in this form includes both the probability of reaction and hopping events. Since the two events are independent of each other, the survival probability can be split into two separate terms: where β 1 = P a /12t d is the average reaction rate and β 2 = 1/t d is the average hopping rate. Therefore, the survival probability of the reaction also follows the exponential function where β = β 1 + β 2 . As a consequence, the survival probability of a reactive pair at short time t after step n follows the Poisson distribution: S n (t) = (βt) n n! exp(−βt), for n ∈ N, (D. 73) where S 0 (t) = exp(−βt).

Rate coefficient at long times
Here, we study the time-dependent kinetics of the diffusion-influenced scheme. We start with the definition of continuous rebinding-time probability density, and use it to express the timedependent rate coefficient. Denoting the continuous rebinding-time probability density after (n + 1) steps as g n+1 (t) = β 1 S n (s 1 |s 1 ; t), for n ∈ N, (D. 74) where S n (s 1 |s 1 ; t) = δ n,0 S 0 (s 1 |s 1 ; t) +ˆt 0 dt n j=0 S n−j (s 1 |s 1 ; t − t )F j (s 1 |s 1 ; t), (D. 75) is the survival time probability density of a particle that started and ended at s 1 on the nth step. The first term on the right-hand side of Eq. (D.75) is the initial probability density S 0 (s 1 |s 1 ; t) = exp(−βt), while the last term involves two convolutions: the continuous time convolution and the discrete step convolution nested inside the time convolution. F n (s 1 |s 1 ; t) is the first-passage time density at the nth step, defined as F n (s 1 |s 1 ; t) = F n (s 1 |s 1 ) δ n,1 β 2 exp(−βt) (D. 76) for n ∈ Z + . Intuitively, the first term describes the first-passage time distribution for single step while the second term accounts for the convolution of the probability of time required for the n − 1 steps after the first step. The continuous rebinding-time probability density is related to the rate coefficient of the particle-pair formalism through: F n (s 1 |s 1 ) β 2 s + β 2 n = β 2 s + β F 1 (s 1 |s 1 ) s s + β 2 + ∞ n=0 F n (s 1 |s 1 ) β 2 s + β 2 n = β 2 s + β F 1 (s 1 |s 1 ) s s + β 2 + F s 1 |s 1 ; z = β 2 s + β 2 , where F (s 1 |s 1 ; z) is the generating function, ∞ n=0 F (s 1 |s 1 )z n is as defined in (D.56). Hence, we have where [B] is the concentration of B and k rad (t) is the irreversible rate coefficient according to the radiation boundary condition. Since A is removed from the system via the bimolecular reaction, the concentration of A will eventually reach a steady-state. The corresponding mean lifetime of the decay τ is used to define the steady-state rate coefficient k ss [52]: where k on = 4πDR rad ef f is the macroscopic rate constant, R rad ef f = k a R/(k a +4πRD) is the effective radius and k a is the intrinsic reaction rate constant. The equilibrium concentration of A is then . (E.4)