Concentration dependence of diffusion-limited reaction rates and its consequences

Diffusion-limited association reactions are ubiquitous in nature. They are particularly important for biological reactions, where the reaction rates are often determined by the diffusive transport of the molecules on two-dimensional surfaces, such as the cell membrane. The peculiarities of diffusion on two-dimensional surfaces may lead to nontrivial reaction kinetics, such as concentration dependent rate of association between two molecules. However, traditionally, the kinetics of biomolecular association reactions has been modeled using the law of mass action, which assumes that the rate of reaction is a concentration independent constant. In this paper, using multiscale molecular simulation, we investigate the concentration dependence of diffusion-limited association reactions on 2D surfaces. In particular, we quantify the influence of short-ranged pair interactions on the concentration dependence of the reaction rates and codify it in an empirical law. Using this law in a chemical kinetic model, we find that the the steady state behaviors of simple chemical systems are drastically modified by the presence of concentration dependent rates. In particular, we find that it leads to suppression of intrinsic noise in dimerization reaction and destabilizes robust oscillation in Lotka-Volterra predator-prey systems. In fact, we see a transition from robust to fine-tuned behavior in the latter. In addition, we show that concentration dependent reaction rates arise naturally in stochastic predator-prey systems due to intrinsic noise. We comment on the consequences of these results and discuss their implications in the modeling of complex chemical and biological systems. In particular, we comment on the range of validity of the law of mass action, which is a staple in all theoretical modeling of these systems.


I. INTRODUCTION
Association reactions are a type of elementary reactions, in which two or more reactant molecules form one or more product molecules. Formation of a dimer from two monomers is an example of an association reaction, as is the formation of water from hydrogen and oxygen [1]. Due to its elementary nature, association reactions play a central role in many physicochemical processes, including pattern formation [2,3], aggregation [4], and cell signaling [5][6][7][8][9]. A typical association reaction involves two steps [10][11][12][13]. In the first step, two molecules are transported near each other through some transport processes. Once the two molecules encounter each other, in the second step, the molecules interact with each other with an intrinsic association rate κ I to form the product molecule. It is usually assumed that the second step is much slower compared to the first step, such that an association reaction occurs after many encounter events. These assumptions have two repercussions: (a) one can assume that the reactant molecules are well mixed, so that the reaction rate is solely determined by the interaction step and (b) the formation of the product molecule follows a Poisson process, such that the rate of association reaction is a time and concentration independent constant [14]. When combined, these two observations lead to the celebrated Law of Mass Action (LMA) [1,15], which states that the propensity of an association reaction is equal to the product of the constant reaction rate, k 0 , and the mass action Φ, where the latter is the total number of possible reactant pairs. In particular, if we consider the following association reaction then LMA states that the propensity, r, of the association reaction is: where [.] denotes the concentration of a molecule, and Φ = [A] [B] is the mass action. LMA is strictly valid when the rate k 0 is a concentration independent constant. For example, LMA works well for dilute solutions, where this "rate law" was originally developed. However, its range of applicability did not stay confined within just the purview of dilute solutions. Its simplicity and the ubiquity of association reactions have led to its application in disparate problems [16] with varying amount of success. However, it is unclear whether some of the assumptions underlying LMA are still valid for these systems. For example, most cell-signaling reactions occur on two dimensional cell membranes, where, often, the intrinsic reaction rates are higher or comparable to the diffusive encounter rate, such that the reaction kinetics depends crucially on the diffusive transport of the molecules. In particular, the peculiarities of diffusion in two dimension (2D), including nonzero probability of encounter between all particles [17,18], can result in concentration dependent diffusive encounter rates, which, in turn, results in a concentration dependent k 0 , which invalidates the application of LMA in such situations [19,20]. Furthermore, a fundamental assumption of LMA is that the molecules are essentially point particles that interact with each other only when they collide with each other [21]. In reality, most molecules interact with each other through finite, albeit short-ranged, interactions, which may also influence the diffusion-limited reaction rates. In fact, in such situations also, the diffusion limited reaction rate, κ, may become concentration dependent [14,19,22] and violate the assumptions of LMA.
A practical workaround to this challenge is to use an adaptive, concentration dependent, rate of association that varies with time [10,17,18,23,24]. Use of such adaptive reaction rates in chemical kinetic models drastically improves the prediction of the transient kinetics. Unfortunately, the functional form of the adaptive rates used in these studies are complex and are not immediately conducive to analytical treatments. Hence, in this paper, we offer a complementary empirical formulation of the concentration dependence, which allows us to do analytical computations of the steady state properties. The feasibility of the analytical treatment offers practical advantages to explore the consequences of the concentration dependence for system parameters that are difficult to explore through simulation.
To theoretically study the concentration dependence of the diffusion-limited reaction rates and to construct an empirical law at concentrations relevant to most applications, one has to simulate the transport and interaction of the molecules in large spatially heterogeneous systems. The main challenge to such line of enquiry is that diffusion is difficult to investigate through molecular simulation. In this paper, we overcome this challenge by using a recently developed multiscale simulation framework, called the Green's Function Reaction Dynamics with Brownian Dynamics or BD-GFRD [13,[25][26][27][28]. We combine BD-GFRD with chemical kinetic models to construct a hierarchical multiscale simulation framework (Eq. 3). In the first level of hierarchy of this framework, using BD-GFRD, we measure the concentration dependence of the rates and codify it in a functional form, κ(Φ) (Fig. 1A-i,ii). In the next level, we use κ(Φ) in a chemical kinetic model to study the behavior of the chemical systems at timescales that are not reachable through BD-GFRD( Fig. 1A-iii). We ask, can we still use LMA even when the reaction rates are concentration dependent? In the two model systems that we study here, we find that steady state properties obtained from concentration dependent rates are qualitatively different from LMA, but under some special conditions this difference is negligible. Because chemical kinetic models are staple in many lines of scientific enquiry [29], our results may provide useful guidelines and design principles for these models.

A. Molecular simulation of dimerization reactions
To study how short-ranged pairwise interaction leads to concentration dependent reaction rates, we use a hierarchical simulation framework, which interfaces a recently developed multiscale molecular simulation method called BD-GFRD [13,[25][26][27][28] with chemical kinetic models to study association reaction kinetics at wide spatiotemporal scales and broad concentration levels ( Fig. 1A-D and Eq. 3). Please check appendix B to find the detailed description of the methods.
To measure the impact of short-ranged interactions (Fig. 1C) on diffusion-limited reaction rates, we investigate the concentration dependence of the rates of the following reactions: We measure the diffusive encounter rates of these reactions in an ensemble in which the number of each species of molecules is conserved. Using this set-up in the BD-GFRD [13,[25][26][27][28] simulations (see appendix B), we measured the reaction rates by varying the concentration of the monomers in the range of 2 − 300/µm 2 . These concentrations are similar to the concentration of integral and peripheral membrane proteins on cell-membranes.
We measured probability distribution, ψ D (t D ), of the time interval, t D , between two consecutive homodimerization reactions. It has three distinct regimes separated by two timescales: the rebinding timescale τ rebind and the timescale above which the reaction events are described by an exponential decay: τ bulk (Fig. 2 A). For t D < τ rebind , ψ D (t D ) decays as t −1 D . Such dependence occurs due to rapid rebinding of the monomers to form dimers [13]. For τ rebind < t D < τ bulk , ψ D (t D ) decays as t −0.2 D . We find that this region is present in all the monomer concentrations considered here (Fig. 2 B). As the concentration of the monomers decreases, the width of this region also decreases, which suggests that this unusual scaling results from reactions that occur before the particles have lost their memory of the previous encounter. Beyond this observation, we do not understand this scaling well and plan to investigate it in a future paper. For timescales above τ bulk , ψ D (t D ) decays exponentially ( Fig. 2A). In particular, we find that rescaling the time by τ bulk collapses all distributions on to a single master curve (Fig. 2C), which implies that the propensity of the reactions r is inversely proportional to τ bulk . τ bulk depends on the concentration of the monomers (Fig. S1) and at low Φ is given by: where L = 1000nm is the system size, a = 2r A , D A = 1µm 2 s −1 is the diffusion constant of A, and  Figure 1. Hierarchical multiscale simulation framework: (A) We combine molecular simulation with chemical kinetic model to study association reaction kinetics at broad spatiotemporal scales. (i) Using BD-GFRD, we study how short-ranged pairwise interaction modifies the rate of association reaction between two molecules on a 2D periodic box of size 1µm × 1µm. (ii) The measured rate, κ is concentration dependent and its dependence on concentration is represented by the function κ(Φ), where Φ is the product of the concentrations of the reactants, a.k.a the mass action. (iii) κ(Φ) is used in a spatially homogeneous chemical kinetic model to study concentration variations at timescales not reachable through molecular simulations (> 1s). (B) The molecules are represented as disks of radius rA and rB, and (C) they interact through either Lennard-Jones (LJ) or WCA interaction. (D) The molecules react as soon as they are closer than the reaction radius ra. If ra < rmin (see C and D for definition of rmin ), there is strong excluded volume interaction between the reactant molecules. The light blue region surrounding a molecule denotes the region over which excluded volume interaction is felt by another molecule. More precisely, when two light blue regions touch, the distance between the two molecules is rmin. [18]. We observe the same behavior for the heterodimerization reactions as well (not shown). Crucially, we observe identical ψ D (t D ) in a set of simulations in which [A] + 2[A 2 ] was kept fixed, but concentrations of [A] and [A 2 ] could vary, which consolidates our observation (Fig. S2). One should note that, although ψ D (t D ) is described by an exponential decay in this regime, it does not imply that this process is Markovian. In 2D diffusion, due to the finite probability of reencounter between two reactants, there is not a single well-defined rate of reaction for all concentrations. Instead, the rate depends on the concentrations of the reactants, which can be used in a manner similar to the Markovian rate constant [10,17,18,23,24].

B. Diffusion-limited bulk dimerization rate
The propensities, r, of the dimerization reactions separated by t D > τ bulk can be computed by fitting the tail of ψ D (t D ) with an exponential function, from which we can determine κ(Φ), the "bulk" concentration dependent diffusion-limited reaction rate, using the following formula: For LMA to hold, κ(Φ) should be independent of Φ. However, we find that κ varies with Φ for all pair interactions considered here (Fig. 3A). Irrespective of the type of interaction potential, the functional form of κ(Φ) is nearly identical for A 2 , except its value near Φ = 0. The heterodimerization reaction, i.e. the formation of AB, also has similar concentration dependence. However, κ(Φ) varies weakly compared to A 2 . Finally, κ(Φ) also depends on the reaction radius r a [30,31] (Fig. 1D). We found that the lower the r a , the more slowly κ(Φ) increases with Φ. In fact, when r a = 0.89r min ≈ σ, κ(Φ) does not change appreciably for the Φ values considered in our simulations (Fig. 3B). In fact, for all the parameters considered here, as Φ → 0, κ(Φ) converges to a concentration independent value κ 0 , which is consistent with previous works [14,18,32,33].

C. Empirical rate law
The preceding discussion implies that, in general, κ(Φ) can be expressed in terms of κ 0 , the value of κ(Φ) as Φ → 0, and a function of the mass action Φ. The latter can be obtained by measuring the variation of κ(Φ) − κ 0 with Φ. Remarkably, we find that κ(Φ)−κ 0 displays universal variation across different reactions studied here. In particular, for the concentrations of the monomers studied here, we find that κ(Φ) − κ 0 ∼ Φ , such that κ(Φ) ≈ κ 0 + κ 1 Φ , where = 3/8 = 0.375 (Fig. 4A). Furthermore, we find that there is a mass-action value Φ 0 such that κ 1 Φ 0 = κ 0 , i.e., Φ 0 measures how quickly a reaction deviates from the concentration independent behavior; the larger the value of Φ 0 , the slower is the deviation. Combining the above definitions, κ(Φ) can be succinctly written as: This form of κ(Φ) is valid for very low packing fractions of the reactants. For example, the packing fraction of the molecules at 300/µm 2 , the highest concentration of molecules studied here, is approximately 0.3%, which is much lower than the packing fraction of molecules in a cellular environment (up to 40% [20]). Therefore, we expect that this approximate form of κ(Φ) will not hold at such high packing fractions. Instead, because of the increase in viscosity observed at high packing fractions, we expect the diffusive encounter rate κ(Φ) to saturate to a maximum value, κ max [34]. It is possible to obtain the value of κ max by doing direct simulation of the molecules at high density. Although we plan to study such a system, we do not attempt to do it in this paper. Instead, we use Smoluchowski's theory of diffusion limited reactions [18] to obtain κ max , which for the parameters of our system is approximately 20µm 2 /s (appendix E). Furthermore, because κ(Φ) is a convex function for small Φ and saturates to a maximum value at higher Φ values, we expect its most general form to be sigmoidal. Therefore, we postulate the following empirical law that describes κ(Φ) for concentrations beyond the values simulated in this paper: It is easy to check that this sigmoidal κ(Φ) has the desired behavior for different limits of Φ. Also, it fits the observed rates well (Fig. 4B). In fact, the functional form that we obtain here is very similar to what has been described earlier [18]. The advantage of our approach is that the concentration dependence of the encounter rate is entirely codified by the flux Φ, whereas everything else are concentration independent parameters, and can, in principle, be obtained by using self-consistent theories of association reaction kinetics [10,17,18,23,24]. Empirical law in this simple form is not only useful as an input to the numerical methods, such as Gillespie algorithm, it is also conducive to analytical treatments, such as the calculation of the steady state concentrations and their stabilities, as described in the next section.
In the rest of this paper, we study the consequences of this sigmoidal rate law using chemical kinetic models (see appendix B). For the sake of brevity, we shall refer to this rate law as the "law of concentration dependent rates" or LCDR, ala LMA. We study two model chemical systems using both LMA and LCDR and compare the steady state behavior of these models under the two different rate laws. We should point out that ours is by no means the first such attempt at using chemical kinetic model to study the long-term behavior of chemical systems. In fact, there has been several papers that have pointed out that chemical kinetic equations are not accurate at low concentrations [18,24]. However, in spite of this limitation, the chemical kinetic models with concentration dependent rate constants have been shown to reproduce kinetics of association reactions with nearly exact transient kinetics [18,24]. Therefore, we are justified in using such an approach to study the consequences  of the concentration dependent rates.

A. Suppression of intrinsic noise
We studied the consequences of the concentration dependent reaction rates on the simplest association reaction possible: the homodimerization reaction.
where k on = κ(Φ) is the association rate; κ max = 20µm 2 s −1 and k of f is the dissociation rate of the dimer. We study the steady state properties of this reaction using a chemical kinetic model. To do so, we measure the propensities of the association reactions using Eq. is given by Eq. 9. We study the steady state properties of this system by varying the dissociation rate, k of f . To obtain equilibrium steady states, the rates have to obey detailed balance, i.e., k of f has to be equal to k on k D , where k D is the equilibrium constant. To compute the mean equilibrium concentration and fluctuation at different k D values, we simulated the dimerization reactions using Gillespie algorithm. The mean equilibrium concentrations, < [A] ss >, and the variance of the fluctuations in the steady state, σ 2 [A] , are identical between LMA and LCDR ( Fig. 5A-B), which is what we expect. To verify the correctness of the simulated results, we also calculated the steady state concentrations analytically for both LMA and LCDR (see appendix C). We found that the mean monomer concentration < [A] ss > obtained  from analytical calculation and simulation are identical (Fig. 5A).
To study the same quantities in a nonequilibrium steady state, we use k of f = k D . For LMA, for which κ(Φ) is a concentration independent constant, this choice of the dissociation rate does not break detailed balance and we reproduce the equilibrium mean concentration and fluctuation profile. In contrast, for LCDR, this functional form of k of f is sufficient to break detailed balance and obtain a nonequilibrium steady state. In addition, we discover that for same steady state concentrations, the variance of the steady state fluctuations are smaller in the nonequilibrium steady state of LCDR than in equilibrium state under both LMA and LCDR kinetics (Fig. 5B). This is a remarkable result and can be explained through a simple analogy. To do so, we first write κ(Φ) in its ap-proximate form for Φ near Φ 0 : Written in this form r = κ(Φ)Φ can be interpreted as a cost function, which computes the cost associated with maintaining the system at a particular mass action value Φ. LMA, for which κ 1 is zero, imposes a penalty κ 0 ∆Φ for a concentration fluctuation that changes Φ to Φ+∆Φ. On the other hand, LCDR imposes an additional penalty of (1 + )κ 1 ∆ΦΦ . In equilibrium, this additional cost is balanced by the reverse reaction, whose propensity is Hence, the fluctuation profiles are identical between LMA and LCDR in equilibrium. In contrast, in the nonequilibrium state, at a given value of Φ, the reverse reaction imposes a constant cost, given by k D [A 2 ]. Therefore, for all Φ > 0, fluctuations cost more in LCDR than in LMA. In fact, for LCDR, the cost of fluctuation increases with increasing Φ, but, for LMA, the cost of fluctuation is independent of Φ. That is, the difference in variance should increase with increasing monomer concentration, which is what we observe in Fig. 5 until < [A] ss >≈ 150µm −2 , beyond which the finiteness of the system becomes important.
We must stress that the attainment of the nonequilibrium steady state was possible due to the concentration dependence of the reaction rates (see appendix C). As we have shown, in the absence of the concentration dependence, as in LMA, the system always reaches chemical equilibrium, unless the reverse reaction is timedependent. However, the reader should note that our study of the effect of concentration dependent reaction rate on the dimerization reaction is a toy model through which we have decided to expose the impact of the concentration dependent rate. To assess the practical relevance of our results, we still need to verify whether, in real systems, detailed balance is broken in the manner described in this paper. However, as the preceding analogy shows, the underlying reason for the noise suppression is the additional cost of fluctuation due to the concentration dependence, which does not depend on the choice of the nonequilibrium steady state. Therefore, the suppression of intrinsic noise by concentration dependent rates may be a useful strategy for controlling noise.
Control of noise is an important aspect of many biological functions [35][36][37][38][39][40]. It is possible that cellular systems may control noise by modulating the local concentration of reactants on 2D surfaces, such as membranes, which can increase or decrease the amount of noise. Indeed, liquid-liquid phase separation of signaling proteins has been postulated to reduce noise in the processing of the signals [41,42]. Therefore, we anticipate our results will provide important insights in deciphering the nature of noise control through such biological mechanisms.

Lotka-Volterra model in the presence of concentration dependent rates
Next, we consider the consequences of the pairwise interaction on the behavior of chemical oscillators. We study the Lotka-Volterra (LV) predator-prey equations [16], which describes the dynamics of the following set of reactions: where X is the prey and Y is the predator, and φ represents death. The most general LV equation can be written in the following way: which admits two steady state solutions: x = y = 0 and ,y) , where x and y are the concentrations of X and Y , and β(x, y) = β(xy) is the concentration dependent predation rate. We assume that β(xy) has the same functional form as κ(Φ) with Φ = xy and Φ 0 = 10 5 µm −4 .The concentration of the nonzero steady state depends on the kinetics used. For LMA, i.e. β(x, y) = β 0 , the steady state concentrations increase linearly with α and γ. In contrast, for LCDR, the steady state concentrations are highly nonlinear when either α or γ is large, but linear behavior is restored when both α and γ are small or when the steady state mass action Φ ss >> Φ 0 (Fig. 6A) .
While the x = y = 0 steady state is a saddle point for both kinetics, the stability of the nonzero steady state depends on the underlying kinetics. For LMA, i.e. β(x, y) = β 0 , the steady state is oscillatory with frequency √ αγ for all values of α and γ. That is, the oscillation of the prey and the predator population is robust to the variation of the parameters. In LCDR, we do see similar behavior when α, γ < 10 −2 or when α, γ > 10 5 , where κ(Φ) is effectively Φ independent. However, outside these parameter ranges, the stability of the steady state depends on α and γ, as can be seen from the eigenvalues at different parameter values (Fig. 6B). In particular, when α > γ, the steady state is a stable fixed point. When α < γ, the steady state is an unstable fixed point and admits exponential rise of the peak predator population and exponential fall of the prey population. Finally, only when α = γ, the steady state solution is periodic (Fig. 6C). Therefore, in LCDR, the parameters have to be fine tuned to achieve sustained oscillation. The frequency of the oscillating state is higher in LCDR than in LMA. A perturbative analysis of the LV equations supports this observation (appendix D). In fact, its prediction agrees well for a range of values of α (Fig. 6D). This analysis also predicts that the oscillation in LV predator-prey model is extremely sensitive to concentration dependence of rates: even weak concentration dependence is sufficient to break the robustness of the oscillation. This result indicates that very special conditions are required for oscillatory reactions to remain stable, such as finite intrinsic reaction rate, κ I , or presence of intrinsic or extrinsic noise.

Stability of steady states in the presence of finite intrinsic reaction rates
So far, we have only considered the effect of the concentration dependent reaction rate in the diffusion limited regime. In this regime, the intrinsic reaction rate, . Phase diagram of LV oscillators when the intrinsic reaction rate is finite for different intrinsic rates κI at D = 1µm 2 /s. We have used the classification scheme provided in [18] to designate different regimes. We observe that finetuned oscillation persists even when the reactions are limited by the intrinsic reaction rates, although the corresponding parameter range narrows as κI is decreased. These results suggest that effect of diffusion may be felt even when the reaction rate is controlled almost entirely by the intrinsic rate κI .
κ I → ∞. However, in biological or physicochemical systems, we often encounter situations in which κ I is finite. Under such circumstances, the reaction kinetics may be entirely determined by the intrinsic reaction rates, which in turn, will restore the robust oscillatory behavior of the LV model. Therefore, we sought to understand how finite κ I modifies the phase diagram obtained in the diffusion limited regime (Fig. 6C). To do so, we define the effective predation rate β ef f : where β(x, y) is the concentration dependent encounter rate. We use this effective predation rate to calculate the steady state solutions of the LV equations and the stability of the solutions. Following the classification described in [18], we investigated four different κ I values for which the transient kinetics are diffusion controlled (100µm 2 /s), diffusion influenced (10µm 2 /s), weakly diffusion influenced (0.1µm 2 /s), and intrinsic reaction rate limited (0.001µm 2 /s). We find that fine-tuned oscillation persists in all four cases, even when the transient kinetics is rate limited and diffusion plays a negligible role (Fig. 7). It happens because the stability of the steady state behavior is not only determined by the effective predation rate, but also its derivative at the steady state (see appendix D). The derivative has non-negligible, albeit diminishing, influence even when κ I << D, which leads to fine-tuned oscillation even when β ef f is essentially concentration independent. Because of the diminishing influence of the derivative, the parameter space for fine-tuned oscillation shrinks with decreasing κ I . This result suggests that the influence of the diffusive transport process can impact steady state behavior even when it has negligible impact on the transient kinetics. In particular, in the context of our results, we find that fine-tuned oscillation is a robust behavior that persists even when the intrinsic reaction rates are several orders of magnitude smaller than the diffusion coefficient. The presence of intrinsic or extrinsic noise is an important determinant of the steady state behavior, which we have not considered so far in this analysis. It is wellknown that noise can significantly change the behavior of chemical oscillators [43][44][45][46][47]. Its impact on the canonical LV systems, systems without concentration dependent rates, have been widely studied. In particular, applying a chemical master equation for the birth-death processes to the LV reactions (Eq. 12), it can be shown [48] that in the presence of intrinsic noise the mean concentrations of the predator and the prey follow modified deterministic equations: where G(x(t), y(t)) = cov(x(t), y(t)) is the (timedependent) covariance between the prey and the predator population. This equation has the same form as Eq. 13, with β 0 xy + β 0 G(x, y) acting as the concentration dependent predation propensity, β(x, y)xy. In particular, we find that for typical stochastic trajectories obtained from the canonical LV model (Fig. 8 A), G(x, y) has nonlinear dependence on Φ. To see this, we define g(x, y) = G(x, y)/Φ, which plays the same role in this context as κ(Φ) − κ 0 does in the context of diffusioninduced concentration dependence. As Fig. 8B shows, g(x, y) has a complex dependence on the mass action Φ, which suggests that intrinsic noise may lead to the emergence of concentration dependent reaction rates. We must stress that the concentration dependence of g(x, y) arose even when the predation rate β 0 was concentration independent. That is, the emergence of concentration dependent predation rate is not limited to dimensions less than three, and we can expect to observe its effect even in 3D. Moreover, even though the functional form of g(x, y) is very different than β(x, y), we expect the same drastic shrinkage of the parameter space available for stable oscillation of the concentrations, since presence of any concentration dependence leads to finetuning of the oscillation (appendix D). This observation may explain why it is difficult to obtain stable oscillation in stochastic LV predator-prey model. However, more careful studies are needed to establish this claim.
Finally, here we have shown the effect of intrinsic noise only for cases when β is not concentration dependent. It will be interesting to study the combined impact of the concentration dependent predation rate β(x, y) and the emergent rate g(x, y) on the steady state stability of LV model and other dynamical systems. Moreover, it is likely that similar concentration dependent rates will arise in biologically relevant situations also, which are noisy chemical systems. Therefore, a careful classification of the biological systems in terms of concentration dependence of the reaction rates is necessary. A classification scheme, in terms of accuracy of well-mixed reactions, has been suggested in ref [18]. However, as we have shown, we also need to consider the impact of diffusion and noise on such a classification.

IV. DISCUSSION
Concentration dependent reaction rate In this paper, we have used a hierarchical multiscale simulation framework to identify the origin of concentration dependence of the rate of association reactions and its consequences.
We find that diffusion in 2D leads to the observed concentration dependence, and this dependence is universal across the different types of association reactions and the pair potentials considered here. In particular, the concentration dependent rate κ(Φ) can be written in a simple empirical form, which leads to drastic changes in the steady state properties of model chemical systems. κ(Φ) is characterized by four parameters κ 0 , κ max , , and Φ 0 . When Φ → 0, κ(Φ) → κ 0 , which is a concentration independent constant. κ 0 depends on the interaction potential, the reaction radius, and the diffusion coefficient. However, the diffusion coefficient merely rescales the value of κ 0 (Fig. S3), whereas the interaction potential and the reaction radius has more substantial impact. κ 0 depends on the interaction potential: κ 0 is higher for attractive LJ interaction than repulsive WCA interaction (Fig. 3B) for otherwise identical systems. Furthermore, we find that the lower the reaction radius the lower the value of κ 0 , which is what we expect (Fig. 3B). In particular, when r a = 0.89r min ≈ σ, the value of κ 0 is comparable to its value computed using Smoluchowski's theory of diffusion limited association reaction, which uses purely collision-based interaction between molecules [14,18] (appendix E). For all the analysis involving κ(Φ), we have assumed that D = 1µm 2 /s and r a = 1.2r min , such that all κ 0 values are around 4µm 2 /s. While the empirical laws fitted in Fig. 4B do depend on the value of κ 0 , it does not make it less general. In fact, computation expense permitted, we can repeat the same set of analysis using r a = 0.89r min , which will have a different κ 0 and Φ 0 values, but same values of κ max and . κ max is the value of κ(Φ) in the other extreme, when Φ → ∞. As we have discussed in section II.C, the value of κ max can be measured by running simulations at higher packing fractions of the molecules, which we will do in a future paper. At higher packing fractions, the diffusive encounter rate is affected by the dynamic variation of the viscosity of the reactant solution. The variation of the viscosity depends on the interaction potential [34,49]. Therefore, we believe, the choice of attractive or repulsive interaction may lead to different κ max values. However, in this paper, we ignore this difference. Instead, irrespective of the interaction potential, we have used a meanfield formula derived using Smoluchowski's theory to estimate the value of κ max , which is around 20µm 2 /s (appendix E).
Between κ 0 and κ max , κ(Φ) has a sigmoidal variation that is characterized by an exponent and a mass action scale Φ 0 . Through data analysis, we find that is equal to 3/8 = 0.375. While we cannot provide a fundamental reason behind this exponent, the meanfield theory (appendix E) predicts that for low Φ values κ(Φ)−κ 0 indeed scales as Φ 0.4 , which is remarkably close to what we have found through our analysis. Therefore, it is likely that the Φ scaling emerges due to the properties of diffusion in 2D.
The scale Φ 0 , in effect, measures the contribution of the concentration dependent rates on the reaction kinetics.
If Φ values are much lower than Φ 0 , the concentration dependent rate is subdominant and the behavior of the chemical system in this regime is well approximated by LMA. However, for Φ values near or higher than Φ 0 , the reaction rates are strongly concentration dependent until they reach κ max , beyond which the rates become concentration independent again. Crucially, Φ 0 is system dependent and higher values of Φ 0 imply broader range of validity of LMA. For example, for the heterodimerization reaction, i.e., A + B → AB, Φ 0 value is higher than that for the homodimerization reaction. Therefore, we expect LMA to describe association reactions between unlike species much better than between like species. We expect Φ 0 to depend nontrivially on the form of the interaction potentials. In this paper, we have only considered isotropic interactions of identical strength and measured their Φ 0 values by fitting the functional form of κ(Φ) to the observed data. However, real systems often interact with each other through anisotropic interactions of varying strength. Understanding how κ(Φ) and Φ 0 behaves in these systems will help establish the range of validity of LMA and render chemical kinetic models more accurate. In particular, it may be possible to estimate the values of κ 0 and Φ 0 using self-consistent theories of association reaction in 2D [10,17,18,23,24].
Dimensionality of concentration dependence In this paper, the concentration dependence of the reaction rates stems from the peculiarities of diffusion in 2D. The reentrant nature of diffusion in dimensions lower than three leads to the observed concentration dependence, which, as we have exposed, leads to drastic changes in the steady state properties of dynamical systems. However, concentration dependence can originate from other factors as well, which we have not considered here. For example, viscosity can change dynamically at molecular packing fractions that are quintessentially found inside a cell. Such dynamic changes in viscosity can dramatically reduce both translational and rotational diffusion coefficients [34], which can give rise to another form of concentration dependence of the diffusion limited reaction rates. Such crowding induced concentration dependence is not limited by the dimension of the problem and can be observed even in 3D [20]. Furthermore, as we have shown, albeit non-rigorously, intrinsic noise can also lead to emergent concentration dependence of the association reaction rates, even when the original reaction rates are concentration independent constants. It is well known that diffusion limited reaction rates are concentration independent constants only for dimensions greater or equal to three. Therefore, intrinsic noise offers another mechanism to obtain concentration dependent reaction rates in 3D. We expect that these three different flavors of concentration dependence will have different impact on the transient kinetics and may even have different steady state behaviors. However, all of these claims warrant rigorous examination, which we hope will inspire many future investigations of the origin and the consequences of the concentration dependent rates.
Impact on biomolecular systems Concentration dependence of the reaction rates leads to drastic changes in the behavior of simple chemical systems, which may have serious repercussions on the behavior of more complex chemical systems, such as the cell-signaling or metabolic reactions. We find that concentration dependent reaction rates lead to reduction in intrinsic noise in dimerization reaction. Crucially, they destabilize robust oscillation in the LV predator-prey system and renders it fine-tuned to the parameter values. Although LV model is less relevant for biomolecular systems, the results derived from this model system can be directly applied to study feedback driven oscillatory systems on cell membranes or other 2D surfaces inside the cell. For example, signaling circuits in cell growth and development do show oscillatory behavior [8,9,29,50,51]. These circuits are usually membrane bound. Therefore, we expect our results to be directly applicable to such systems.
The concentration dependent rates have another important implication in the context of modeling of biomolecular systems. The systems biology models used to study complex biological functions, such as cell signaling or metabolism, often use reaction rates measured from experiments. More often than not, the measured reaction rates vary wildly over several orders of magnitude [52]. Although such variations are often attributed to experimental imprecision, our results offer a plausible alternate hypothesis. As we have seen, the concentration dependent rates can vary across several orders of magnitude. Also, it is well-known that cells modulate the concentration of the biomolecules to achieve different tasks. It is possible that the rates measured by different experiments had different concentration of the reactant molecules, which led to the broad variation of the measured rates. This hypothesis has two implications. Firstly, it implies that the larger the range of reaction rates, the larger is the concentration fluctuation of the reacting molecules. Secondly, it will be possible to fit a sigmoidal κ(Φ) using the measured rates and use that concentration dependence to construct a class of systems biology models. Given the concentration dependent reaction rates have such drastic impact on the simple systems, it will lead to very interesting and potentially undiscovered phenomena in complex chemical systems.
Acknowledgement This work was funded by an LDRD grant (XX01) from LANL and has the following LA-UR number:LA-UR-20-24037. Computations used resources provided by the LANL Institutional Computing Program, which is supported by the U.S. DOE National Nuclear Security Administration under Contract No. DE-AC52-06NA25396. The author thanks Prof. Angel E. Garcia, Dr. Van A. Ngo, Prof. Yair Shokef, and Dr. Sumit Majumder for useful discussion and critical reading of the manuscript. The author also thanks two anonymous referees whose inputs have significantly improved the quality of this paper.
The basic tenet of GFRD is that tagged-particles, such as proteins, in biomolecular systems can be partitioned into two groups: (a) isolated particles, which freely diffuse and (b) interacting particles, which interact with other particles. In a special form of GFRD, the isolated particles are propagated in an event-driven manner using the Green's function of diffusion until the particle encounters another particle and it can no longer be treated as an isolated particle, whereas the interacting particles are propagated using a molecular mechanics algorithm, such as molecular dynamics, dissipative particle dynamics or Markov state models [25]. In the present work, we use overdamped Brownian dynamics (BD) as our molecular mechanics algorithm. This updated form is called BD-GFRD [25,26].
We assumed that the particles of type A are spheres of radius 2 nm and particles of type B are of radius 2.9 nm (Fig. 1B) and they interact with each other through short-ranged isotropic interactions, such as the 6-12 Lennard-Jones or WCA interaction, (Fig. 1C) of strength 6k B T , where k B is the Boltzmann's constant and T = 310K is the temperature of the heat bath. In addition, the diffusion costants of particles of type A and type B were 1µm 2 /s and 0.69µm 2 /s, respectively. Because we assume overdamped dynamics, their mass is irrelevant for the simulation. These parameter values correspond to typical protein-protein interaction parameters, e.g. Ras-Raf interaction [32]. Furthermore, we assume that the reactions happen on a 2D surface, such as the plasma membrane of a cell. Hence, all of our BD-GFRD simulations were done on a 1µm × 1µm plane with periodic boundary conditions. Using this setup we ran the simulations for 96 cpu-hours or 10 4 dimeriza-tion events, whichever occured first. We could reach a timescale of about 100 ms for 300/µm 2 , the largest concentration considered here.
While using spherical particles is necessary to use current BD-GFRD framework, the interaction potential need not be isotropic. However, for the ease of exposition we have used isotropic potential in this paper. We have studied reaction kinetics under both attractive and repulsive potential. We have modeled the former using a LJ potential with cutoff at 2.5σ and the latter with a WCA interaction [53] (Fig. 1C). In this paper we report results for cases when the interaction strength for both potentials were 6k B T . However, we have checked our results for various interaction strengths up to 10k B T and have found that variation of interaction energies within this range does not affect our results. The lengthscale of the interaction potentials, σ, was chosen to be equal to the sum of the radius of the two reactants.
We consider two different types of association reactions: homodimerization and heterodimerization. For the former, we consider reactions of type A + A → A 2 , while for the latter, we consider reactions of the form A + B → AB. To simulate association reaction we used a special case of Doi's volume reaction model [30,31]. We assume that the particles react as soon as their distance becomes less than or equal to the reaction distance r a (Fig. 1D). The reaction radius can be larger or smaller than r min , the location of the minimum of the interaction potential. Based on our previous work, unless otherwise stated, we assume that r a = 1.2r min . We have found that the concentration dependence does not depend qualitatively on r a , but it has quantitative impact, which we have discussed in this paper.

Measurement of κ(Φ)
We measured the reaction rates in an ensemble, in which, the total number of particles of all species remained constant. We started with all monomers and no homo-or heterodimers. Hence, we removed a product molecule from the simulation as soon as it formed and replaced it with the reactant molecules that it was formed from. We placed the reactant molecules randomly on the simulation box to avoid introduction of unwanted correlation. Because we consider only dilute systems, such random replacements very rarely encounter another molecule after the replacement and, hence, does not break detailed balance. Furthermore, if there is such an encounter, we reject the random replacement and repeat the random placement until there is no encounter. In fact, this process is equivalent to starting with a new initial condition after each dimerization event. For denser systems, which we have not considered in this paper, we have to use sampling techniques, such as the metropolis algorithm, that preserve detailed balance.
To measure the reaction rates, we computed the probability distribution function (PDF) of the time interval between two product formation reactions. The reaction propensities, r, were then calculated by fitting the exponential tail of the PDFs. The exponential tail results from the reactions, whose rates are lower than the diffusive encounter rate. Therefore, we can use well-mixed approximation, i.e., Eq. 7 to calculate κ(Φ).

Chemical kinetic model
We use κ(Φ) to construct chemical kinetic models that were solved using Gillespie algorithm or analytical calculation. For the dimerization reactions, each stochastic simulation were run until the simulation time reached t max = 100s or 100/k of f , whichever was shorter. The steady state values were measured by averaging the observables between t max /2 and t max over 100 different replicates. The system reached steady state within t max /10 for all parameter values explored. The chemical kinetic equation for the product concentration, P , is given by: In this equation, k on = r(Φ) = κ(Φ)Φ is the propensity of the concentration dependent diffusion-influenced reaction and Φ is the mass action of the reaction. For LMA, k on is Φ-independent constant, whereas for LCDR, k on depends on Φ. k of f is the dissociation rate. For reversible reactions approaching chemical equilibrium, k of f is given by the product of the association rate k on and the equilibrium constant k D , which ensures detailed balence. For nonequilibrium steady states, k of f were varied independently from k on .
Appendix C: Steady state concentration of molecules in dimerization reaction

Equation for steady state concentration
Let's consider the following dimerization reaction.
The concentrations of the molecules obey the following set of kinetic equations.
In steady state the following relation holds: In addition, we assume that the total number of monomers, N , is fixed, such that: Combining Eqns. C7 and C8, we get the following equation, whose roots are the steady state concentrations of the monomers.
where x = [A] is the steady state concentration of the monomer and Φ = x(x − 1)/2 is the mass action.

Solution
The solution of Eq. C9 depends on the nature of the steady state. For example, if the steady state is at chemical equilibrium, then k of f has to be equal to k on k D at all times, where k D is the equilibrium constant. As a result, for both LMA and LCDR, Eq. C9 simplifies to which has a simple solution: However, when the steady state is a nonequilibrium steady state, k of f does not have to be exactly equal to k on k D ; we can control it independently from k on . For the sake of simplicity let's assume that k of f = k D . Furthermore, we assume that under LMA k on = k 0 , where k 0 is a constant, from which we find that the solution at the nonequilibrium steady state is: where k sc D is equal to k D /k 0 . It is evident from this solution that, under LMA, it is not possible to attain a nonequilibrium steady state by manipulating k of f . Doing so, merely scales k D to a different value and the steady state fluctuations are equilibrium fluctuations that obey detailed balance.
In contrast, when LCDR is used it is possible to break detailed balance by manipulating k of f . To show this we compute the steady state solution using k on = κ(Φ) and k of f = k D . Unfortunately, closed form solution is not easy to obtain. So, we use Newton-Raphson root finding algorithm to find the solutions of Eq. C9. Steady state solutions for a few Φ 0 and N are shown in Fig. 9. Clearly, the nonequilibrium steady state solutions are different for LCDR than the equilibrium solution that we obtain for both LMA and LCDR. In general, the higher the Φ 0 , the more similar are the LCDR and LMA solutions. Physically, it implies that the validity of LMA increases as the difficulty for an association reaction increases. If the association reaction is a rare event, then the reaction events will follow Poisson distribution. Hence, the reaction rate will be a concentration independent constant.
The Jacobian is: which in steady state has the simple form The eigenvalues of this Jacobian satisfy the following characteristic equation.
When α > γ, the real part of the eigenvalues are negative and the steady state is a stable fixed point. On the other hand, when γ > α, the real part of the eigenvalues is greater than zero and the steady state is unstable to perturbations and the predator population grows exponentially and the prey populations shrinks exponentially, eventually resulting in extinction of the both prey and predator population. However, one should note that true collapse of the population is not possible in ODE based model, because x = 0, y = 0 is a saddle point, which means that the population can recover as long as it is greater than 0.
When α = γ, the population may reach its steady state through oscillation or without oscillation. The crossover between these two regions happen when the imaginary part of the eigenvalues become 0. That is when, The crossover region is defined by a line with slope m = η±1 η∓1 . However, we note that δ increases with Φ. Hence, m changes with Φ. At large enough Φ, m ≈ 1, whereas when Φ, hence δ, is small, the slope m → ∞ or → 0, depending on whether α < γ or α > γ, respectively. We can combine the preceding analysis to find the steady state solution of the oscillatory state (Fig. 10) and construct a phase diagram, which is shown in Fig. 11 Figure 11. Phase diagram obtained using the perturbation theory. Perturbation theory correctly predicts the phase boundary between stable fixed points, unstable fixed points, and sustained oscillation (α = γ, yellow line). However, it incorrectly predicts the phase boundaries that distinguish oscillating from non-oscillating approach to the fixed points (broken blue and crimson lines) at small α and γ. A key drawback of the perturbation theory is that it does not account for the sigmoidal nature of κ(Φ) and its saturation at κmax, which removes the non-oscillating regions (Explosion and Collapse & Fixed point) from the phase diagram.

Appendix E: Computation of κmax using Smoluchowski theory
The concentration dependent reaction rate in the presence of collisional-interaction is given by: where D is the diffusion coefficient, σ is the sum of the radius of two reacting molecules, and b is the average radius of a circular region where only one reaction is possible [18]. Clearly b is a concentration dependent quantity and the higher the concentration of the reactants the lower the value of b. In fact, we can estimate b as a function of concentration using a meanfield approximation. If the mass action is Φ and the area of the simulation box is A, then on an average, the area per reaction is A/Φ and b is given by [18]: where r 1 and r 2 are the radius of the two interacting molecules. Plugging this expression in Eq. E1 for A 2 , we get figure 12, which shows that κ max ≈ 20 and is approximately 0.4 for small Φ values. Furthermore, κ 0 can be estimated by extrapolating Eq. E1 to the Φ → 0 limit. Practically, it is estimated by finding the value of κ/D at Φ = 1, which is approximately 1.2 (Fig. 12A). The value of the exponent derived from the meanfield theory is remarkably close to the values found through the analysis of the simulation data, but it is not exactly equal to what we found. This difference stems from the fact that A/Φ underestimates the area available for each reaction. In reality, area available varies as A/Φ α , where α < 1. We found that as α is lowered from 1, computed using Eq. E1 also decreases. Therefore, for a suitable α, the value of will be 0.375. However, we need careful investigation of the diffusion limited reactions to identify such an α.