# Reaction-diffusion kinetics on lattice at the microscopic scale

Wei-Xiang Chew1,2 1 Laboratory for Biologically Inspired Computing, RIKEN
Center for Biosystems Dynamics Research, Suita, Osaka, Japan
2 Department of Physics, Faculty of Science, University of Malaya, 50603,
Kuala Lumpur, Malaysia.
∗E-mail<EMAIL_ADDRESS>Kazunari Kaizu1 1 Laboratory for Biologically
Inspired Computing, RIKEN Center for Biosystems Dynamics Research, Suita,
Osaka, Japan
2 Department of Physics, Faculty of Science, University of Malaya, 50603,
Kuala Lumpur, Malaysia.
∗E-mail<EMAIL_ADDRESS>Masaki Watabe1 1 Laboratory for Biologically Inspired
Computing, RIKEN Center for Biosystems Dynamics Research, Suita, Osaka, Japan
2 Department of Physics, Faculty of Science, University of Malaya, 50603,
Kuala Lumpur, Malaysia.
∗E-mail<EMAIL_ADDRESS>Sithi V. Muniandy2 1 Laboratory for Biologically
Inspired Computing, RIKEN Center for Biosystems Dynamics Research, Suita,
Osaka, Japan
2 Department of Physics, Faculty of Science, University of Malaya, 50603,
Kuala Lumpur, Malaysia.
∗E-mail<EMAIL_ADDRESS>Koichi Takahashi1 and Satya N. V. Arjunan1∗ 1
Laboratory for Biologically Inspired Computing, RIKEN Center for Biosystems
Dynamics Research, Suita, Osaka, Japan
2 Department of Physics, Faculty of Science, University of Malaya, 50603,
Kuala Lumpur, Malaysia.
∗E-mail<EMAIL_ADDRESS>
## Abstract

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. The
theoretical framework of 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 observed 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-dephosphorylation cycle. In addition,
when a reaction compartment is populated with volume-excluding obstacles, MLM
captures the non-classical reaction kinetics caused by anomalous diffusion of
reacting molecules.

##  missingum@section Introduction

#####

In the intracellular environment, macromolecules can be heterogeneously
distributed in space and react stochastically at low concentrations. The
conventional mass action-based approach is insufficient to describe the
reaction-diffusion (RD) behavior of the macromolecules and thus, it is
necessary to incorporate space and stochasticity into the model [1, 2, 3, 4,
5, 6]. Generally, we can represent space as a continuum (off-lattice) or a
discretized lattice model. In the former, each molecule is represented as a
point or a hard-body sphere that propagates via Brownian motion in continuous
space [7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23] .
Bimolecular reaction is often modeled as a collision-based interaction [7, 11,
12, 14, 19, 21] according to the Smoluchowski model of diffusion-influenced
reactions [24, 25]. In some models, the finite size of molecules is taken into
account in the reaction and therefore, the effects of volume-exclusion by
molecules can be reproduced [11, 12, 14, 15, 13, 16, 18, 17, 21]. Although
continuous space-time models are physically consistent, the cost of
computation becomes significant when simulating non-dilute and crowded
conditions in the cell [26].

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 well-mixed
condition, there is a limit to the size of the subvolume [34, 28, 35], which
in turn imposes a limit to the spatial resolution. Diffusion of molecules
across subvolumes is modeled as a first-order 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 well-mixed limit in RDME models and to bridge
mesoscopic and microscopic scales [37, 34, 35, 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, 42, 43, 44, 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 activation-limited. 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.

##  missingum@section 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 continuum-
based framework

Consider an irreversible bimolecular reaction involving two distinct species:

$A+B\xrightarrow{}C,$ (1)

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)$:

$\frac{d[A](t)}{dt}=\frac{d[B](t)}{dt}=-k_{irr}(t)[A](t)[B](t).$ (2)

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}$:

$k_{irr}(t)=\frac{k_{D}k_{a}}{k_{D}+k_{a}}\left[1+\frac{k_{a}}{k_{D}}\Phi\left(\frac{k_{a}}{k_{D}}\sqrt{\frac{t}{\tau^{\prime}}}\right)\right].$
(3)

Here, $k_{D}=4\pi RD$ is the collision rate,
$\Phi(x)=\exp({x^{2}})\text{erfc}(x)$ and
$\tau^{\prime}=(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)\simeq k_{eff}\left[1+\frac{k_{a}}{k_{a}+k_{D}}\frac{R}{\sqrt{\pi
Dt}}\right],$ (4)

at long-time. $k_{eff}$ is the steady-state or the effective reaction rate
constant given by [53]:

$k_{eff}\coloneqq k_{irr}(t\to\infty)=\frac{k_{a}k_{D}}{k_{a}+k_{D}}.$ (5)

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),$ (6)

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)

$p_{reb}(R,t|R,0)=\left(\frac{k_{a}}{4\pi
R^{3}}\right)\left(\frac{k_{a}}{k_{D}}+1\right)\left(\frac{1}{\sqrt{\pi\tau}}-\exp(\tau)\text{erfc}(\sqrt{\tau})\right),$
(7)

with $\tau=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

$k_{irr}(t)=k_{a}\left[1-\int_{0}^{t}p_{reb}(R,\tau|R,0)d\tau\right].$ (8)

At long-time, we then have

$\displaystyle k_{eff}\coloneqq k_{irr}(t\to\infty)$
$\displaystyle=k_{a}\left[1-\int_{0}^{\infty}p_{reb}(R,\tau|R,0)d\tau\right],$
(9)

where the integrated term gives the total rebinding probability:

$P_{reb}=\int_{0}^{\infty}p_{reb}(R,t|R,0)dt=\frac{1}{1+\frac{k_{D}}{k_{a}}}.$
(10)

Therefore, the effective rate constant (5) can also be written in terms of the
total rebinding probability:

$k_{eff}=k_{a}(1-P_{reb}).$ (11)

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$. $t_{d}=l^{2}/6D_{x}$, where $D_{x}$ is the
diffusion coefficient of molecule $x$. 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}=\Delta
N_{C}/Z$. Here $\Delta N_{C}=k_{a}N_{A}N_{B}t_{d}/V$ 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]

$P_{a}=\frac{k_{a}}{3\sqrt{2}(D_{A}+D_{B})l}.$ (12)

The above relation is applicable when the reaction is activation-limited
($k_{a}\ll k_{D}$). For diffusion-influenced reactions ($k_{a}\gg k_{D}$), the
collision rate $Z$ is reduced relative to the production rate $\Delta N_{C}$.
The acceptance probability $P_{a}=\Delta N_{C}/Z$ would then have the issue of
exceeding unity when $\Delta N_{C}>Z$. The Spatiocyte scheme overcomes this
issue by reducing the simulation interval by a factor of $\alpha$ to
$t^{\prime}=t_{d}\alpha$. 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}=\alpha$ and
$P_{a}^{\prime}=P_{a}\alpha$, respectively. Algorithm 1 describes how $\alpha$
is set. In summary, the Spatiocyte scheme operates with $\alpha=1$ when
$P_{a}\leq 1$ (activation-limited case) and with $\alpha<1$ when $P_{a}>1$
(diffusion-influenced case).

Initialization:

$t_{sim}\leftarrow$0, scheduler $S\leftarrow\\{\\}$

for  _each species $x$_ do

$\rho_{x}=$max$\\{P_{axy}\\}$, where $xy$ denotes the pair of reactive species
$x$ and $y$;

$S\leftarrow t^{\prime}=t_{d}\alpha$, where
$\alpha=\left\\{\begin{array}[]{lr}1/\rho_{x},&\text{for }\rho_{x}>1\\\
1,&\text{for }\rho_{x}\leq 1\end{array}\right.$;

reaction acceptance probability $P^{\prime}_{axy}=P_{axy}\alpha$;

step acceptance probability $P_{wx}=\alpha$;

end for

Main loop:

while  _$S\neq\\{\\}$ and $t_{sim}<t_{end}$_ do

$t_{sim}\leftarrow\tau_{x}=$ next event in $S$;

get species identity $x$;

get current voxel location $s_{0}$;

reschedule next event, $\tau_{x}=\tau_{x}+t^{\prime}$

for  _each molecule of species $x$_ do

choose a random target voxel $s_{1}\in$ $\\{\text{nearest neighbor of
}s_{0}\\}$;

if  _$s_{1}$ is vacant_ then

draw $rand$;

if  _$rand\leq P_{wx}$_ then walk succeeded, $s_{0}\leftarrow s_{1}$;

else

walk rejected, $s_{0}\leftarrow s_{0}$;

end if

else if  _$s_{1}$ contains reactant species $y$_ then

draw $rand$;

if  _$rand\leq P^{\prime}_{axy}$_ then

reaction $xy$ accepted, $s_{0}\leftarrow s_{1}$

else

reaction failed and walk rejected, $s_{0}\leftarrow s_{0}$

end if

else

walk rejected, $s_{0}\leftarrow s_{0}$;

end if

end for

end while

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 discrete-space version with a step-dependent rate coefficient on
lattice as [39, 58]

$k_{m}=k_{a}^{\prime}\left[1-\sum_{n=0}^{m}H_{n}\right],\ \text{for
}m,n\in\mathbb{N},$ (13)

where $m$ is the simulation step, which is related to the simulation time by
$6D_{x}t^{\prime}=ml^{2}$, $k_{a}^{\prime}$ 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):

$k_{eff}^{\prime}=\lim_{m\to\infty}k_{m}=k_{a}^{\prime}\left[1-\sum_{n=0}^{\infty}H_{n}\right],$
(14)

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}\ll k_{D}$, $\alpha=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 $n$th 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

$p(s_{0}\rightarrow s_{1})F_{n}(s_{0}|s_{1})=F_{n+1}(s_{0}|s_{0}),$ (15)

where $p(s_{0}\rightarrow 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}\to s_{1}\to 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)=\sum_{n=0}^{\infty}F_{n}(s_{0}|s_{0})z^{n}$ [59] (see
Appendix D).

As for $P_{a}\leq 1$, the trajectories that have undergone failed reaction
attempts before step $n$ are included in the rebinding-time probability
distribution:

$H_{n}(s_{0}|s_{1})=P_{a}\sum_{j=1}^{n}F_{n+j}^{j}(s_{0}|s_{0})(1-P_{a})^{j-1},\text{
for }n\in\mathbb{N},\ j\in\mathbb{Z}^{+},$ (16)

where $F^{j}_{n}(s_{0}|s_{0})$ is the probability to reach the origin for the
$j$th time at the $n$th step (I.1.9 in [60]):

$F_{n}^{j}(s_{0}|s_{0})=\sum_{i=1}^{n}F_{n-i}^{j-1}(s_{0}|s_{0})F_{i}(s_{0}|s_{0}),\text{
for }\ j\in\mathbb{Z}^{+},$ (17)

where $F_{n}^{1}(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):

$H(s_{0}|s_{1};z)=\frac{P_{a}F(s_{0}|s_{0};z)}{F(s_{0}|s_{0};z)(P_{a}-1)+z}.$
(18)

By taking the limit $z\to 1$ on $H(s_{0}|s_{1};z)$, we obtain the total
rebinding probability on lattice as

$H_{reb}=\lim_{z\to 1}H(s_{0}|s_{1};z)=\frac{P_{a}}{P_{a}+\frac{1}{F(1)}-1},$
(19)

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)\approx
0.256318$ ([62], p. 153) for HCP lattice.

Finally, if we set the initial rate $k_{a}^{\prime}=3\sqrt{2}lDP_{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

$\displaystyle k_{eff}^{\prime}$
$\displaystyle=3\sqrt{2}Dl\left(\frac{1}{F(1)}-1\right)\frac{P_{a}}{P_{a}+\frac{1}{F(1)}-1}\
.$ (20)

####  2 Diffusion-influenced case ($k_{a}\gg k_{D}$, $\alpha<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}\to s_{0}),\text{ for
}n\in\mathbb{N},$ (21)

where

$p(s_{1}\to s_{0})=\frac{P_{a}\alpha
P_{1}(s_{0}|s_{1})}{1-(1-P_{w})(1-P_{1}(s_{0}|s_{1}))}$ (22)

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}$
($=\frac{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)=\frac{P_{a}\alpha/12}{1-11(1-\alpha)/12}S(s_{1}|s_{1};z),$
(23)

where $S(s_{1}|s_{1};z)$ is the probability generating function of
$S_{n}(s_{1}|s_{1})$ .

Taking the limit $z\to 1$, we get the total rebinding probability as (Appendix
B )

$G_{reb}=\lim_{z\to 1}G(s_{0}|s_{1};z)=\frac{P_{a}}{P_{a}+\frac{1}{F(1)}-1},$
(24)

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 diffusion-influenced case as

$\displaystyle k_{eff}^{\prime}$ $\displaystyle=k_{a}^{\prime}[1-G_{reb}],$
(25)

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}^{\prime}=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

$P_{a}=\frac{k_{a}}{3\sqrt{2}Dl}.$ (26)

Employing the $G_{reb}=P_{reb}$ relation, the voxel size is found to be about
$2\%$ greater than the molecule size:

$l=\frac{4\pi R}{3\sqrt{2}(\frac{1}{F(1)}-1)}\approx 1.0209R.$ (27)

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:

$\displaystyle k_{eff}^{\prime}$
$\displaystyle=\frac{k_{a}^{\prime}k_{D}^{\prime}}{k_{a}^{\prime}+k_{D}^{\prime}}=k_{D}^{\prime}G_{reb},$
(28) $\displaystyle G_{reb}$
$\displaystyle=\frac{1}{1+k_{D}^{\prime}/k_{a}^{\prime}},$

where $k_{D}^{\prime}=3\sqrt{2}lD({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

$P_{a}=(1/F(1)-1)\left[\frac{3\sqrt{2}(k_{a}+k_{D})(1/F(1)-1)}{4\pi
k_{a}}-1\right]^{-1}.$ (29)

However, this is done at the expense of losing accuracy in the total rebinding
probability, since $G_{reb}\neq P_{reb}$.

For the reversible reaction, $A+B\underset{k_{d}}{\stackrel{{\scriptstyle
k_{a}}}{{\rightleftharpoons}}}C$, the local detailed balance on lattice is
achieved by choosing a lattice dissociation rate constant $k_{d}^{\prime}$
from the following equilibrium constant relation:

$K_{eq}=\frac{k_{a}^{\prime}}{k_{d}^{\prime}}=\frac{k_{a}}{k_{d}}.$ (30)

The MLM method can simulate the dissociation reaction as a first-order process
with rate $k_{d}^{\prime}$ 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. The Python
notebooks used to generate the simulation results reported here are available
at https://github.com/wxchew/MLM. The performance benchmark models for all
tested methods are included in Spatiocyte package (http://spatiocyte.org) as
examples.

##  missingum@sectionI 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}\ll
1$), the expected first rebinding probability at $n$th step is obtained using
Eq. (16), whereas in the diffusion-influenced case ($k_{a}/k_{D}\geq 1$), the
probability is calculated from the generating function $G(s_{0}|s_{1};z)$

$G_{n}(s_{0}|s_{1})=\left.\left[\frac{1}{n!}\frac{d^{n}}{dz^{n}}G(s_{0}|s_{1};z)\right]\right|_{z=0}.$
(31)

Table (1) shows the simulated and the expected theoretical values for
$n\in[1,5]$ steps. The simulation results agree well with the expected values,
with discrepancies never exceeding $0.1\%$. 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 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).

Table 1: Theoretical and simulated rebinding-time probabilities on lattice for activation-limited and diffusion-influenced cases. Simulation parameters: $l=0.01\ \mu\mathrm{m}$, volume = $(100\ l)^{3}$ with periodic boundary, runs = $1\times 10^{9}$. |  $H_{n}(s_{0}|s_{1})$, $P_{a}=0.5$ |  $G_{n}(s_{0}|s_{1})$, $P^{\prime}_{a}=2\alpha$, $\alpha=1/2$
---|---|---
n | Theory | Simulation | Error (%) | Theory | Simulation | Error (%)
1 | 0.0416666 | 0.0416586 | 0.019 | 0.1538461 | 0.1538326 | 0.009
2 | 0.0156250 | 0.0156228 | 0.014 | 0.0473373 | 0.0473431 | 0.012
3 | 0.0107784 | 0.0107779 | 0.005 | 0.0313306 | 0.0313317 | 0.004
4 | 0.0074297 | 0.0074274 | 0.031 | 0.0200584 | 0.0200534 | 0.026
5 | 0.0056802 | 0.0056773 | 0.049 | 0.0147588 | 0.0147496 | 0.062
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\ \mu\mathrm{m}$, $D=1\ \mu\mathrm{m^{2}/s}$, volume = $(10000\ l)^{3}$ with periodic boundary, runs = $1\times 10^{9}$. $\alpha=1/P_{a}$ for diffusion-influenced cases. $k_{a}/k_{D}$ | 0.01 | 0.1 | 1 | 10 | 100
---|---|---|---|---|---
Lattice theory | 0.0062657 | 0.05973879 | 0.397486 | 0.874154 | 0.985988
Simulation | 0.0062672 | 0.05973410 | 0.397459 | 0.874126 | 0.986000
Discrepancy ($\%$) | 0.025 | 0.0078 | 0.0068 | 0.0032 | 0.0012

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)\propto
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.

Figure 2:  The rebinding time of a reactive pair that is initially in-contact.
The rebinding time is sampled from simulations with $k_{a}/k_{D}$ = 0.1, 1 and
100. Markers show the simulation results of Spatiocyte while solid lines
depict the analytical results from the continuum-based scheme (7). The
vertical dashed line marks the characteristic diffusion time step, $t_{d}$.
Simulation parameters: $l=0.01\ \mu\mathrm{m}$, volume = $(10000\ l)^{3}$ with
periodic boundary, runs = $10^{4}$, $D_{A}=1\mu\mathrm{m^{2}s^{-1}}$,
$D_{B}=0$, $\alpha=1/P_{a}$ for the diffusion-influenced case.

Note that in the diffusion-influenced case ($k_{a}/k_{D}=1\text{ 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}$.

####  2 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 $[A]$ and
$[B]$, respectively. We recorded the surviving fraction of $A$ molecules at
each time step. Figure 3(a) displays the survival probability of $A$ and the
expected theoretical curve
$S_{A}(t)=\exp[-[B]\int_{0}^{t}k_{irr}(t^{\prime})dt^{\prime}]$ (Eq. 2.35 in
[65]). From the survival probability, we calculated the time-dependent
reaction rate coefficient using (Eq. 2.1 in [65])

$k_{irr}(t)=-\frac{1}{[B]S_{A}(t)}\frac{dS_{A}(t)}{dt}.$ (32)

We adopted the following discretization scheme for the time derivative to get
the discrete rate coefficient:

$k_{n+1}=-\frac{S_{n+2}-S_{n}}{[B]S_{n+1}\ (t_{n+2}-t_{n})},\text{ for
}n\in\mathbb{Z}^{+},$ (33)

where $n$ is the index of the discretized $S_{A}$ and $t$. The boundary cases
are computed as

$k_{1}=-\frac{S_{2}-S_{1}}{[B]S_{1}\ (t_{2}-t_{1})},\
k_{N}=-\frac{S_{N}-S_{N-1}}{[B]S_{N}\ (t_{N}-t_{N-1})},$ (34)

where $N$ denotes the final time step. The reaction rate coefficient obtained
for various $k_{a}/k_{D}$ ratios are shown in Figure 3(b) along with their
corresponding theoretical curves from Eq. (4).

(a)

(b)

Figure 3: Survival probability and time-dependent rate coefficient. (a)
Survival probability of $A$ in $A+B\xrightarrow{}B$ with $k_{a}/k_{D}$ = 0.1,
1 and 100. (b) Simulated time-dependent rate coefficients of the reaction and
the corresponding long-time approximation of Collins-Kimball (C-K) theory in
Eq. (4). Simulation parameters: volume = $(3.5\ \mu\mathrm{m})^{3}$ with
periodic boundary, $R=0.01\ \mu\mathrm{m}$, $l=0.01\times 1.0209\
\mu\mathrm{m}$, $D_{A}=0$, $D_{B}=1\ \mu\mathrm{m^{2}s^{-1}}$,
$N_{a}=N_{b}=4000$, duration = 0.05 s, runs = $3\times 10^{4}$,
$\alpha=1/P_{a}$ for the diffusion-influenced case.

Recall that the long-time asymptotic variant of the Collins-Kimball theory (4)
has the form

$k_{irr}(t)\simeq C_{1}\left(1+\frac{C_{2}}{\sqrt{t}}\right),$ (35)

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 non-steady-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}$.

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\times 1.0209\ \mu\mathrm{m}$, $D_{A}=0$, $D_{B}=1\ \mu\mathrm{m^{2}s^{-1}}$, volume = $(350\ l)^{3}$ with periodic boundary, $N_{A}=N_{B}=4000$, duration = $0.05\ \mathrm{s}$, runs = $3\times 10^{4}$, $\alpha=1/P_{a}$ for the diffusion-influenced case. $k_{a}/k_{D}$ | 0.1 | 1 | 100
---|---|---|---
Theoretical $C_{1}$ ($\mu\mathrm{m^{3}s^{-1}}$) | 0.011424 | 0.062832 | 0.124420
Simulation |  0.011423$\pm$0.0012 |  0.062848$\pm$0.0029 |  0.124459$\pm$0.0046
Discrepancy ($\%$) | 0.011 | 0.026 | 0.032
Theoretical $C_{2}$ ($\mathrm{s^{1/2}}$) | 0.00051 | 0.00282 | 0.00559
Simulation |  0.00054$\pm$0.01 |  0.00279$\pm$0.0052 |  0.00563$\pm$0.004
Discrepancy ($\%$) | 5.5 | 1.04 | 0.77
m.s.e. of fit |  3.4$\times 10^{-7}$ |  2.2$\times 10^{-6}$ |  4.2$\times 10^{-6}$

###  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 4(a)). Spatiocyte
achieves comparable or better performance than eGFRD in the typical
concentration range of cytoplasmic macromolecules (0.1 to 10 $\mu$M). For
example at 6 $\mu$M in volume $30\ \mathrm{\mu 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\ \mathrm{\mu m}^{3}$ vs.
$3\ \mathrm{\mu m}^{3}$). The simulation times of Spatiocyte scale almost
linearly with the number of molecules ($T\propto N$), which is not apparent
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 4(b) 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.

(a)

(b)

Figure 4: 3D diffusion performance of particle-based methods. Vertical axis,
$T$ shows the run times to diffuse molecules with diffusion coefficient
$D_{x}$ = 1 $\mathrm{\mu m}^{2}\mathrm{s}^{-1}$ in volume $V$ for 10 s. Bottom
horizontal axis, $N$ represents the number of diffusing molecules, while the
top axis shows the corresponding concentration at $V=30\ \mathrm{\mu m}^{3}$.
(a) Molecules are represented as hard-sphere particles with radius $r=l/2$ =
2.5 nm. (b) Molecules are dimensionless point particles that can overlap one
another. eGFRD does not support point particle diffusion and conversely, fast
BD here can only diffuse point particles. Smoldyn simulation interval is set
to the step interval $t_{d}$ (4.17 $\mu$s) of Spatiocyte and fast BD for
comparison. The eGFRD algorithm uses variable time steps. Each model was
simulated for a predefined run time, $t_{r}$ and the resulting simulated time,
$t_{s}$ was recorded. We calculated $T$, the run time in seconds for 10 s of
simulated time with $T=10t_{r}/t_{s}$. $t_{r}$ was set such that at least
hundreds of simulation steps have been completed. The resulting range of
$t_{r}$ was between 1 hour to several days. Solid lines depict the ideal
scaling for Spatiocyte. Vertical dashed lines indicate the typical
concentration range of proteins in the cytoplasm (0.1 to 10 $\mu$M). All
simulations were executed on the same server with Intel Xeon Platinum 8180 2.5
GHz (max 3.80 GHz) CPU, 768 GB memory and Ubuntu 18.04 LTS operating system.

####  2 Reaction

Recently, Andrews [68] benchmarked the performance of Smoldyn, MCell [69],
eGFRD, SpringSaLaD [18] and ReaDDy [16] particle simulators when running the
well-known Michaelis-Menten enzymatic reaction. Smoldyn required the least
amount of time to complete the benchmark. Running the model on our hardware
(see Figure 5 for specifications) with the same 1 ms simulation interval,
Spatiocyte took 113 s, whereas Smoldyn required 31 s. Since it would take too
long for eGFRD to complete the simulation of the original model [68], we
decreased the number of molecules, diffusion coefficients and reaction rates.
The execution times of Spatiocyte, Smoldyn and eGFRD when running the model
with the new parameters are presented in Figure 5. The simulators generated
almost identical results. Spatiocyte and Smoldyn had similar run times,
whereas eGFRD required about one to two orders of magnitude longer. Although
Spatiocyte is about four times slower than Smoldyn when executing the original
model, both had very similar times with the new parameters. Our results
indicate that the relative performance of Spatiocyte and Smoldyn depends on
the model parameters.

Figure 5: Particle simulation performance of the Michaelis-Menten reaction.
Original benchmark model from [66, 68] was modified with volume ($V$) = 90.9
$\mu\mathrm{m}^{3}$, diffusion coefficient ($D_{x}$) = 1
$\mu\mathrm{m}^{2}\mathrm{s}^{-1}$, $k_{1}$ = 0.01
$\mu\mathrm{m}^{3}\mathrm{s}^{-1}$, $k_{2}$ = $k_{3}$ = 0.1 $\mathrm{s}^{-1}$.
Molecule or voxel radius ($r$), simulation or diffusion step interval ($\Delta
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.

###  C Application Examples

We applied MLM to model two fundamental RD systems of intracellular signaling,
the production-degradation process, previously studied using lattice-based
methods [70, 51, 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.

####  1 Production-degradation process

Consider the production and degradation processes of protein $A$ represented
by a zero-order production coupled with a second-order degradation:

$\emptyset{\stackrel{{\scriptstyle k_{1}}}{{\xrightarrow{}}}}A,\
A+B{\stackrel{{\scriptstyle k_{2}}}{{\xrightarrow{}}}}B.$ (36)

The concentration of $A$ will go through an initial transient state before
settling down at a steady-state 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, 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 6(a), 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.

(a)

(b)

(c)

Figure 6: Production-degradation response of A. (a) Time-series profile of $A$
in Eq. (36) simulated with Spatiocyte (using intrinsic rate, $k_{a}$), eGFRD
and RDME. $D_{A}=D_{B}=0.1\ \mu\mathrm{m^{2}s^{-1}}$, molecule radius
$r_{A}=r_{B}=r\in\\{0.005,0.02,0.05\\}\ \mu$m. Note that $r$ represents half
of the subvolume size in RDME, and the actual molecule radius in Spatiocyte
and eGFRD. For comparison, solid line shows the well-mixed model. (b) Mean
equilibrium concentration of $A$ from Spatiocyte, eGFRD, RDME and RDME with
modified propensity (RDMEm) with $D_{A}=D_{B}=D_{x}\in\\{0.01,0.02,0.1\\}\
\mu\mathrm{m^{2}s^{-1}}$. Solid and dashed lines represent expected results
according to the well-mixed model and the microscopic theory respectively. (c)
Steady-state distribution of $A$ from Spatiocyte and eGFRD with $r=0.05\ \mu$m
and $D_{A}=D_{B}=D_{x}\in\\{0.1,0.02\\}\ \mu\mathrm{m^{2}s^{-1}}$. RDMEm
simulated with $r=1$ and $D=0.02$ is also shown for comparison. The frequency
is normalized such that the sum over the bin is unity. Dotted line represents
the well-mixed model simulated using the Gillespie method. Simulation
parameters: $k_{1}=0.1\ \mu\mathrm{m^{-3}s^{-1}}$, $k_{2}=0.02\
\mu\mathrm{m^{3}s^{-1}}$, $[B]=1\ \mu\mathrm{m^{-3}}$, runs = 700, duration
$>10^{4}$ s to achieve steady-state, volume = $100\ \mu\mathrm{m^{3}}$ with
periodic boundary.

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 6(a) and Table 4, there was no discrepancy when the intrinsic rate
$k_{a}$ was used to compute the reaction-acceptance-probability (26). The
well-known relation between $k_{a}$ and $k_{2}$ is given by Eq. (5), wherein
$k_{2}$ is represented by $k_{eff}$. 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 [37, 34, 38].

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\ \mu\mathrm{m^{-3}}$. Discrepancy (%) from the well-mixed model is shown in parenthesis. Simulation parameters: production rate, $k_{1}=0.1\ \mu\mathrm{m^{-3}s^{-1}}$; degradation rate, $k_{2}=0.02\ \mu\mathrm{m^{3}s^{-1}}$; $[B]=1\ \mu\mathrm{m^{-3}}$; volume $=100\ \mu\mathrm{m^{3}}$; $D_{A}=D_{B}=0.1\ \mu\mathrm{m^{2}s^{-1}}$; runs = 600. Simulation scheme | $l=0.01$ | $l=0.04$ | $l=0.1$
---|---|---|---
Spatiocyte with $k_{2}$ | 9.014 (80.28) | 6.023 (20.46) | 5.393 (7.86)
Spatiocyte with $k_{a}$ | 5.009 (0.18) | 4.984 (0.32) | 4.990 (0.2)
eGFRD | 4.968 (0.64) | 4.975 (0.5) | 4.950 (1)

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 diffusion-influenced. In Figure 6(b), 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 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 6(b) 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 6(b)) 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 [70, 37, 34, 38, 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 6(b)). 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 6(c). 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 reaction
$\emptyset\underset{k_{2}}{\stackrel{{\scriptstyle
k_{1}}}{{\rightleftharpoons}}}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.

####  2 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 7(a). The substrate MAPK (K in Figure 7(a)) 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 $\tau_{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 [80, 14].

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 7(b) displays the steady-state response curves of
Spatiocyte and reference theoretical models. Note that since the reactivation
time $\tau_{rel}$ is equal to or less than the diffusion time step $t_{d}$
(given in Figure 7(b)), 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\ \mu\mathrm{m^{2}s^{-1}}$),
whereas at much slower diffusion ($D_{x}=0.06\ \mu\mathrm{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.

(a)

(b)

Figure 7: Effects of rebinding in dual phosphorylation cycle. (a) Reaction
model showing MAPK (K) is first activated into Kp and then Kpp by MAPKK (KK)
in two phosphorylation steps. Kpp is also deactivated by phosphatase (P) in
two dephosphorylation steps to become K again. Enzymes KK and P become
inactive immediately after reacting with their respective substrates and then
relax back to the active state after some delay $\tau_{rel}$. (b) Fraction of
Kpp in response to MAPKK/phosphatase ratio at steady-state. Circle and square
markers denote simulation result using Spatiocyte with $D_{x}=4\
\mu\mathrm{m^{2}s^{-1}}$ and $D_{x}=0.06\ \mu\mathrm{m^{2}s^{-1}}$,
respectively. Dashed and solid lines represent distributive and processive
mechanism models, respectively. Cross and plus markers show the results from
the original Spatiocyte scheme, wherein the voxel and molecule sizes are
exactly the same. We used a short reactivation time, $\tau_{rel}=1\ \mu$s,
relative to $t_{d}$ (for comparison $t_{d}\approx 1\ \mu$s when $D_{x}=4\
\mu\mathrm{m^{2}s^{-1}}$, $t_{d}\ \approx 70\ \mu$s when $D_{x}=0.06\
\mu\mathrm{m^{2}s^{-1}}$) with the total number of substrates,
$\mathrm{K_{total}}=120$. Hysteresis responses from mean-field distributive
model with five-fold substrate concentration ($\mathrm{K_{total}}=600$) are
indicated by dotted and dash-dotted lines with initial conditions
$\mathrm{[K_{pp}]/[K]_{total}}=1$ and $[\mathrm{K_{pp}]/[K]_{total}}=0$,
respectively. Diamond and triangle markers represent Spatiocyte responses with
five-fold substrate concentration at the indicated diffusion coefficient,
$D_{x}$. Simulation parameters: molecule size $l=0.0025\times 1.0209\ \mu$m,
diffusion coefficient $D_{x}$, $\mathrm{[KK]+[P]}=60$, duration = 200 s,
volume = $1\ \mu\mathrm{m^{3}}$ with periodic boundary.

The parameter ranges examined so far have a stable steady-state as
demonstrated by the response curves in Figure 7(b). When the total
concentration of the substrate is increased five-fold, the mean-field theory
generates hysteresis, shown by the dotted and dash-dotted lines. The dotted
line represents the response when initialized with
$\mathrm{[K_{pp}]/[K]_{total}}=1$, whereas the dash-dotted line has the
initial condition $\mathrm{[K_{pp}]/[K]_{total}}=0$. MLM produced similar
responses when the diffusion is fast $(D_{x}=4)$ (diamond markers in Figure
7(b)). 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 7(b). To be fully
consistent with the continuum-based 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].

####  3 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\xrightarrow{}\emptyset$ in the
presence of uniformly distributed immobile obstacles. $E$ and $S$ have the
radius $5\ \mathrm{nm}$ and diffusion coefficient, $D_{0}=1\
\mu\mathrm{m^{2}s^{-1}}$. 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.
Excluded volume is quantified by the lattice occupancy of the obstacles,
$\phi=N_{o}/N_{v}$, where $N_{o}$ and $N_{v}$ are the numbers of obstacles and
total voxels, respectively. Simulation is carried out in a periodic cubic
compartment with length $L=1\ \mu\mathrm{m}$ for a duration of $1000t_{d}$.
Reactants have dilute concentrations, $[S]=5[E]=0.001N_{v}$ and are placed
randomly at the beginning of simulation.

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 8(a) 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 [82, 81],

$D^{\prime}=D_{0}(1-\phi/\phi_{p}),$ (37)

where $\phi_{p}\approx 0.77$ is the percolation threshold for HCP lattice. We
confirmed that the long-time diffusion coefficients obtained for $\phi$ in
Figure 8(a) (dashed lines) are consistent with $D^{\prime}$ in Eq. (37).

Figure 8(b) shows that the survival probability of $E$ decays slower when the
volume occupancy, $\phi$ 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 ($\phi=0$) in Figure 8(c),
there is a good agreement for the simulated $k(t)$ with the Collins-Kimball
rate coefficient (3). As $\phi$ 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 8(c),
$\phi=0.7$), the kinetics begins to deviate from the Collins-Kimball theory.
The deviation is strongest at $\phi=0.8$, which is beyond the percolation
threshold. Note that at lower volume occupany ($\phi=0.3,0.5$), the anomalous
to normal diffusion crossover time in Figure 8(a) is faster than the
observation time in Figure 8(c). 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 8(a), $\phi=0.7$), the effects of anomalous diffusion is visible in
the kinetics (Figure 8(c), $\phi=0.7$). At above the percolation threshold
($\phi=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}\ll 1$) but non-classical
kinetics is observed when it is diffusion-influenced ($k_{a}/k_{D}\gg 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 [83, 43, 49] or by extending the Smoluchowski and Collins-
Kimball framework using a generalized diffusion equation [84, 85].

(a)

(b)

(c)

Figure 8: Diffusion and bimolecular reaction kinetics in crowded compartment.
(a) Time-dependent diffusion coefficient of tracer molecules in the presence
of immobile obstacles at volume occupancy, $\phi$. The diffusion coefficient
at a time point is determined from the mean-squared displacement of simulated
particle trajectories. Dashed lines denote the diffusion coefficient at long-
time as predicted by $D^{\prime}=D_{0}(1-\phi/\phi_{p})$. (b) Survival
probability of $E$ in $E+S\xrightarrow{}\emptyset$ at $\phi$. (c) The
corresponding time-dependent reaction rates (dashed lines) at $\phi$. Solid
lines represent Collins-Kimball theory with the long-time diffusion
coefficient calculated in (a). Simulation parameters: compartment volume =
$(1\ \mu\mathrm{m})^{3}$ with periodic boundary, $R=0.01\ \mu\mathrm{m}$,
$l=0.01\times 1.0209\ \mu\mathrm{m}$, $D_{E}=D_{S}=1\
\mu\mathrm{m^{2}s^{-1}}$, $k_{a}=10k_{D}$, $[S]=5[E]=0.001N_{v}$, duration =
0.02 s.

##  missingum@section 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.

##  missingum@section Author Contributions

W.-X.C., K.K., K.T. and S.N.V.A. designed research; W.-X.C. performed
research; W.-X.C., K.K., M.W., S.V.M., and S.N.V.A. analyzed data; and
W.-X.C., K.K. and S.N.V.A wrote the manuscript. All authors read and commented
on the manuscript.

##  missingum@section 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 A Rebinding probability distribution

The rebinding probability distribution is defined as (Eq. 3.10 in [56] and Eq.
S27 in [14])

$p_{reb}(R,t;R,0)=k_{a}p(R,t;R,0),$ (A.1)

where $p(r,t;r_{0},0)$ is the Green’s function in the diffusion equation:

$\frac{\partial p(r,t;r_{0},0)}{\partial t}=D\nabla^{2}p(r,t;r_{0},0),$ (A.2)

subjected to initial condition

$p(r,0)=\frac{\delta(r-r_{0})}{4\pi r^{2}},$ (A.3)

and boundary conditions such that

$p(r,t)\to 0\text{ as }r\to\infty.$ (A.4)

and

$4\pi R^{2}D\left.\frac{\partial p(r,t;r_{0},0)}{\partial
r}\right|_{r=R}=k_{a}p(R,t;r_{0},0).$ (A.5)

The latter condition is known as the radiation boundary condition. The Green’s
function $p(r,t;r_{0},0)$ has been solved in (p. 368 in [93]) to be

$\displaystyle p(r,t;r_{0},0)$ $\displaystyle=\frac{1}{8\pi
rr_{0}}\frac{1}{\sqrt{\pi
Dt}}\left[\exp{[-(r-r_{0})^{2}/4Dt]}+\exp{[-(r+r_{0}-2R)^{2}/4Dt]}\right.$
(A.6) $\displaystyle\left.-2B\sqrt{\pi
Dt}\exp{\left[B^{2}Dt+B(r+r_{0}-2R)\right]}\text{erfc}\left(B\sqrt{Dt}\right)\right],$

where $B=(1+k_{a}/k_{D})/R$.
For $r=r_{0}=R$, we thus have

$\displaystyle p(R,t;R,0)$ $\displaystyle=\frac{1}{4\pi
R^{2}}\frac{1}{\sqrt{\pi Dt}}\left[1-B\sqrt{\pi
Dt}\exp{\left(B^{2}Dt\right)}\text{erfc}\left(B\sqrt{Dt}\right)\right].$ (A.7)

Finally by substituting Eq. (A.7) into Eq. (A.1), we obtain the probability
distribution

$p_{reb}(R,t;R,0)=\left(\frac{k_{a}}{4\pi
R^{3}}\right)\left(\frac{k_{a}}{k_{D}}+1\right)\left(\frac{1}{\sqrt{\pi\tau}}-\exp(\tau)\text{erfc}(\sqrt{\tau})\right),$
(A.8)

where $\tau=tD(1+k_{a}/k_{D})^{2}/{R^{2}}$.

##  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^{\prime}$, can be estimated using the rate equation as

$k_{a}^{\prime}=\frac{\Delta N_{c}V}{N_{A}N_{B}t^{\prime}},$ (B.1)

where $N_{i}$ denotes the number of molecules of species $i$, $\Delta N_{c}$
denotes the change in $N_{c}$ and $V$ is the compartment volume. The number of
successful reactions in a single step $t^{\prime}$ can be crudely estimated as
$\Delta N_{C}=ZP_{a}^{\prime}$ where $Z=N_{B}N_{A}/N_{v}$ is the average
number of encounter, $N_{v}=\sqrt{2}V/l^{3}$ is the total number of voxels in
a compartment volume $V$, and $P_{a}^{\prime}=P_{a}\alpha$ is the actual
reaction acceptance probability during the encounter.

For the activation-limited scheme, where $t^{\prime}=t_{d}$ and
$P_{a}^{\prime}=P_{a}$, The initial reaction rate is then given by

$\displaystyle k_{a}^{\prime}$
$\displaystyle=\frac{P_{a}^{\prime}l^{3}}{\sqrt{2}t^{\prime}}$ (B.2)
$\displaystyle=\frac{P_{a}l^{3}}{\sqrt{2}t_{d}}$
$\displaystyle=3\sqrt{2}P_{a}Dl,$

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^{\prime}=t_{d}\alpha$
and $P_{a}^{\prime}=P_{a}\alpha$, we have

$\displaystyle k_{a}^{\prime}$ $\displaystyle=\frac{\alpha
P_{a}l^{3}}{\sqrt{2}\alpha t_{d}}$ (B.3) $\displaystyle=3\sqrt{2}P_{a}Dl.$

Also note that the physical dimension of $k_{a}^{\prime}$ satisfies
$cm^{3}s^{-1}$.
The above derivation for HCP lattice can be generalized to other lattice
arrangements:

$k_{a}^{\prime}=\frac{\pi P_{a}Dl}{d},$ (B.4)

where $d$ is the packing density of the lattice (e.g. $d=\pi/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=\frac{4\pi R}{3\sqrt{2}\left(\frac{1}{F(1)}-1\right)},$ (C.1)

where $R$ is the molecule size and $F(1)\approx 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

$l=\frac{4d}{\left.\frac{1}{F(1)}-1\right.}R.$ (C.2)

For example, for the simple cubic lattice we have the voxel length:

$l=\frac{4\pi/6}{\frac{1}{0.340537}-1}R=1.081515R,$ (C.3)

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\in\mathbb{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 $n$th 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])

$F(s_{0}|s_{0};z)=\sum_{n=0}^{\infty}F_{n}(s_{0}|s_{0})z^{n}=1-\frac{1}{P(s_{0}|s_{0};z)},$
(D.1)

where $P(s_{0}|s_{0};z)=\sum_{n=0}^{\infty}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]:

$P(s_{0}|s_{0};z)=\left[\frac{2(1+3\xi^{2})}{\pi(1-\xi)(1+3\xi)}\right]^{2}K(k_{+})K(k_{-}),$
(D.2) $k_{+}^{2}=\frac{16\xi}{(1-\xi)(1+3\xi)^{3}},$ (D.3)
$k_{-}^{2}=\frac{16\xi^{3}}{(1-\xi)^{3}(1+3\xi)},$ (D.4)
$\xi=\frac{-1+\sqrt{1+z/3}}{1+\sqrt{1-z}},$ (D.5)

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})=\frac{1}{12^{n}}\sum_{j=0}^{n}\binom{n}{j}(-4)^{n-j}b_{j},\text{
for }\ j\in\mathbb{N}$ (D.6)

where

$b_{j}=\sum_{k=0}^{j}\binom{j}{k}^{2}\binom{2k}{k}\binom{2j-2k}{j-k}.$ (D.7)

The first-passage time distribution is related to the voxel occupation
probability recursively via:

$F_{n}(s_{0}|s_{0})=P_{n}(s_{0}|s_{0})-\sum_{j=1}^{n-1}P_{n-j}(s_{0}|s_{0})F_{j}(s_{0}|s_{0}),\text{
for }\ j\in\mathbb{Z}^{+}.$ (D.8)

###  A Activation-limited case ($k_{a}\ll k_{D}$, $\alpha=1$)

For $P_{a}=1$, the rebinding-time probability distribution
$F_{n}(s_{0}|s_{1})$ is equivalent to the first-passage 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

$H_{n}(s_{0}|s_{1})=P_{a}F^{1}_{n+1}(s_{0}|s_{0})+P_{a}(1-P_{a})F^{2}_{n+2}(s_{0}|s_{0})+P_{a}(1-P_{a})^{2}F^{3}_{n+3}(s_{0}|s_{0})+...,$
(D.9)

wherein $F^{j}_{n}(s_{0}|s_{0})$ is the probability of reaching the origin for
the $j$th time at $n$th step (I.1.9 in [60]):

$F_{n}^{j}(s_{0}|s_{0})=\sum_{i=1}^{n}F_{n-i}^{j-1}(s_{0}|s_{0})F_{i}(s_{0}|s_{0}),\text{
for }\ j\in\mathbb{Z}^{+},$ (D.10)

with $F_{n}^{1}(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}\sum_{j=1}^{n}F_{n+j}^{j}(s_{0}|s_{0})(1-P_{a})^{j-1},\text{
for }\ j\in\mathbb{Z}^{+},\ n\in\mathbb{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})$:

$\displaystyle H(s_{0}|s_{1};z)$
$\displaystyle=\sum_{n=0}^{\infty}H_{n}(s_{0}|s_{1})z^{n}$ (D.12)
$\displaystyle=P_{a}\sum_{n=0}^{\infty}\sum_{j=1}^{\infty}F_{n+j}^{j}(s_{0}|s_{0})\
(1-P_{a})^{j-1}z^{n}$
$\displaystyle=P_{a}\sum_{j=1}^{\infty}(1-P_{a})^{j-1}z^{-j}\sum_{n=0}^{\infty}F_{n+j}^{j}(s_{0}|s_{0})\
z^{n+j}$
$\displaystyle=P_{a}\sum_{j=1}^{\infty}(1-P_{a})^{j-1}z^{-j}\sum_{n=0}^{\infty}F_{n}^{j}(s_{0}|s_{0})\
z^{n}$

where in the last step we have $\sum_{k=1}^{j-1}F_{k}^{j}(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]):

$\sum_{n=0}^{\infty}F_{n}^{j}(s_{0}|s_{0})\ z^{n}=F(s_{0}|s_{0};z)^{j},$
(D.13)

in Eq. (D.12) we then have:

$\displaystyle H(s_{0}|s_{1};z)$
$\displaystyle=P_{a}\sum_{j=1}^{\infty}(1-P_{a})^{j-1}z^{-j}F(s_{0}|s_{0};z)^{j}$
(D.14)
$\displaystyle=\frac{P_{a}F(s_{0}|s_{0};z)}{F(s_{0}|s_{0};z)(P_{a}-1)+z}.$

Finally the total rebinding probability of an in-contact pair on lattice is
obtained by taking the limit $z\to 1$:

$H_{reb}=\lim_{z\to 1}H(s_{0}|s_{1};z)=\frac{P_{a}}{P_{a}+\frac{1}{F(1)}-1},$
(D.15)

where $F(1)=F(s_{0}|s_{0};z=1)\approx 0.256318$ (p. 153 in [62]) is the return
probability on HCP lattice.

####  1 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])

$P(s|s;z)\approx P(1)-c_{1}\sqrt{1-z}+O(1-z),$ (D.16)

where $P(1)=P(s|s;z=1)\approx 1.344661$ and $c_{1}=3^{3/2}/2\pi$.
The corresponding expansion of the generating function of $F_{n}(s|s)$ is then

$\displaystyle F(s|s;z)$ $\displaystyle=1-\frac{1}{P(s|s;z)}$ (D.17)
$\displaystyle\approx 1-\frac{1}{P(1)-c_{1}\sqrt{1-z}}$ $\displaystyle\approx
1-\frac{1}{P(1)}-\frac{c_{1}}{P(1)^{2}}\sqrt{1-z},$

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-limited case:

$\displaystyle H(s_{0}|s_{1};z)$
$\displaystyle=\frac{P_{a}F(s_{0}|s_{0};z)}{z+F(s_{0}|s_{0};z)(P_{a}-1)}$
(D.18)
$\displaystyle=\frac{P_{a}F(s_{0}|s_{0};z)}{z[1-F(s_{0}|s_{0};z)(1-P_{a})/z]}.$

By the expansion of the denominator we have

$\displaystyle H(s_{0}|s_{1};z)$
$\displaystyle=\frac{P_{a}F(s_{0}|s_{0};z)}{z}\left\\{1+\frac{(1-P_{a})F(s_{0}|s_{0};z)}{z}+\left[\frac{(1-P_{a})F(s_{0}|s_{0};z)}{z}\right]^{2}+...\right\\}$
(D.19)
$\displaystyle=\frac{P_{a}}{z}\left\\{F(s_{0}|s_{0};z)+\frac{(1-P_{a})}{z}F(s_{0}|s_{0};z)^{2}+\left[\frac{(1-P_{a})}{z}\right]^{2}F(s_{0}|s_{0};z)^{3}+...\right\\}.$

Substituting Eq. (D.17) into $H(s_{0}|s_{1};z)$ and collecting the leading
terms gives

$H(s_{0}|s_{1};z)\approx w\sqrt{1-z}+O(1-z),$ (D.20)

where

$\displaystyle w$
$\displaystyle=-\frac{c_{1}P_{a}}{zP(1)^{2}}\left\\{1-2\frac{(P_{a}-1)(P(1)-1)}{zP(1)}+3\left[\frac{(P_{a}-1)(P(1)-1)}{zP(1)}\right]^{2}+...\right\\}$
(D.21)
$\displaystyle=-\frac{c_{1}P_{a}}{zP(1)^{2}}\sum_{n=1}^{\infty}n(-1)^{n+1}\left[\frac{(P_{a}-1)(P(1)-1)}{zP(1)}\right]^{n-1}$
$\displaystyle=-\frac{c_{1}P_{a}}{zP(1)^{2}}\left[1+\frac{(P_{a}-1)(P(1)-1)}{zP(1)}\right]^{-2}$
$\displaystyle=-\frac{c_{1}P_{a}}{z\left[1+P_{a}(P(1)-1)+P(1)(z-1)\right]^{2}}$
$\displaystyle=-\frac{c_{1}P_{a}}{\left[1+P_{a}(P(1)-1)\right]^{2}}.$

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\to\infty$ is therefore

$H_{n}(s_{0}|s_{1})\approx-\frac{w}{2\sqrt{\pi}}n^{-3/2}+O(n^{-5/2}).$ (D.22)

####  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:

$\displaystyle k_{m}$
$\displaystyle=k_{a}^{\prime}\left[1-\sum_{n=0}^{m}H_{n}(s_{0}|s_{1})\right],\
\text{for }m,n\in\mathbb{N},$ (D.23)

which can be rewritten as

$\displaystyle k_{m}$
$\displaystyle=k_{a}^{\prime}\left[1-\sum_{n=0}^{\infty}H_{n}(s_{0}|s_{1})+\sum_{n=m}^{\infty}H_{n}(s_{0}|s_{1})\right].$
(D.24)

The first summation term is the total rebinding probability while the second
term can be evaluated using the Euler-Maclaurin formula:

$\displaystyle\sum_{n=m}^{\infty}H_{n}(s_{0}|s_{1})$
$\displaystyle\approx\int_{m}^{\infty}dn\frac{w}{2\sqrt{\pi}}n^{-3/2}$ (D.25)
$\displaystyle\approx\frac{w}{\sqrt{\pi n}}$
$\displaystyle\approx\frac{lw}{\sqrt{6D\pi t}},$

where we have used the definition $nl^{2}=6Dt$ in the last step.
Now we have the asymptotic reaction rate as

$\displaystyle\lim_{t\to\infty}k(t)$ $\displaystyle\approx
k_{a}^{\prime}\left[1-H_{reb}+\frac{lw}{\sqrt{6D\pi t}}\right].$ (D.26)

After rearrangement we have

$\displaystyle\lim_{t\to\infty}k(t)$ $\displaystyle\approx
k_{a}^{\prime}(1-H_{reb})\left[1+\frac{lw}{(1-H_{reb})\sqrt{6D\pi t}}\right]$
(D.27) $\displaystyle\approx
k_{a}^{\prime}(1-H_{reb})\left[1+\frac{c_{1}P_{a}l}{(1+(P(1)-1)P_{a})\sqrt{6D\pi
t}}\right].$

Using the definition $k_{a}^{\prime}(1-H_{reb})=k_{eff}^{\prime}$, 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

$\lim_{t\to\infty}k(t)\approx
k_{eff}^{\prime}\left[1+\frac{k_{a}R}{(k_{a}+k_{D})\sqrt{\pi Dt}}\right],$
(D.28)

which has the exact same form as the continuum case.

###  B Diffusion-influenced case ($k_{a}\gg k_{D}$, $\alpha<1$).

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}=\alpha$. 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^{\prime}=nt_{d}\alpha$, 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}\to s_{0}),\text{ for
}n\in\mathbb{N}$ (D.29)

where $p(s_{1}\to 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

$\displaystyle p(s_{1}\to s_{0})$ $\displaystyle=P_{a}\alpha
P_{1}(s_{0}|s_{1})\sum_{k=0}^{\infty}\left\\{[1-P_{1}(s_{0}|s_{1})](1-P_{w})\right\\}^{k}$
(D.30) $\displaystyle=\frac{P_{a}\alpha
P_{1}(s_{0}|s_{1})}{1-(1-P_{w})(1-P_{1}(s_{0}|s_{1}))},$

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\in\\{\text{adjacent voxel of }s_{1}\\}\setminus s_{0}$ at the previous
simulation step111there are $k$ simulation steps in between each diffusion
step $n$. When $P_{w}=1$ and $\alpha=1$ as in the activation-limited case,
then the reaction probability becomes $p(s_{1}\to
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

$\displaystyle F_{n+1}(s_{1}|s_{1})$
$\displaystyle=\sum_{s}P_{1}(s|s_{1})F_{n}(s_{1}|s),\ \text{for
}n\in\mathbb{N}$ (D.31)
$\displaystyle=P_{1}(s_{0}|s_{1})\delta_{n,1}+P_{1}(s_{1}|s_{1})\delta_{n,0}+P_{1}(s_{2}|s_{1})F_{n}(s_{1}|s_{2}),$

where the first term on the right-hand side relates to the failed reaction
attempt $s_{1}\to s_{0}\to s_{1}$222only $s_{1}\to s_{0}$ is considered as a
diffusion step, whereas the rejection $s_{0}\to s_{1}$ is not, the second term
describes the hop from $s_{1}\to s_{1}$, and the last term accounts for the
trajectory $s_{1}\to s_{2}$, which is continued by a series of $n$ steps that
have ended up in $s_{1}$ again.

From Eq. (D.31), we obtain the generating function of $F_{n}(s_{1}|s_{1})$ as

$F(s_{1}|s_{1};z)=z^{2}P_{1}(s_{0}|s_{1})+zP_{1}(s_{1}|s_{1})+zP_{1}(s_{2}|s_{1})F(s_{1}|s_{2};z).$
(D.32)

Thus we obtain

$F(s_{1}|s_{2};z)=\frac{F(s_{1}|s_{1};z)-z^{2}P_{1}(s_{0}|s_{1})-zP_{1}(s_{1}|s_{1})}{zP_{1}(s_{2}|s_{1})},$
(D.33)

where

$F(s_{1}|s_{1};z)=1-\frac{z^{2}P_{1}(s_{0}|s_{1})}{P(s_{0}|s_{0};z)-1}$ (D.34)

is given in terms of the generating function of $P_{n}(s_{0}|s_{0})$ (the
detailed derivation of Eq. (D.34) is given in 1).

Now, we define the probability that a particle is in-contact after $n$-step
as:

$S_{n}(s_{1}|s_{1})=\gamma_{1}S_{n-1}(s_{1}|s_{1})+\sum^{n-1}_{m=0}\gamma_{2}S_{m}(s_{1}|s_{1})\
\bar{F}_{n-m-1}(s_{1}|s_{2})+\delta_{n,0}S_{0}(s_{1}|s_{1}),\text{for
}n\in\mathbb{N},$ (D.35)

where the first term accounts for the trajectories $s_{1}\to s_{0}\to s_{1}$
and $s_{1}\to s_{1}$, the second term represents the trajectories $s_{1}\to
s_{2}\to s_{1}$ and the last term accounts for the initial condition. In
detail, the coefficient

$\displaystyle\gamma_{1}$
$\displaystyle=\left[(1-P_{a}\alpha)P_{1}(s_{0}|s_{1})+P_{w}P_{1}(s_{1}|s_{1})\right]\sum_{k=0}^{\infty}\left\\{[1-P_{1}(s_{0}|s_{1})](1-P_{w})\right\\}^{k}$
(D.36)
$\displaystyle=\frac{(1-P_{a}\alpha)P_{1}(s_{0}|s_{1})+P_{w}P_{1}(s_{1}|s_{1})}{1-[1-P_{1}(s_{0}|s_{1})](1-P_{w})},$

accounts for the total probability of arrival at $s_{1}$ from a rejected
reaction attempt (first sub-term) or from the adjacent neighbor $s_{1}$
(second sub-term) given that there was no successful escape to
$s\in\\{\text{adjacent voxel of }s_{1}\\}\setminus s_{0}$ at the last
simulation step $k$ before the arrival,
while the coefficient

$\displaystyle\gamma_{2}$
$\displaystyle=P_{w}P_{1}(s_{2}|s_{1})\sum_{k=0}^{\infty}\left\\{[1-P_{1}(s_{0}|s_{1})](1-P_{w})\right\\}^{k}$
(D.37)
$\displaystyle=\frac{P_{w}P_{1}(s_{2}|s_{1})}{1-[1-P_{1}(s_{0}|s_{1})](1-P_{w})},$

accounts for the total probability of arriving at $s_{2}$ from $s_{1}$ given
that there was no successful escape to $s\in\\{\text{adjacent voxel of
}s_{1}\\}\setminus s_{0}$ at the last simulation step $k$ before the arrival,
and finally $\bar{F}_{n}(s_{1}|s_{2})=F_{n}(s_{1}|s_{2})$ denotes the first-
passage time distribution of the scheme with step-acceptance probability
$P_{w}=\alpha$ (proof given in 2).

We then multiply Eq. (D.35) with $z^{n}$:

$S_{n}(s_{1}|s_{1})\ z^{n}=\gamma_{1}zS_{n-1}(s_{1}|s_{1})\
z^{n-1}+\gamma_{2}z\sum^{n-1}_{m=0}S_{m}(s_{1}|s_{1})\
z^{m}F_{n-m-1}(s_{1}|s_{2})z^{n-m-1}+\delta_{n,0}S_{0}(s_{1}|s_{1})\ z^{n},$
(D.38)

and take the sum to infinity to obtain

$S(s_{1}|s_{1};z)=\gamma_{1}zS(s_{1}|s_{1};z)+\gamma_{2}zS(s_{1}|s_{1};z)F(s_{1}|s_{2};z)+S_{0}(s_{1}|s_{1}).$
(D.39)

After collecting the terms, we obtain the generating function of
$S_{n}(s_{1}|s_{1})$:

$S(s_{1}|s_{1};z)=\frac{S_{0}(s_{1}|s_{1})}{1-\gamma_{1}z-\gamma_{2}zF(s_{1}|s_{2};z)}.$
(D.40)

Substituting Eq. (D.30) and Eq. (D.40) into Eq. (D.29) then gives the
rebinding-time probability distribution:

$G_{n+1}(s_{0}|s_{1})=\frac{P_{a}\alpha
P_{1}(s_{0}|s_{1})}{1-[1-P_{1}(s_{0}|s_{1})](1-P_{w})}S_{n}(s_{1}|s_{1}),\text{
for }n\in\mathbb{N},$ (D.41)

with the corresponding probability generating function

$\displaystyle G(s_{0}|s_{1};z)$ $\displaystyle=\frac{P_{a}\alpha
P_{1}(s_{0}|s_{1})}{1-[1-P_{1}(s_{0}|s_{1})](1-P_{w})}S(s_{1}|s_{1};z).$
(D.42)

In the diffusion-influenced scheme of Spatiocyte, we have
$P_{1}(s_{0}|s_{1})=1/12$, $P_{1}(s_{1}|s_{1})=4/12$,
$P_{1}(s_{2}|s_{1})=7/12$ and $P_{w}=\alpha$. Using these parameters we then
have the following quantities:

$\displaystyle\gamma_{1}$
$\displaystyle=\frac{(1-P_{a}\alpha)+4\alpha}{12\left[1-11(1-\alpha)/12\right]},$
(D.43) $\displaystyle\gamma_{2}$
$\displaystyle=\frac{7\alpha}{12\left[1-11(1-\alpha)/12\right]},$ (D.44)
$F(s_{1}|s_{2};z)=\frac{F(s_{1}|s_{1};z)-z^{2}/12-4z/12}{7z/12},$ (D.45)
$\displaystyle F(s_{1}|s_{2};z=1)$
$\displaystyle=\frac{F(s_{1}|s_{1};z=1)-1/12-4/12}{7/12}$ (D.46)
$\displaystyle=\frac{8-1/F(1)}{7},$

where we have used definition Eq. (D.58) in Eq. (D.46).
Using Eq. (D.46), we obtain the limit of Eq. (D.40) as:

$\displaystyle S(s_{1}|s_{1};z=1)$
$\displaystyle=\left[1-\frac{(1-P_{a}\alpha)+4\alpha}{12\left[1-11(1-\alpha)/12\right]}-\frac{7\alpha}{12\left[1-11(1-\alpha)/12\right]}F(s_{1}|s_{2};z=1)\right]^{-1}$
(D.47)
$\displaystyle=\left[1-\frac{(1-P_{a}\alpha)+4\alpha}{12\left[1-11(1-\alpha)/12\right]}-\frac{7\alpha}{12\left[1-11(1-\alpha)/12\right]}\frac{8-1/F(1)}{7}\right]^{-1}$
$\displaystyle=\frac{12\left[1-11(1-\alpha)/12\right]}{P_{a}\alpha-\alpha+\alpha/F(1)}.$

Finally, we substitute Eq. (D.47) into Eq. (D.42) to obtain

$\displaystyle G(s_{0}|s_{1};1)$
$\displaystyle=\frac{P_{a}\alpha/12}{1-11(1-\alpha)/12}\frac{12\left[1-11(1-\alpha)/12\right]}{P_{a}\alpha-\alpha+\alpha/F(1)}$
(D.48) $\displaystyle=\frac{P_{a}}{P_{a}+\frac{1}{F(1)}-1}$

Therefore, we have the total rebinding probability as:

$G_{reb}=G(s_{0}|s_{1};z=1)=\frac{P_{a}}{P_{a}+\frac{1}{F(1)}-1}.$ (D.49)

####  1 Return probability $F_{n}(s_{1}|s_{1})$

We denote $P_{n}(s|s_{0})$ as the voxel occupation transition probability from
$s_{0}$ to $s$. It is related to $F_{n}(s|s_{0})$ via the convolution relation
([62], p. 121)

$\displaystyle
P_{n}(s|s_{0})=\delta_{ss_{0}}\delta_{n,0}+\sum_{j=1}^{n}F_{j}(s|s_{0})P_{n-j}(s|s),\text{
for }n\in\mathbb{N}.$ (D.50)

If a random walker started at $s_{0}$, it must go through $s_{1}$ before
reaching the destination voxel $s$. Then we have

$P_{n}(s|s_{0})=\delta_{ss_{0}}\delta_{n,0}+\delta_{s_{0}s_{1}}\delta_{n,1}P_{1}(s|s_{1})+\sum_{j=1}^{n}F_{j}(s_{1}|s_{0})P_{n-j}(s|s_{1}).$
(D.51)

Note that $P_{n}(s|s_{1})=P_{n+1}(s|s).$ Thus, with $s_{0}=s_{1}$, we have

$P_{n+1}(s|s)=\delta_{n1}P_{2}(s|s)+\sum_{j=1}^{n}F_{j}(s_{1}|s_{1})P_{n-j+1}(s|s).$
(D.52)

Multiplying both sides with $z^{n+1}$ gives

$z^{n+1}P_{n+1}(s|s)=\delta_{n1}z^{n+1}P_{2}(s|s)+\sum_{j=1}^{n}z^{j}F_{j}(s_{1}|s_{1})z^{n-j+1}P_{n-j+1}(s|s).$
(D.53)

Then taking the sum of both sides from $n=0$ to infinity gives

$P(s|s;z)-P_{0}(s|s)=z^{2}P_{2}(s|s)+F(s_{1}|s_{1};z)[P(s|s;z)-P_{0}(s|s)],$
(D.54)

where

$P(s|s;z)=\sum_{n=0}^{\infty}z^{n}P_{n}(s|s),\
F(s_{1}|s_{1};z)=\sum_{n=0}^{\infty}z^{n}F_{n}(s_{1}|s_{1})\text{ and
}F_{0}(s_{1}|s_{1})=0.$ (D.55)

As such, we have

$F(s_{1}|s_{1};z)=1-\frac{z^{2}P_{2}(s|s)}{P(s|s;z)-1}.$ (D.56)

The total return probability to $s_{1}$ from $s_{1}$ is then

$\displaystyle\sum_{n=0}^{\infty}F_{n}(s_{1}|s_{1})$
$\displaystyle=\lim_{z\to{1^{-}}}F(s_{1}|s_{1};z)$ (D.57)
$\displaystyle=1-\frac{P_{2}(s_{0}|s_{0})}{P(s_{0}|s_{0};1^{-})-1}.$

Using definition Eq. (D.1) and $P_{2}(s_{0}|s_{0})=1/12$, finally we have

$\sum_{n=0}^{\infty}F_{n}(s_{1}|s_{1})=1-\frac{1/F(1)-1}{12}.$ (D.58)

####  2 Return probability $\bar{F}_{n}(s_{1}|s_{2})$

If we increment the step count $n$ for every successful step to a new voxel,
then the first-passage time distribution from $s_{2}$ to $s_{1}$ at step $n$
is given by

$\bar{F}_{n}(s_{1}|s_{2})=\sum^{\infty}_{m=n}\binom{m-1}{n-1}P_{w}^{n}(1-P_{w})^{m-n}F_{n}(s_{1}|s_{2}),\text{
for }n\in\mathbb{Z}^{+},$ (D.59)

where $P_{w}=\alpha$,is the step acceptance probability. It can be shown that

$\displaystyle\bar{F}_{n}(s_{1}|s_{2})$
$\displaystyle=P_{w}^{n}F_{n}(s_{1}|s_{2})\sum^{\infty}_{m=n}\binom{m-1}{n-1}(1-P_{w})^{m-n}$
(D.60) $\displaystyle=P_{w}^{n}F_{n}(s_{1}|s_{2})\frac{1}{P_{w}^{n}}$
$\displaystyle=F_{n}(s_{1}|s_{2}).$

###  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^{\prime}=t_{d}\alpha$ to resolve fast reaction events. We show that as
$\alpha$ 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.

####  1 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

$P_{h}(t=m)=P_{w}(1-P_{w})^{m-1},\text{ for }m\in\mathbb{Z}^{+},$ (D.61)

The survival probability (no hopping) until $m$th trial is then

$\displaystyle P_{h}(t>m)$
$\displaystyle=\sum_{m}^{\infty}P_{w}(1-P_{w})^{m-1}$ (D.62)
$\displaystyle=(1-P_{w})^{m-1}.$

If we perform the trial every $\delta$ sec such that $P_{w}=\beta_{1}\delta$,
where $\beta_{1}=t_{d}^{-1}$ is the average hopping rate per second. The
survival probability becomes

$\displaystyle P_{h}(t>m\delta)$ $\displaystyle=P_{h}(t>t^{\prime})$ (D.63)
$\displaystyle=(1-\beta_{1}\delta)^{m-1},$

where $t^{\prime}=m\delta$.
Similarly, we have

$\displaystyle P_{h}(t>t^{\prime})$
$\displaystyle=(1-\beta_{1}\delta)^{\frac{t^{\prime}}{\delta}-1}$ (D.64)
$\displaystyle=\frac{(1-\beta_{1}\delta)^{\frac{t^{\prime}}{\delta}}}{(1-\beta_{1}\delta)}.$

Taking the limit of small $\delta$, we then have

$\displaystyle\lim_{\delta\to 0}P_{h}(t>t^{\prime})$
$\displaystyle=\left[\lim_{\delta\to
0}(1-\beta_{1}\delta)^{1/\delta}\right]^{t^{\prime}}$ (D.65)
$\displaystyle=\exp(-\beta_{1}t^{\prime}).$

Since $P_{w}=\alpha$, when $\alpha$ is small enough, the hopping time
distribution of a particle approximates the exponential distribution

$\psi_{h}(t)=\exp(-\beta_{1}t),$ (D.66)

with $\beta_{1}=t_{d}^{-1}$.

####  2 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},\text{ for }n\in\mathbb{N},$ (D.67)

where $P_{r}=P_{a}\alpha/12$ is the reaction probability and
$P_{e}=11P_{w}/12=11\alpha/12$ is the escape probability. Let the simulation
trial performed at infinitesimal time $\delta$, such that
$t^{\prime}=n\delta=t_{d}\alpha$.
The survival probability as a function of time is then

$\displaystyle S(t^{\prime})$ $\displaystyle=\lim_{\delta\to 0}S_{n}$ (D.68)
$\displaystyle=\lim_{\delta\to
0}\left[1-\frac{\alpha}{12}(P_{a}+11)\right]^{n}$
$\displaystyle=\lim_{\delta\to
0}\left[1-\frac{\delta}{t_{d}}\frac{(P_{a}+11)}{12}\right]^{t^{\prime}/\delta}$
$\displaystyle=\left[\lim_{\delta\to
0}\left(1-\beta\delta\right)^{1/\delta}\right]^{t^{\prime}}$
$\displaystyle=\exp(-\beta t^{\prime}).$

where $\beta=({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:

$S(t^{\prime})=\exp(-\beta_{1}t^{\prime})\exp\left(-\frac{11\beta_{2}t^{\prime}}{12}\right),$
(D.69)

where $\beta_{1}=P_{a}/12t_{d}$ is the average reaction rate and
$\beta_{2}=1/t_{d}$ is the average hopping rate. Therefore, the survival
probability of the reaction also follows the exponential function

$\psi_{r}(t^{\prime})=\exp(-\beta_{1}t^{\prime}).$ (D.70)

####  3 Time dependent survival probability

In summary, the survival probability of the reaction and hopping events are
(from Eq. (D.66) and Eq. (D.70))

$\displaystyle\psi_{r}(t)$ $\displaystyle=\exp(-\beta_{1}t)\ \text{, where}\
\beta_{1}=\frac{P_{a}}{12t_{d}},$ (D.71) $\displaystyle\psi_{h}(t)$
$\displaystyle=\exp(-\beta_{2}t)\ \text{, where}\ \beta_{2}=\frac{1}{t_{d}}.$

Thus, the survival probability after one step is

$\psi(t)=\psi_{r}(t)\psi_{h}(t)=\exp\\{-\beta t\\},$ (D.72)

where $\beta=\beta_{1}+\beta_{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)=\frac{(\beta t)^{n}}{n!}\exp(-\beta t),\text{ for }n\in\mathbb{N},$
(D.73)

where $S_{0}(t)=\exp(-\beta t)$.

####  4 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 time-dependent rate coefficient.
Denoting the continuous rebinding-time probability density after $(n+1)$ steps
as

$g_{n+1}(t)=\beta_{1}S_{n}(s_{1}|s_{1};t),\text{ for }n\in\mathbb{N},$ (D.74)

where

$S_{n}(s_{1}|s_{1};t)=\delta_{n,0}S_{0}(s_{1}|s_{1};t)+\int_{0}^{t}dt^{\prime}\sum_{j=0}^{n}S_{n-j}(s_{1}|s_{1};t-t^{\prime})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 $n$th 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(-\beta
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
$n$th step, defined as

$\displaystyle F_{n}(s_{1}|s_{1};t)$
$\displaystyle=F_{n}(s_{1}|s_{1})\left[\vphantom{\int}\delta_{n,1}\beta_{2}\exp(-\beta
t)\right.$ (D.76)
$\displaystyle\left.+(1-\delta_{n,1})\int_{0}^{t}dt^{\prime}\beta_{2}\exp(-\beta
t^{\prime})\frac{\beta_{2}^{n}(t-t^{\prime})^{n-1}}{(n-1)!}\exp(-\beta_{2}(t-t^{\prime}))\right],$

for $n\in\mathbb{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:

$k(t)=k_{a}^{\prime}\left[1-\int_{0}^{t}dt^{\prime}g(t^{\prime})\right],$
(D.77)

as shown in the main text. We then take the Laplace transform of $k(t)$ which
is easier to work with:

$s\hat{k}(s)=k_{a}^{\prime}[1-\hat{g}(s)].$ (D.78)

Note that $\hat{g}(s)$ is related to the rebinding-time and survival-time
probability densities via:

$\hat{g}(s)=\sum_{n=1}^{\infty}\hat{g}_{n}(s_{1}|s_{1};s)=\beta_{1}\sum_{n=1}^{\infty}\hat{S}_{n}(s_{1}|s_{1};s).$
(D.79)

The corresponding Laplace transform of Eq. (D.75) is given as

$\hat{S}_{n}(s_{1}|s_{1};s)=\frac{\delta_{n,0}}{s+\beta}+\sum_{j=0}^{n}\hat{S}_{n-j}(s_{1}|s_{1};s)\hat{F}_{j}(s_{1}|s_{1};s),$
(D.80)

where

$\hat{F}_{n}(s_{1}|s_{1};s)=F_{n}(s_{1}|s_{1})\frac{\beta_{2}}{s+\beta}\left[\delta_{n,1}+(1-\delta_{n,1})\left(\frac{\beta_{2}}{s+\beta_{2}}\right)^{n}\right].$
(D.81)

The infinite sum of Eqs.(D.80) and (D.81) are

$\displaystyle\sum_{n=0}^{\infty}\hat{S}_{n}(s_{1}|s_{1};s)$
$\displaystyle=\frac{1}{s+\beta}+\sum_{n=0}^{\infty}\sum_{j=0}^{n}\hat{S}_{n-j}(s_{1}|s_{1};s)\hat{F}_{j}(s_{1}|s_{1};s),$
(D.82)
$\displaystyle=\frac{1}{s+\beta}+\sum_{n=0}^{\infty}\hat{S}_{n}(s_{1}|s_{1};s)\sum_{n=0}^{\infty}\hat{F}_{n}(s_{1}|s_{1};s),$
$\displaystyle=\frac{1}{s+\beta}\left[1-\sum_{n=0}^{\infty}\hat{F}_{n}(s_{1}|s_{1};s)\right]^{-1},$
$\displaystyle\sum_{n=0}^{\infty}\hat{F}_{n}(s_{1}|s_{1};s)$
$\displaystyle=\frac{\beta_{2}}{s+\beta}\left[F_{1}(s_{1}|s_{1})+\sum_{n=2}^{\infty}F_{n}(s_{1}|s_{1})\left(\frac{\beta_{2}}{s+\beta_{2}}\right)^{n}\right]$
(D.83)
$\displaystyle=\frac{\beta_{2}}{s+\beta}\left[F_{1}(s_{1}|s_{1})\frac{s}{s+\beta_{2}}+\sum_{n=0}^{\infty}F_{n}(s_{1}|s_{1})\left(\frac{\beta_{2}}{s+\beta_{2}}\right)^{n}\right]$
$\displaystyle=\frac{\beta_{2}}{s+\beta}\left[F_{1}(s_{1}|s_{1})\frac{s}{s+\beta_{2}}+F\left(s_{1}|s_{1};z=\frac{\beta_{2}}{s+\beta_{2}}\right)\right],$

where $F(s_{1}|s_{1};z)$ is the generating function,
$\sum_{n=0}^{\infty}F(s_{1}|s_{1})z^{n}$ is as defined in (D.56).
Hence, we have

$\sum_{n=0}^{\infty}\hat{S}_{n}(s_{1}|s_{1};s)=[s+\beta-
sF_{1}(s_{1}|s_{1})z-\beta_{2}F(s_{1}|s_{1};z)]^{-1},$ (D.84)

where $z=\beta_{2}/(s+\beta_{2}).$

Substituting Eq. (D.84) into Eq. (D.79) and by the final value theorem, we
obtain the long-time behavior of $k(t)$ by taking the limit $s\to 0$ in Eq.
(D.78).
Assuming the asymptotic Laplace form of the rate coefficient on lattice ([52]
Eq. 2.37a):

$s\hat{k}(s)\approx k_{eff}^{\prime}(1+a_{eff}^{\prime}\sqrt{s/D}+...).$
(D.85)

We then set $s=0$ to obtain the effective lattice reaction rate constant:

$\displaystyle k_{eff}^{\prime}$ $\displaystyle=k_{a}^{\prime}[1-\hat{g}(0)].$
(D.86)

Evaluating $\hat{g}(0)$ by referring to Eq. (D.79), we then get

$\displaystyle k_{eff}^{\prime}$
$\displaystyle=\frac{3\sqrt{2}P_{a}Dl}{1+P_{a}/(1/F(1)-1)},$ (D.87)

which is consistent with the result shown in main text.
The second order term of Eq. (D.85) is evaluated by expanding $s\hat{k}(s)$
around $s=0$:

$\displaystyle\lim_{s\to 0}\frac{d}{d\sqrt{s}}s\hat{k}(s)$
$\displaystyle=\lim_{s\to 0}\frac{d}{d\sqrt{s}}(-k_{a}^{\prime}\hat{g}(s))$
(D.88)
$\displaystyle=\frac{-k_{a}^{\prime}\beta_{1}}{\sqrt{\beta_{2}}}\lim_{q\to
0}\frac{d}{dq}\sum_{n=0}^{\infty}\hat{S}_{n}(s_{1}|s_{1};s)$
$\displaystyle=k_{a}^{\prime}\beta_{1}\sqrt{\beta_{2}}\left[\beta-\beta_{2}F(1)\right]^{-2}\lim_{q\to
0}\frac{d}{dq}F(q)$
$\displaystyle=\frac{2k_{a}^{\prime}P_{a}}{\sqrt{\beta_{2}}}\left[1+\frac{P_{a}}{1/F(1)-1}\right]^{-2}\lim_{q\to
0}q\frac{d}{dz}P(s_{0}|s_{0};z),$

where $q=\sqrt{s/\beta_{2}}$ and $z=1/(q^{2}+1)$.
Thus, by comparing the terms we obtain

$\displaystyle a_{eff}^{\prime}$
$\displaystyle=\frac{\sqrt{D}}{k_{eff}^{\prime}}\frac{2k_{a}^{\prime}P_{a}b}{\sqrt{\beta_{2}}}[P_{a}\\{P(s_{0}|s_{0};1)-1\\}+1]^{-2}$
(D.89) $\displaystyle=\frac{\sqrt{2/3}bP_{a}l}{1+P_{a}/(1/F(1)-1)},$

where

$b=\lim_{q\to 0}q\frac{d}{dz}P(s_{0}|s_{0};z)=\frac{3\sqrt{3}}{4\pi}.$ (D.90)

Applying the definitions of reaction acceptance probability (B.2) and voxel
size (C.1), we obtain

$a_{eff}^{\prime}=\frac{k_{a}R}{k_{a}+k_{D}}.$ (D.91)

Note that the corresponding time domain form of Eq. (D.85) is given as

$k(t)\approx k_{eff}^{\prime}\left[1+a^{\prime}_{eff}\sqrt{\pi
Dt}+...\right].$ (D.92)

Hence, the long-time behavior of the lattice rate coefficient follows the same
form as in the continuum case:

$k(t)\approx k_{eff}^{\prime}\left[1+\frac{k_{a}R}{(k_{a}+k_{D})\sqrt{\pi
Dt}}\right].$ (D.93)

##  Appendix E Production-degradation process

In the coupled reactions $\emptyset{\stackrel{{\scriptstyle
k_{1}}}{{\xrightarrow{}}}}A$, $A+B{\stackrel{{\scriptstyle
k_{2}}}{{\xrightarrow{}}}}B$, the survival probability of a newly produced $A$
molecule in an equilibrated pool of $B$ is

$S_{rad}(t|eq)=\exp\left[-[B]\int_{0}^{t}k_{rad}(t^{\prime})dt^{\prime}\right].$
(E.1)

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 $\tau$ is used to define the steady-state rate coefficient $k_{ss}$
[52]:

$\tau=\frac{1}{[B]k_{ss}}\equiv\int_{0}^{\infty}S_{rad}(t|eq)=\hat{S}_{rad}(0|eq),$
(E.2)

where the hat denotes Laplace transform.
For small $[B]$, $k_{ss}$ is given by (Eq. 4.5 in [52])

$k_{ss}\approx
k_{on}\left\\{1+\left[4\pi(R_{eff}^{rad})^{3}[B]\right]^{1/2}\right\\},$ (E.3)

where $k_{on}=4\pi DR_{eff}^{rad}$ is the macroscopic rate constant,
$R_{eff}^{rad}=k_{a}R/(k_{a}+4\pi RD)$ is the effective radius and $k_{a}$ is
the intrinsic reaction rate constant.
The equilibrium concentration of $A$ is then

$[A]=\frac{k_{1}}{k_{ss}[B]}.$ (E.4)

## References

  * [1] Martin Howard and Andrew D Rutenberg.  Pattern formation inside bacteria: Fluctuations due to the low copy number of proteins.  Phys. Rev. Lett., 90(12):128102, 2003.
  * [2] Michael Klann and Heinz Koeppl.  Spatial simulations in systems biology: From molecules to cells.  Int. J. Mol. Sci., 13(6):7798–7827, 2012.
  * [3] Stefan Engblom, Andreas Hellander, and Per Lötstedt.  Multiscale simulation of stochastic reaction-diffusion networks.  In Stochastic Processes, Multiscale Modeling, and Numerical Methods for Computational Cellular Biology, pages 55–79. Springer, 2017.
  * [4] Ximing Li and William R. Holmes.  Biophysical attributes that affect CaMKII activation deduced with a novel spatial stochastic simulation approach.  PLoS Comput. Biol., 14(2):1–32, 2018.
  * [5] Tyler Earnest, John Andrew Cole, and Zaida Luthey-Schulten.  Simulating biological processes: Stochastic physics from whole cells to colonies.  Rep. Prog. Phys., 2018.
  * [6] Stephen Smith and Ramon Grima.  Spatial stochastic intracellular kinetics: A review of modelling approaches.  Bull. Math. Biol., pages 1–50, 2018.
  * [7] Joel R Stiles, Thomas M Bartol Jr, Edwin E Salpeter, and Miriam M Salpeter.  Monte Carlo simulation of neuro-transmitter release using MCell, a general simulator of cellular physiological processes.  In Computational neuroscience, pages 279–284. Springer, 1998.
  * [8] Steven J Plimpton and Alex Slepoy.  Microbial cell modeling via reacting diffusive particles.  J. Phys. : Conf. Ser., 16(1):305, 2005.
  * [9] Jeroen S van Zon and Pieter Rein 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.
  * [10] Chris Sanford, Matthew L K Yip, Carl White, and John Parkinson.  Cell++-simulating biochemical pathways.  Bioinformatics, 22(23):2918–2925, 2006.
  * [11] Douglas Ridgway, Gordon Broderick, Ana Lopez-Campistrous, Melania Ru’aini, Philip Winter, Matthew Hamilton, Pierre Boulanger, Andriy Kovalenko, and Michael J Ellison.  Coarse-grained molecular simulation of diffusion and reaction kinetics in a crowded virtual cytoplasm.  Biophys. J., 94(10):3748–3759, 2008.
  * [12] Michael J Byrne, M Neal Waxham, and Yoshihisa Kubota.  Cellular dynamic simulator: An event driven molecular simulation environment for cellular physiology.  Neuroinformatics, 8(2):63–82, 2010.
  * [13] Gerd Gruenert, Bashar Ibrahim, Thorsten Lenser, Maiko Lohel, Thomas Hinze, and Peter Dittrich.  Rule-based spatial modeling with diffusing, geometrically constrained molecules.  BMC Bioinf., 11(1):307, 2010.
  * [14] Koichi Takahashi, Sorin Tănase-Nicola, and Pieter Rein Ten Wolde.  Spatio-temporal correlations can drastically change the response of a MAPK pathway.  Proc. Natl. Acad. Sci. USA, 107(6):2473–2478, 2010.
  * [15] Michael Klann, Heinz Koeppl, and Matthias Reuss.  Spatial modeling of vesicle transport and the cytoskeleton: The challenge of hitting the right road.  PLoS One, 7(1):e29645, 2012.
  * [16] Johannes Schöneberg and Frank Noé.  ReaDDy-a software for particle-based reaction-diffusion dynamics in crowded cellular environments.  PLoS One, 8(9):e74261, 2013.
  * [17] Steven S Andrews.  Smoldyn: Particle-based simulation with rule-based modeling, improved molecular interaction and a library interface.  Bioinformatics, 33(5):710–717, 2016.
  * [18] Paul J. Michalski and Leslie M. Loew.  SpringSaLaD: A spatial, particle-based biochemical simulation platform with excluded volume.  Biophys. J., 110(3):523–529, 2016.
  * [19] Teresa Lehnert and Marc Thilo Figge.  Dimensionality of motion and binding valency govern receptor-ligand kinetics as revealed by agent-based modeling.  Front. Immunol., 8:1692, 2017.
  * [20] A. T. Bittig and A. M. Uhrmacher.  ML-Space: Hybrid spatial Gillespie and particle simulation of multi-level rule-based models in cell biology.  IEEE/ACM Trans. Comput. Biol. Bioinformatics, 14(6):1339–1349, November 2017.
  * [21] Zahedeh Bashardanesh and Per Lötstedt.  Efficient Green’s function reaction dynamics (GFRD) simulations for diffusion-limited, reversible reactions.  J. Comput. Phys., 2017.
  * [22] Aleksandar Donev, Chiao-Yu Yang, and Changho Kim.  Efficient reactive Brownian dynamics.  J. Chem. Phys., 148(3):034103, 2018.
  * [23] Alireza Sayyidmousavi and Katrin Rohlf.  Reactive multi-particle collision dynamics with reactive boundary conditions.  Phys. Biol., 15(4):046007, 2018.
  * [24] MV Smoluchowski.  Mathematical theory of the kinetics of the coagulation of colloidal solutions.  Z. Phys. Chemie, 92:129–168, 1917.
  * [25] Frank C Collins and George E Kimball.  Diffusion-controlled reaction rates.  J. Colloid Sci., 4(4):425–437, 1949.
  * [26] Kouichi Takahashi, Satya Nanda Vel Arjunan, and Masaru Tomita.  Space in systems biology of signaling pathways–towards intracellular molecular crowding in silico.  FEBS Lett., 579(8):1783–8, March 2005.
  * [27] Florence Baras and M Malek Mansour.  Reaction-diffusion master equation: A comparison with microscopic simulations.  Phys. Rev. E, 54(6):6139, 1996.
  * [28] Johan Hattne, David Fange, and Johan Elf.  Stochastic reaction-diffusion simulation with MesoRD.  Bioinformatics, 21(12):2923–2924, 2005.
  * [29] Brian Drawert, Stefan Engblom, and Andreas Hellander.  URDME: A modular framework for stochastic simulation of reaction-transport processes in complex geometries.  BMC Syst. Biol., 6(1):76, 2012.
  * [30] Ann E Cowan, Ion I Moraru, James C Schaff, Boris M Slepchenko, and Leslie M Loew.  Spatial modeling of cell signaling networks.  Methods Cell Biol., 110:195, 2012.
  * [31] Iain Hepburn, Weiliang Chen, Stefan Wils, and Erik De Schutter.  STEPS: Efficient simulation of stochastic reaction–diffusion models in realistic morphologies.  BMC Syst. Biol., 6(1):36, 2012.
  * [32] Elijah Roberts, John E Stone, and Zaida Luthey-Schulten.  Lattice microbes: High-performance stochastic simulation method for the reaction-diffusion master equation.  J. Comput. Chem., 34(3):245–255, 2013.
  * [33] Brian Drawert, Andreas Hellander, Ben Bales, Debjani Banerjee, Giovanni Bellesia, Bernie J Daigle Jr, Geoffrey Douglas, Mengyuan Gu, Anand Gupta, Stefan Hellander, et al.  Stochastic simulation service: Bridging the gap between the computational expert and the biologist.  PLoS Comput. Biol., 12(12):e1005220, 2016.
  * [34] Stefan Hellander, Andreas Hellander, and Linda Petzold.  Reaction-diffusion master equation in the microscopic limit.  Phys. Rev. E, 85(4):1–5, 2012.
  * [35] Stefan Hellander, Andreas Hellander, and Linda Petzold.  Reaction rates for mesoscopic reaction-diffusion kinetics.  Phys. Rev. E, 91(2):023312, 2015.
  * [36] Daniel T Gillespie.  A general method for numerically simulating the stochastic time evolution of coupled chemical reactions.  J. Comput. Phys., 22(4):403–434, 1976.
  * [37] David Fange, Otto G Berg, Paul Sjöberg, and Johan Elf.  Stochastic reaction-diffusion kinetics in the microscopic limit.  Proc. Natl. Acad. Sci. USA, 107(46):19820–5, 2010.
  * [38] Samuel A. Isaacson.  A convergent reaction-diffusion master equation.  J. Chem. Phys., 139(5), 2013.
  * [39] David C Torney and Harden M McConnell.  Diffusion-limited reaction rate theory for two-dimensional systems.  In Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, volume 387, pages 147–170. The Royal Society, 1983.
  * [40] F Montalenti and R Ferrando.  Probability of dimer reassociation in two dimensions.  Phys. Rev. E, 61(4):3411, 2000.
  * [41] Michael J Saxton.  Chemically limited reactions on a percolation cluster.  J. Chem. Phys., 116(1):203–208, 2002.
  * [42] Hugues Berry.  Monte carlo simulations of enzyme reactions in two dimensions: Fractal kinetics and spatial segregation.  Biophys. J., 83(4):1891–1901, October 2002.
  * [43] S. Schnell and T. E. Turner.  Reaction kinetics in intracellular environments with macromolecular crowding: Simulations and rate laws.  Prog. Biophys. Mol. Biol., 85(2-3):235–260, 2004.
  * [44] R. Grima and S. Schnell.  A systematic investigation of the rate laws valid in intracellular environments.  Biophys. Chem., 124(1):1–10, 2006.
  * [45] Laurier Boulianne, Sevin Al Assaad, Michel Dumontier, and Warren Gross.  GridCell: A stochastic particle-based biological system simulator.  BMC Syst. Biol., 2(1):66, July 2008.
  * [46] Jeremy D Schmit, Ercan Kamber, and Jané Kondev.  Lattice model of diffusion-limited bimolecular chemical reactions in confined environments.  Phys. Rev. Lett., 102(21):218302, 2009.
  * [47] Satya Nanda Vel Arjunan and Masaru Tomita.  A new multicompartmental reaction-diffusion modeling method links transient membrane attachment of E. coli MinE to E-ring formation.  Syst. Synth. Biol., 4(1):35–53, 2010.
  * [48] Daniel T. Gillespie, Effrosyni Seitaridou, and Carol A. Gillespie.  The small-voxel tracking algorithm for simulating chemical reactions among diffusing molecules.  J. Chem. Phys., 141(23), 2014.
  * [49] Laura Pitulice, Eudald Vilaseca, Isabel Pastor, Sergio Madurga, Josep LluÍs Garcés, Adriana Isvoran, and Francesc Mas.  Monte Carlo simulations of enzymatic reactions in crowded media. effect of the enzyme-obstacle relative size.  Math. Biosci., 251:72–82, 2014.
  * [50] David Gomez and Stefan Klumpp.  Biochemical reactions in crowded environments: Revisiting the effects of volume exclusion with simulations.  Front. Phys., 3:45, 2015.
  * [51] Marc Sturrock.  Stochastic reaction-diffusion algorithms for macromolecular crowding.  Phys. Biol., 13(3):036010, 2016.
  * [52] Noam Agmon and Attila Szabo.  Theory of reversible diffusion-influenced reactions.  J. Chem. Phys., 92(9):5270–5284, 1990.
  * [53] Richard M Noyes.  Effects of diffusion rates on chemical kinetics.  Prog. React. Kinet. Mech., 1:129–160, 1961.
  * [54] Richard M. Noyes.  A treatment of chemical kinetics with special applicability to diffusion controlled reactions.  J. Chem. Phys., 22(8):1349, 1954.
  * [55] Otto G Berg.  On diffusion-controlled dissociation.  Chem. Phys., 31(1):47–57, 1978.
  * [56] K Razi Naqvi, K J Mork, and S Waldenstrom.  Diffusion-controlled reaction kinetics. Equivalence of the particle pair approach of Noyes and the concentration gradient approach of Collins and Kimball.  J. Phys. Chem., 84(11):1315–1319, 1980.
  * [57] Steven S Andrews and Dennis Bray.  Stochastic simulation of chemical reactions with spatial resolution and single molecule detail.  Phys. Biol., 1(3):137, 2004.
  * [58] Elliott W. Montroll.  Random walks on lattices. III. calculation of first passage times with application to exciton trapping on photosynthetic units.  J. Math. Phys., 10(4):753–765, 1969.
  * [59] G S Joyce.  On the cubic modular transformation and the cubic lattice Green functions.  J. Phys. A: Math. Gen., 31(22):5105, 1998.
  * [60] Elliott W Montroll and George H Weiss.  Random walks on lattices. II.  J. Math. Phys., 6(2):167–181, 1965.
  * [61] Shunya Ishioka and Masahiro Koiwa.  Random walks on diamond and hexagonal close packed lattices.  Philos. Mag. A, 37(4):517–533, 1978.
  * [62] Barry D Hughes.  Random walks and random environments.  Clarendon Press; Oxford University Press, 1995.
  * [63] Kazunari Kaizu, Kozo Nishida, Yuki Sakamoto, Suguru Kato, Toru Niina, Naoki Nishida, Moriyoshi Koizumi, Naohiro Aota, and Koichi Takahashi.  ecell/ecell4: V4.1.3.  https://github.com/ecell/ecell4, 2018.
  * [64] Wolfgang Pfluegl and Robert J Silbey.  Long-time properties of random walks with a single trap.  Phys. Rev. E, 58(4):4128, 1998.
  * [65] Attila Szabo.  Theory of diffusion-influenced fluorescence quenching.  J. Phys. Chem., 93(19):6929–6939, 1989.
  * [66] Steven S Andrews, Nathan J Addy, Roger Brent, and Adam P Arkin.  Detailed simulations of cell biology with Smoldyn 2.1.  PLoS Comput. Biol., 6(3):e1000705, 2010.
  * [67] Stephen Smith and Ramon Grima.  Fast simulation of Brownian dynamics in a crowded environment.  J. Chem. Phys., 146(2):024105, 2017.
  * [68] Steven S Andrews.  Particle-Based Stochastic Simulators.  In Encyclopedia of Computational Neuroscience. Springer, 2018.
  * [69] Rex A. Kerr, Thomas M. Bartol, Boris Kaminsky, Markus Dittrich, Jen-Chien Jack Chang, Scott B. Baden, Terrence J. Sejnowski, and Joel R. Stiles.  Fast Monte Carlo simulation methods for biological reaction-diffusion systems in solution and on surfaces.  SIAM J. Sci. Comput., 30(6):3126–3149, 2008.
  * [70] Radek Erban and S Jonathan Chapman.  Stochastic modelling of reaction-diffusion processes: Algorithms for bimolecular reactions.  Phys. Biol., 6(4):046001, 2009.
  * [71] Claudia Cianci, Stephen Smith, and Ramon Grima.  Molecular finite-size effects in stochastic models of equilibrium chemical systems.  J. Chem. Phys., 144(8):084101, 2016.
  * [72] Lufen Chang and Michael Karin.  Mammalian MAP kinase signalling cascades.  Nature, 410(6824):37, 2001.
  * [73] James E Ferrell and Ramesh R Bhatt.  Mechanistic studies of the dual phosphorylation of mitogen-activated protein kinase.  J. Biol. Chem., 272(30):19008–19016, 1997.
  * [74] Kazuhiro Aoki, Masashi Yamada, Katsuyuki Kunida, Shuhei Yasuda, and Michiyuki Matsuda.  Processive phosphorylation of ERK MAP kinase in mammalian cells.  Proc. Natl. Acad. Sci. USA, 108(31):12675–12680, 2011.
  * [75] Stephen Smith, Claudia Cianci, and Ramon Grima.  Analytical approximations for spatial stochastic gene expression in single cells and tissues.  J. R. Soc. Interface, 13(118):20151051, 2016.
  * [76] George G Szpiro.  Kepler’s conjecture.  Hoboken, NJ: Wiley, 2003.
  * [77] Steven B Zimmerman and Stefan O Trach.  Estimation of macromolecule concentrations and excluded volume effects for the cytoplasm of Escherichia coli.  J. Mol. Biol., 222(3):599–620, 1991.
  * [78] Satya N V Arjunan and Koichi Takahashi.  Multi-algorithm particle simulations with Spatiocyte.  In: Kihara D. (eds) Protein Function Prediction. Methods in Molecular Biology, vol 1611. Humana Press, New York, NY, pages 219–236, 2017\.
  * [79] Michael A Gibson and Jehoshua Bruck.  Efficient exact stochastic simulation of chemical systems with many species and many channels.  J. Phys. Chem. A, 104(9):1876–1889, 2000.
  * [80] Johan Elf and Måns Ehrenberg.  Spontaneous separation of bi-stable biochemical systems into spatial domains of opposite phases.  Syst. Biol., 1(2):230–236, 2004.
  * [81] Eudald Vilaseca, Adriana Isvoran, Sergio Madurga, Isabel Pastor, Josep LluÍs Garcés, and Francesc Mas.  New insights into diffusion in 3D crowded media by Monte Carlo simulations: Effect of size, mobility and spatial distribution of obstacles.  Phys. Chem. Chem. Phys., 13(16):7396–7407, 2011.
  * [82] Michael J Saxton.  Lateral diffusion in an archipelago. distance dependence of the diffusion coefficient.  Biophys. J., 56(3):615–622, 1989.
  * [83] Raoul Kopelman.  Fractal reaction kinetics.  Science, 241(4873):1620–1626, 1988.
  * [84] AV Barzykin and M Tachiya.  Diffusion-influenced reaction kinetics on fractal structures.  J. Chem. Phys., 99(12):9591–9597, 1993.
  * [85] Jaeyoung Sung, Hwangseo Park, and Sangyoub Lee.  Green’s function approach to nonclassical reaction kinetics in fractal media.  Phys. Rev. Lett., 86(11):2463, 2001.
  * [86] S S Andrews, S N V Arjunan, G Balbo, A T Bittig, J Feret, K Kaizu, and F Liu.  Simulating macromolecular crowding with particle and lattice-based methods (Team 3).  In: Gilbert, D., Heiner, M., Takahashi, K. and Uhrmacher, A. M. (eds). Multiscale Spatial Computational Systems Biology, pages 170–187, 2015\.
  * [87] Masaki Watabe, Satya N V Arjunan, Seiya Fukushima, Kazunari Iwamoto, Jun Kozuka, Satomi Matsuoka, Yuki Shindo, Masahiro Ueda, and Koichi Takahashi.  A computational framework for bioimaging simulation.  PloS One, 10(7):e0130089, 2015.
  * [88] Hanae Shimo, Satya Nanda Vel Arjunan, Hiroaki Machiyama, Taiko Nishino, Makoto Suematsu, Hideaki Fujita, Masaru Tomita, and Koichi Takahashi.  Particle simulation of oxidation induced band 3 clustering in human erythrocytes.  PLoS Comput. Biol., 11(6):e1004210, 2015.
  * [89] Lina Meinecke and Markus Eriksson.  Excluded volume effects in on-and off-lattice reaction-diffusion models.  IET Syst. Biol., 11(2):55–64, 2016.
  * [90] Ramon Grima and Santiago Schnell.  A mesoscopic simulation approach for modeling intracellular reactions.  J. Stat. Phys., 128(1-2):139–164, 2007.
  * [91] Claudia Cianci, Stephen Smith, and Ramon Grima.  Capturing Brownian dynamics with an on-lattice model of hard-sphere diffusion.  Phys. Rev. E, 95(5):052118, 2017.
  * [92] Anthony E Fernando, Kerry A Landman, and Matthew J Simpson.  Nonlinear diffusion and exclusion processes with contact interactions.  Phys. Rev. E, 81(1):011903, 2010.
  * [93] Horatio Scott Carslaw and John Conrad Jaeger.  Conduction of heat in solids.  Oxford: Clarendon Press, 1959, 2nd ed., 1959.
  * [94] Philippe Flajolet and Andrew Odlyzko.  Singularity analysis of generating functions.  SIAM J. Discrete Math., 3(2):216–240, 1990.