Surface reaction-diffusion kinetics on lattice at the microscopic scale

Microscopic models of reaction-diffusion processes on the cell membrane can link local spatiotemporal effects to macroscopic self-organized patterns often observed on the membrane. Simulation schemes based on the microscopic lattice method (MLM) can model these processes at the microscopic scale by tracking individual molecules, represented as hard-spheres, on fine lattice voxels. Although MLM is simple to implement and is generally less computationally demanding than off-lattice approaches, its accuracy and consistency in modeling surface reactions have not been fully verifed. Using the Spatiocyte scheme, we study the accuracy of MLM in diffusion-influenced surface reactions. We derive the lattice-based bimolecular association rates for two-dimensional surface-surface reaction and one-dimensional volume-surface adsorption according to the Smoluchowski-Collins-Kimball model and random walk theory. We match the time-dependent rates on lattice with off-lattice counterparts to obtain the correct expressions for MLM parameters in terms of physical constants. The expressions indicate that the voxel size needs to be at least 0.6% larger than the molecule to accurately simulate surface reactions on triangular lattice. On square lattice, the minimum voxel size should be even larger, at 5%. We also demonstrate the ability of MLM-based schemes such as Spatiocyte to simulate a reaction-diffusion model that involves all dimensions: three-dimensional diffusion in the cytoplasm, two-dimensional diffusion on the cell membrane and one-dimensional cytoplasm-membrane adsorption. With the model, we examine the contribution of the 2D reaction pathway to the overall reaction rate at different reactant diffusivity, reactivity and concentrations.


I INTRODUCTION
Many essential and intriguing intracellular biochemical systems are mediated by the cell membrane.These systems include cell polarity establishment [1], symmetrical cell division [2], modulation of signal transduction [3] and directed cell migration [4].Spatiotemporal patterns arising from protein self-organization on the membrane [5], play a central role in these systems.The proteins self-organize primarily by reaction and diffusion processes.Membrane interactions can be classified as surface-surface reactions, where membrane-bound molecules react with each other; and volume-surface reactions, where cytosolic molecules react with membrane lipids or membrane-bound molecules.
To uncover the mechanisms underlying these systems, reaction-diffusion modeling approaches have been useful [6][7][8][9][10].In general, the choice of modeling approach depends on the time and length scales of the system [11,12].When the molecule copies are abundant and are well mixed in the surface compartment, macroscopic modeling approaches that apply rate [13] or reaction-diffusion [13,14] equation, are sufficient.If the molecule number is small or when the molecules are not homogeneously distributed in the compartment, mesoscopic methods based on the reaction-diffusion master equation (RDME) [15][16][17][18] can be employed since they account for both the fluctuations from small number of molecules and the spatial inhomogeneity across well mixed compartment subvolumes [19,20].
Although macroscopic and mesoscopic approaches are applicable for large scale simulations, the well mixed assumption imposes several limitations.These approaches for example, cannot explicitly capture the effects of space at the microscopic scale that arise from the interactions of finite-sized molecules [21][22][23], fast rebinding of reactants [24][25][26] and microscopic surface inhomogeneity such as lipid domains and membrane-associated cytoskeletal structures [21,[27][28][29][30].The spatial effects can alter not only the diffusion behavior [31][32][33], but also the reaction kinetics [22,[28][29][30]34], leading to different physiological outcomes.For example, clustering of membrane receptors changes the response of signaling network [26], fluctuation in protein copy number promotes cell polarization in the absence of spatial cue [6], rapid protein rebinding affects spatiotemporal patterns on the membrane [25] and amplifies noise during ligand interactions [35].Moreover, macroscopic and mesoscopic approaches adopt the macroscopic reaction rate constant for all reactions, which is not sufficient for irreversible bimolecular surface-surface and one-dimensional volume-surface reactions because these time-dependent reaction rates do not reach steady state [36,37].
Microscopic approaches are more suitable to model surface interactions for short timescales when the microscopic spatial effects need to be directly accounted for [38][39][40][41][42].In off-lattice microscopic particle-based methods, diffusion is simulated in continuous space with Gaussian distributed displacement.Bimolecular reactions are executed according to the Smoluchowski-Collins-Kimball (SCK) [43,44] or the Doi [45][46][47] physical model.In the former, the reaction occurs either immediately or with a probability of reflection when the distance between reactants equals to a predefined reaction radius, whereas in the Doi model, the reaction occurs with a fixed probability per unit time when the reactants are closer than the radius.Off-lattice SCK methods that support surface reactions include Smoldyn [48,49], CDS [50] and eGFRD [51], while the Doi model is adopted by MCell [52], ReaDDy [53,54] and SpringSaLaD [55].Smoldyn also recently included the option to support the Doi model.All of these methods except MCell can simulate volume occupying molecules.In a recent performance benchmark that did not include CDS [42], Smoldyn displayed the fastest simulation run time for a simple enzymatic reaction model in volume compartment.
Schemes based on the microscopic lattice method1 (MLM) [39,56] attempt to reduce the cost of resolving molecular collisions by discretizing the space into fine molecule-sized voxels.In the Spatiocyte scheme for example, a molecule only checks its destination voxel for occupancy before performing a bimolecular reaction with the occupying molecule or moving into it if it is vacant [56,57].Consequently, Spatiocyte exhibits better run time and scaling performances than Smoldyn when diffusing volume occupying molecules [56].The run time of Spatiocyte is also comparable to Smoldyn in the benchmark enzymatic reaction model [56].The reduced computational cost and the simplicity of MLM implementation have promoted its applications in both biological [57][58][59][60][61][62] and chemical [63][64][65] surface reactions.Nonetheless, biological surfaces such the cell and nuclear membranes are not arranged as fixed lattice structures.Further, since diffusion and reaction kinetics can be influenced by the lattice arrangement [66][67][68], the accuracy of MLM compared to off-lattice particle-based methods requires careful examination.Notably, a consistent approach is needed to determine MLM parameters such as voxel size and on-lattice reaction probability that can replicate the kinetics in continuous space.
In previous work, the SCK model was used together with the Spatiocyte scheme to construct a general theoretical framework of MLM for simulating reaction and diffusion processes in threedimensional (3D) space [56].Within this framework, the expressions for on-lattice reaction rate constant, reaction and rebinding probabilities, and voxel size were derived to reproduce offlattice reaction kinetics consistently.Here, we extend this framework for two-dimensional (2D) surface-surface and one-dimensional (1D) volume-surface reactions.We also employ the SCK model to derive on-lattice time-dependent rate coefficients for the surface reactions.We then obtain the expressions for the MLM parameters by equating the off-lattice rate equations with the on-lattice counterparts.In section II of this paper, we introduce the existing continuumbased reaction kinetics theory for surface reactions.In section III, we derive the expressions for surface reaction rates on lattice according to the Spatiocyte scheme and verify them using the continuum theory.In section IV, we demonstrate the applicability of the derived expressions for surface reactions that involve all dimensions.We also look at the contribution of 3D and 2D reaction pathways to the overall reaction rate.Finally in section V, we discuss the implications and limitations of this work.

II CONTINUUM-BASED REACTION KINETICS THEORY
Consider a many-body bimolecular reaction, with A and B having radii r A and r B , and diffusion coefficients D A and D B , respectively.
According to the SCK model, when the distance between a pair of A and B molecules is the sum of their radii R = r A + r B , the two will react with an intrinsic rate k a .The fraction of A remaining in the system is described by the survival probability, S irr, where [ ] denotes the concentration.When [B(0)] [A(0)], the survival probability of A is provided in the rate equation [36]: where k(t) represents the time-dependent rate coefficient.The solution for the survival probability requires the integration of the rate coefficient (Eq.2.35 in [36]): According to the particle-pair formalism of SCK model [69][70][71][72], the many-body reaction can be approximated by a simpler two-body problem.The time-dependent rate coefficient can thus be expressed as the product of k a and the survival probability of an in-contact reactant pair: Here, p reb (R, τ |R, 0) specifies the rebinding-time probability distribution for a reactive particlepair separated by a distance R at time τ , given that the pair were initially in contact.
The specific functional form of the rate coefficient depends on the spatial dimension of the reactant diffusion.The dimension is 2D for surface-surface reactions, whereas for volume-surface reactions, it is determined by the target reactant of the cytosolic molecule.The dimension is 3D when the reactant is a membrane-associated molecule or 1D when the cytosolic molecule reacts directly with the lipid membrane.For clarity, we refer to the 1D volume-surface reaction as adsorption.Since we have previously described the theory for 3D reactions [56], in the following subsections we provide the theory for 2D surface-surface reaction and 1D volumesurface adsorption.

A 2D surface-surface reaction 1 Irreversible reaction
The time-dependent rate coefficient for 2D association reaction with radiation boundary condition is given in the Laplace form as [36]: Here, k a2D is the intrinsic rate constant with dimensions of length, L and time, T, given by L 2 T −1 , and g(s) is the Green's function for 2D diffusion defined as [73] g K 0 and K 1 are the modified Bessel functions of the second kind, s = sR 2 /D, and Eq. ( 5) can thus be written as with κ = k a /D.In the limit of small s, we can approximate the modified Bessel functions: and where γ = 0.5772156 is the Euler constant.Using these approximations, the asymptotic expansion of Eq. ( 7) can be expressed as: with The corresponding long-time approximation is given as [74]: (11) where the relative error to the exact form is less than 1% at t = 100R 2 /D.Note that unlike in 3D space [56], the long-time approximation of the 2D rate coefficient does not converge to a constant term as t → ∞.

Steady-state rate constant
It is convenient to describe the 2D association reaction with a single characteristic rate constant.This is possible by defining a steady-state rate constant in terms of the mean lifetime of A, τ m [36]: k ss can be evaluated using the mean field self-consistency condition [36]: Substituting the asymptotic form of k 2D (s), as defined in Eq. ( 10), into Eq.( 13) yields Rewriting some variables in terms of the molecule area fraction, φ = πR 2 [B] and taking the small concentration limit, φ → 0 gives the following approximation Finally, the steady-state rate constant for radiation boundary condition is obtained as Similar to the 3D effective rate constant in 1/k ss3D = 1/k a +1/(4πRD), the 2D steady-state rate constant depends on the intrinsic rate k a and the relative diffusion coefficient D. Interestingly, the 2D rate constant has the additional dependency on the concentration of B.

Reversible reaction
In the SCK model for 2D reversible reaction a bound pair A-B dissociates with the rate constant k d2D (T −1 ) into A and B, separated at distance R. The survival probability of A, defined as S rev,A (t), can be calculated using the first variant of the multiparticle kernel theory (MPK1) [75,76].Although the closed form solution for S rev,A (t) in 2D is not available, it can be evaluated by numerically solving the normalized deviation defined as Here, the term is the diffusion factor function, is the chemical kinetics relaxation rate constant, and [B] = λ/k a2D is the modified concentration.F gem (s) = 1 + k a2D g(s) contains the 2D Green's function term g(s) as given in Eq. ( 6), whereas the function uses the term S irr,A (s; c 0 ), which is the Laplace transform of the irreversible reaction survival probability, S irr,A (t; c 0 ).

B 1D volume-surface adsorption
Before describing the rate for volume-surface adsorption, we first consider the 1D SCK model, where a single immobile B interacts with multiple mobile A on a filament according to Eq. (1).A can collide with B from both sides of B, while there is no self interaction among A molecules.The time-dependent rate coefficient of this reaction with radiation boundary condition is given as [36] k with κ = k a /D denoting the ratio between the intrinsic adsorption rate constant k a = k sa (unit LT −1 ) and the relative diffusion coefficient D. At long time, Eq. ( 21) behaves asymptotically as Next we consider a volume-surface adsorption system that involves an adsorbing plane at x = 0 and bulk molecules at x > 0. Initially, the molecules of concentration c 0 are distributed uniformly in the bulk and are absent on the surface.For surface adsorption process that follows the radiation boundary condition, the number of molecules adsorbed to the surface varies as (Eq.3.37 of [77]): where S is the area of the plane.The corresponding adsorption rate is well-described by the time-dependent adsorption rate coefficient k sa (t ): Note that the adsorption rate coefficient differs from the 1D SCK rate by a factor of two: , because in the latter, it occurs on both sides of the plane.At long time, the adsorption rate coefficient becomes As the bulk molecules are adsorbed to the surface, a spatial concentration gradient develops in the bulk.The spatialtemporal concentration profile of the bulk molecules C(x, t) follows Eq. 3.35 of [77]: When adsorbed molecules can dissociate from the surface with a rate k sd (T −1 ), their number varies according to (Eq.A.12 in [48]) where and

III MICROSCOPIC LATTICE METHOD
In this section, we derive the on-lattice rate coefficients for 2D surface-surface reaction and 1D volume-surface adsorption based on the Spatiocyte scheme [56,57].The particle-pair formalism of SCK model and random walk theory will be used in the derivation.In MLM, the SCK rate coefficient in Eq. ( 4) is discretized into (see Section II.C of [56] for derivation): k a is the initial lattice rate constant (see Appendix A), H n (s 0 |s 1 ) is the rebinding-time probability distribution in diffusion step n, and m is the simulation step.Rebinding here refers to the next reaction event of an in-contact reactant pair on lattice, s 0 denotes the voxel at the origin and s 1 refers to an element from the set of immediate neighbor voxels of s 0 .The rebinding-time probability is a function of random walk quantities such as P n (s a |s b ), the voxel occupation probability from voxel s b to voxel s a , that is, the probability of being at voxel s a after n steps, given that the walk started at voxel s b ; and F n (s a |s b ), the first-passage time distribution from voxel s b to s a , that is the probability of arriving at s b for the first time on the nth step, given that the walk started at voxel s a .These quantities depend on the lattice arrangement, dimension of diffusion and also the simulation scheme.The simulation step m is related to the simulation time t through the relation 2dDt = ml 2 , where d is the dimension of diffusion.In the following subsections of 2D and 1D reactions, we will derive the on-lattice rate coefficient based on the simulation scheme.

A 2D surface-surface reaction 1 Irreversible reaction
The methods presented in this work is generalized for any regular lattice arrangement, but we focus on the triangular lattice since it is used to simulate surface-surface reactions in the Spatiocyte scheme.The derivation of the rate coefficients for activation-limited (k a2D D) and diffusion-influenced (k a2D D) reactions is treated separately because the simulation scheme executes these two types of reaction in distinct manner.
In the activation-limited scheme, the generating function for the rebinding-time probability distribution H n (s 0 |s 1 ) is given as (see Appendix D.1 in [56] for derivation): where P a denotes the reaction probability.In terms of P (s 0 |s 0 ; z), the generating function of P n (s 0 |s 0 ), H(s 0 |s 1 ; z) becomes (see Appendix A for full derivation): The generating function P (s 0 |s 0 ; z) for the triangular lattice in asymptotic form is given as (see Appendix B for details): By substituting Eq. ( 33) into Eq.( 32), we obtain the following approximated form: where E = 12 exp {b 1 (1/P a − 1)} and b 1 = 2π/ √ 3. We then apply singularity analysis (Figure VI.4 of [78]) on Eq. ( 34) to obtain the large n behavior: With Eq. ( 35), we evaluate the discrete sum in Eq. ( 30) using the Euler-Mascheroni formula together with the definition of recurrence in 2D random walk: The solution in high order m terms is given by: Finally, we apply the definition of initial rate for the triangular lattice as given in Eq. (A.5) and the relation ml 2 = 4Dt to obtain the time-dependent rate coefficient: where C l = 48D exp {b 1 (1/P a − 1)} /l 2 and l is the voxel size.
In the derivation of diffusion-influenced scheme, it is convenient to work with the Laplace form of Eq. ( 30): Here Ĝ(s) is the Laplace form of the rebinding-time probability density on lattice, defined as (see Eq.D79 in [56]): where x here refers to the average time interval needed for a molecule with diffusion coefficient D x to hop across one voxel.By applying the final value theorem, we get the asymptotic form for Eq. ( 38) as Finally, by taking the small z expansion together with Eq. (33), we obtain the asymptotic rate coefficient expression: By comparing the lattice and continuum rate coefficient, we found that the asymptotic expression in Eq. ( 42) for the diffusion-influenced scheme is the same as its continuum counterpart shown in Eq. ( 10), while the time domain expression in Eq. ( 37) for the activation-limited scheme is consistent with the continuum counterpart shown in Eq. (11).To match the lattice and continuum rates, we need to impose the equality C l = C c .It then implies that the reaction probability should be chosen as where f = l/R denotes the ratio of voxel to the molecule size.Since probability P a is positive, it sets an additional constraint: To satisfy the last inequality, f = l/R has to be adapted according to the value of κ.Since κ is always positive, we only need to set a lower bound expression for the voxel size: In 3D MLM, accurate reaction kinetics requires the voxel size of HCP lattice to be larger than the molecule by l ≈ 1.02086R [56].If an HCP lattice volume compartment is bounded by a triangular lattice surface, the 3D voxel size condition would still satisfy Eq. (45).Therefore, all surface and volume voxels in the model can adopt the same HCP voxel size.
The accuracy of the lattice theory can be verified by comparing the theoretical values for the rate coefficient k 2D (t) with the simulated values.We obtained the theoretical rate coefficient from the numerical inverse Laplace transform of Eq. (41).We simulated the reaction in Eq. ( 1) with Spatiocyte at both the activation-limited (κ = 0.01 • 4π) and the diffusion-limited (κ = 100 • 4π) regimes.We logged the number of surviving A and used it to calculate the rate coefficient.The discretization of the time derivative in Eq. ( 2) gives the formula for discrete rate coefficient: where j is the index of the discretized S A and t.The boundary cases are computed as where N denotes the final time step.We compared the rate coefficients from the simulations with the theoretical values from Eq. (37). Figure 1a displays good agreements for both at activation-limited (κ = 0.01 • 4π) and diffusion-limited (κ = 100 • 4π) regimes for t t d .Next, we compared the simulated survival probability of the same reaction with the continuum-based theory, where the values are numerically evaluated according to As shown in Figure 1b, the simulated results overlap almost perfectly with the continuum-based theory, thus, confirming the accuracy of MLM.C according to the SCK model needs to satisfy the local detailed balance.This is achieved in MLM by adopting a rate constant k d2D for the dissociation reaction such that the relation is satisfied.We perform numerical simulations to confirm the ability of MLM to correctly reproduce the steady state and time-dependent behaviors in the reversible reaction.Association rates in the activation-limited (κ = 0.01) and diffusion-limited (κ = 100) cases were used in the simulation, while the dissociation rate k d2D is set to be 10 times larger than the association rate.Simulated result is compared with the MPK1 theory in Eq. (18), obtained by numerical Laplace transform.The outcome shown in Figure 2 indicates good agreement between the simulation and theory for time scales ranging from t d until equilibrium.

Generalization of MLM theory for other lattice arrangements
The expression of MLM parameter derived for triangular lattice can be generalized to other lattice arrangements that adopt MLM.In general, the variable C l in Eq. ( 37) takes the form of where b 1 and b 2 are coefficients present in the highest order term of the generating function P (s 0 |s 0 ; z): On the other hand, the reaction probability has the following general form: The expression for the probability has the following constraint on the voxel size: Here as an example, we consider the square lattice, a popular lattice choice to simulate surface reactions.The coefficients for square lattice are given as b 1 = 1/π and b 2 = 8 (Eq.A.187 in [66]).The corresponding reaction probability is with the voxel size constrained by l > 1.04722R.(55) Therefore, to recapitulate the correct continuum rate, the voxel size in square lattice has to be about 5% larger than the molecule size.This voxel size is substantially larger than the 0.6% required by the triangular lattice.The different voxel size requirements reflect the influence of lattice arrangement on the first-passage time behavior and emphasize the importance of choosing the right MLM parameters to generate accurate reaction kinetics.

B 1D volume-surface adsorption
We next formulate the on-lattice 1D rate coefficient according to the SCK model and apply the rate expression to the problem of volume-surface adsorption.
In 1D lattice, the generating function for the voxel occupancy probability from origin to origin is defined as [79] P (s 0 |s 0 The corresponding first passage time distribution, obtained by the relation F (s 0 |s 0 ; z) = 1 − 1/P (s 0 |s 0 ; z) (Eq.(I.18) in [79]), is given by Substituting Eq. ( 57) into Eq.( 31) yields the generating function for the rebinding-time probability distribution where we consider only the highest order term √ 1 − z in the limit of z → 1.The corresponding large n coefficient is obtained from the generating function according to the rule given in Figure VI.4 of [78] as: Applying Eq. ( 59) to the Noyes' rate formula in Eq. ( 30), we obtain the asymptotic form for the 1D rate coefficient: Using the definitions of initial lattice rate constant given in Eq. (A.13) and the 1D simulation step size ml 2 = 2Dt, we have the rate expression as a function of time: Note that Eq. ( 61) shares the same time-dependent form as the continuum-based theory given in Eq. ( 22).For volume-surface adsorption, the definitions for initial adsorption rate constant in Eq. (A.17) and the 3D simulation step size relation nl 2 = 6Dt are used in Eq. ( 60).The resulting adsorption rate coefficient is given as which shares the same long-time scaling behavior with the continuum-based theory in Eq. ( 25) up to the same order.In contrast to the 3D and 2D cases, the long-time expression for the 1D rate coefficient does not depend on the reaction probability and the voxel size.Figure 3: (a) Time series of adsorbed molecules simulated with irreversible (IA, triangle and circle markers) and reversible (RA, plus and square markers) adsorptions.In each case, strong (k sa = 500 µms −1 ) and weak (k sa = 50 µms −1 ) adsorption rates were tested.In the reversible adsorption, the membrane dissociation rates are k sd = 62.5 and 6250 s −1 , corresponding to the association rates k sa = 50 and 500 µms −1 , respectively.Solid and dashed lines represent the continuum-based values according to the irreversible and reversible reaction formulas in Eqs. ( 23) and ( 27), respectively.(b) The concentration profile of cytosolic A along the axis perpendicular to the adsorbing surface at x = 0 for the given time points.The adsorption is irreversible with the rate k sa = 50 µms −1 .Theoretical lines shown are according to the continuum-based theory in Eq. (26).Simulation parameters: l = 0.01 µm, D A = 1 µm 2 s −1 , and initial number of cytosolic molecule N a = 1000.
Since the long-time rate coefficient has the same form in both lattice and continuous spaces, we only need to match the initial lattice rate constant k sa with the adsorption rate constant k sa in continuum.This gives an expression for the reaction probability in terms of the adsorption rate constant, diffusion coefficient and voxel size (derivation is shown in Appendix C): To examine the accuracy of MLM in simulating the adsorption kinetics, we performed Spatiocyte simulations using the derived expression for the reaction probability.We used a large number of cytosolic A molecules in a cuboid compartment with a cross sectional area (1 µm) 2 and length 4 µm.An adsorbing plane is placed in the middle of the cuboid compartment, allowing adsorption from both sides of the surface.The number of adsorbed molecules at each time step is monitored.
Figure 3a shows the time series of A on the adsorbing plane for irreversible (adsorption only) and reversible (adsorption and desorption) reactions.Simulated results agree well with the expected values according to the continuum theories for the irreversible reaction in Eq. ( 23) and reversible reaction in Eq. ( 27).The good fit can be seen at both strongly (k sa = 500 µms −1 ) and weakly (k sa = 50 µms −1 ) adsorbing rates.To examine the spatialtemporal concentration profile, we counted the number of cytosolic molecules near the adsorbing plane in the irreversible adsorption.The resulting concentration profile along the axis perpendicular to the adsorbing plane are shown in Figure .3b.The simulation results coincide very well with the curves of continuum-based theory in Eq. (26).

IV APPLICATION OF SURFACE REACTIONS
A cytosolic molecule can react with a membrane-bound reactant via two possible pathways: it can either perform 3D diffusion in the cytoplasm and then directly react with the membranebound reactant exposed to the cytosol or it can bind first to the membrane and then perform 2D diffusion before reacting with the reactant.Both of these pathways are often adopted simultaneously in the cell.Previous works have investigated how each pathway contributes to the overall process [80][81][82][83].Here we apply the Spatiocyte scheme with the derived MLM expressions to simulate surface reactions comprising all dimensions.We study the contribution of each pathway to the overall reaction rate under the influence of different diffusivity and reactivity.
We consider a cuboid compartment of dimension H × L × L, depicting the cytoplasmic volume.The top surface of the cuboid is reflective, whereas the bottom surface represents an absorbing lipid membrane.Each of these surfaces has the area L × L. Within the system, there are two elementary species, A and B, with radius r = 0.005 µm.A c denotes the cytosolic state of A that diffuses freely in the bulk at a rate of D c .A c can reversibly associate with the membrane to become A m : The ratio of the membrane association constant over the dissociation constant is the equilibrium constant, k sa /k sd = K eq .Upon the adsorption onto the membrane, A m performs 2D diffusion at a rate of D m .On the membrane, B molecules are initialized to be immobile and randomly distributed with concentration [B] 0 .
A can react with B via the 3D pathway: or the 2D pathway: k a{2D,3D} denotes the intrinsic association rate constants for 2D and 3D reactions, whereas k r represents the dissociation rate constant.
To quantify the dominance of the 2D pathway, we measured the fraction of the 2D equilibrium rate in the total reaction rate, as in [82]: k on{2D,3D} represents the macroscopic effective rates for the 2D and 3D association reactions.The k on3D /k on2D ratio is calculated using the simulated equilibrium concentrations according to the formula which is derived by solving the rate equations for Eqs. ( 66) and ( 65) at equilibrium.We examined the dominance of the 2D pathway with changes in D c /D m , [B] 0 , and the association reaction probability, P a{2D,3D} for the 2D and 3D pathways.We fixed other variables such as the sizes of the system and molecule, k r , K eq , k a2D /k r , and the initial concentration [A c ].We used the typical cytosolic rate for D c (10 µm 2 s −1 ) with D c /D m ratio ranging from 1 to 1000.K eq = 0.15 µm and k sd = 10 s −1 are within the biologically realistic values [84,85].
From the simulation results in Figure 4, we can observe the overall decreasing trend of f 2D as the ratio D c /D m increases.The exact value of f 2D depends on the reaction probability and the concentration of reactant, [B] 0 .When the association reaction is diffusion limited (P a2D = P a3D = 1) and the reactant concentration is low ([B] 0 =100 µm −2 ), f 2D becomes more than 50% for D c /D m between 1 and 30.When D c /D m > 30, the 3D pathway becomes dominant instead.At very high [B] 0 (500 µm −2 ), the 3D pathway is dominant for all ratios of D c /D m .When the association reaction is activation-limited (P a2D = P a3D = 0.01), f 2D is still larger than 50% for D c /D m in the range [1,30], and becomes less than 50% when the ratio is higher than 30.However, unlike in the diffusion-limited case, f 2D in activation-limited reaction is less sensitive to the changes in [B] 0 .
In typical cells, membrane-associated molecules diffuse 10 to 100 times slower than their cytosolic counterpart.In such a condition, our simulation results imply the following: the 2D reaction pathway will dominate the overall reaction, provided the concentration of membraneassociated reactant is low, its diffusion on the membrane is fast and the reaction is activationlimited.Conversely, the 3D reaction pathway will become dominant when the diffusion of membrane species is slow or when the membrane-associated reactant is abundant and reacts with high probability upon collision.
V DISCUSSION AND CONCLUSION MLM surface reactions have not been verified in terms of their consistency with continuumbased theory.To address this issue, we used the theoretical framework of MLM [56] to derive the correct expressions for 2D surface-surface reaction and 1D volume-surface adsorption on lattice.By employing the SCK model and the random walk theory, we showed that the 2D lattice reaction exhibits the same long-time behavior as the continuum-based theory.After equating the on-lattice rate expression with that of the continuum theory, we obtained the formula for the reaction probability in terms of physical and lattice parameters.
Furthermore, the positively valued reaction probability imposes an additional constraint on the voxel size: it should be larger than the molecule at least by about 0.6% for the triangular lattice and by 5% for the square lattice.These constraints also meet the minimum voxel size requirement of the corresponding lattice arrangement in 3D [56].If the voxel size is exactly the same as the molecule, the simulated time-dependent reaction kinetics will deviate from the expected behavior in continuum.Such deviations should be carefully considered especially when simulating reactions containing nonlinear terms.
In 2D reversible reaction, we showed that correct equilibrium and time-dependent behaviors can be achieved by dissociating the substrate into an in-contact pair of product molecules, with a rate constant satisfying the local detailed balance.In 1D volume-surface adsorption, the long-time asymptotic behavior of MLM has the same form as in the continuum-based theory.The Spatiocyte scheme also generated spatiotemporal adsorption kinetics that is consistent with continuum theory when the correct expression for reaction probability was used.
Finally, we studied the contribution of a 2D reaction pathway in a surface reaction model with Spatiocyte simulations.We found that the dominant surface reaction pathway can be sensitive to the surface reactant concentration, intrinsic reaction rate and the relative diffusivity of reactants between the bulk and the surface.For example, the 2D reaction pathway would play a significant role in regulating the overall rate for a system that has a sparse membraneassociated reactant with activation-limited rate constants.
The main advantage of MLM when modeling intracellular reaction-diffusion processes is its ability to capture the microscopic properties of molecules directly without incurring high computational cost.As an illustration, it only takes minutes for Spatiocyte to simulate thousands of molecules with a time step of µs for a duration of seconds on a single CPU core (see performance in [56]).Spatiocyte takes physical quantities comprising molecule size, diffusion coefficient and intrinsic reaction rate as input, and generates time-series output such as molecule copy number and trajectory.
At present, Spatiocyte supports surface reactions with various geometries at the cellular scale.It has been successfully used to capture the influence of microscopic effects on the behavior of cells at the macroscopic scale.These include the formation of a high density ring over the entire bacterium cell membrane as a result of transient membrane association and rebinding of proteins [57], the clustering of proteins on the red blood cell membrane from oxidative stress [62] and the oligomerization of receptors and its influence on ligand binding kinetics [86].As the spatiotemporal resolution of imaging techniques continue to advance [87], time-dependent reaction kinetics and molecular trajectories will become more accessible.These high resolution experimental data coupled with efficient microscopic simulation techniques such as MLM will provide a complementary way to investigate mechanisms underlying various biological reactiondiffusion processes.
The uniform voxel size adopted by MLM reduces computational complexity and conse-quently, contributes to its low computational cost.However, in surface systems requiring realistic simulation of distinct-sized molecules with non-spherical structures, additional considerations would be needed for MLM to be applied.One potential solution is to reduce the voxel size and let a single molecule occupy more than one voxel according to its size and shape.Alternatively, we can represent molecules with distinct shapes and sizes off-lattice and perform hybridized simulation with on-lattice molecules.The implementation and accuracy of such schemes compared to fully off-lattice methods would require further examination.Another future milestone for MLM is to establish and verify its consistency in highly crowded environments.The on-lattice rate has to be reformulated to account for the many-particle interaction.
The resulting lattice theory should then be compared and matched with the continuum-based theory.

B Voxel occupancy probability on triangular lattice
The voxel occupancy probability from origin to origin, P n (s 0 |s 0 ) for the triangular lattice is given as [66,88]:

Figure 1 :
Figure 1: Comparison of on-lattice simulations with on-and off-lattice theories for surface-surface reaction A + B − → B. (a) Simulated on-lattice time-dependent rate coefficients (solid lines) compared with on-lattice MLM theory in Eq. (41) (dashed lines).For better visualization of the time-dependent behavior of the two extreme cases, the simulated and theoretical lines are normalized by the initial theoretical value.(b) Simulated on-lattice survival probability of A (points) compared with off-lattice SCK theory in Eq. (48) (solid lines).Activation-limited (κ = 0.01 • 4π) and diffusion-limited (κ = 100•4π) cases are indicated by the top and bottom lines respectively.Simulations were performed with Spatiocyte and the following parameters: Area = (6.5 × 6.5) µm 2 , R = 0.01 µm, l = 0.01 × 1.0209 µm, D A = 1, D B = 0 µm 2 s −1 , N a = N b = 423, duration = 0.2 s, logging interval=50t d