Inﬂuence of neighboring reactive particles on diffusion-limited reactions

Competition between reactive species is commonplace in typical chemical reactions. Speciﬁcally the primary reaction between a substrate and its target enzyme may be altered when interactions with secondary species in the system are substantial. We explore this competition phenomenon for diffusion-limited reactions in the presence of neighboring particles through numerical solution of the diffusion equation. As a general model for globular proteins and small molecules, we consider spherical representations of the reactants and neighboring particles; these neighbors vary in local density, size, distribution, and relative distance from the primary target reaction, as well as their surface reactivity. Modulations of these model variables permit inquiry into the inﬂuence of excluded volume and competition on the primary reaction due to the presence of neighboring particles. We ﬁnd that the surface reactivity effect is long-ranged and a strong determinant of reaction kinetics, whereas the excluded volume effect is relatively short-ranged and less inﬂuential in comparison. As a consequence, the effect of the excluded volume is only modestly dependent on the neighbor distribution and is approximately additive; this additivity permits a linear approximation to the many-body effect on the reaction kinetics. In contrast, the surface reactivity effect is non-additive, and thus it may require higher-order approximations to describe the reaction kinetics. Our model study has broad implications in the general understanding of competition and local crowding on diffusion-limited chemical reactions. © 2013 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution 3.0


I. INTRODUCTION
In diffusion-limited reactions, chemical species diffuse and almost immediately react when within a given radius of the reaction site. Many reactions in solution, such as fluorescence quenching, radical recombination, and growth of colloidal particles, occur through diffusion-limited reaction and because of this prevalence, it has been studied intensively and extensively. 1 Although in some cases, diffusion-limited reactions can be described by a simple diffusion equation, complicated geometries, boundary conditions, and external forces resist analytical solution. Competition between reactive sites may further complicate estimation of reaction kinetics. This difficulty originates from many-body interactions between a reactant and its environment (other reactive species), and indeed the many-body problem has been quite challenging even for simple homogeneous bimolecular reactions in the high concentration limit. [2][3][4] In a bimolecular reaction (A + B → products), calculation of reaction-related quantities like the rate coefficient is non-trivial in the presence of a large number of A and B reactants. For example, to consider all the reaction events in a system, we may consider an ensemble and its time-evolutions describing all the possible reaction pairs between A and B in timely order. Then for each time-evolution, in principle, a) C. Eun and P. M. Kekenes-Huskey contributed equally to this work. b) Authors to whom correspondence should be addressed. Electronic addresses: ceun@ucsd.edu and pkekeneshuskey@ucsd.edu.
we can calculate the contributions of each A-B pair to the rate coefficient and by averaging the individual rates, we can obtain the ensemble-averaged rate coefficient. However, the computational demands of this global approach scale exponentially with increasing numbers of A and B reactants, and the calculation quickly becomes intractable. Alternatively, from a local viewpoint focusing one reaction event, we may simply consider the reaction of a certain A-B pair, but in presence of other reactants, we also have to consider other pairings between A or B and other reactants. Hence, this many-body effect forms the foundation of kinetic competition. Microscopic kinetic theory avoids this complication by focusing on the dilute cases, in which one reactant is in vast excess over the other reactant. 5 The dilute condition can reduce the many-body problem (many A:many B) into a twobody problem, with the ignorance of the interaction between the excess species (many A:one B). The competition effect in diffusion-limited reactions arises when the concentration of reactants is sufficiently high, or the reactants are concentrated in a certain confined region, at which points that interference between the same species becomes non-negligible. 4,[6][7][8][9][10] Examples of the latter situation include reactants with multiple reaction sites [11][12][13][14] and a cell surface with multiple receptors for ligands. 15,16 Considerable effort has been devoted to explain the competition effect in diffusion-limited reactions; general formalisms have been developed to incorporate the competition effect into diffusion-limited reaction theories 3, 4 but because of intrinsic complication from the many-body effect, approximations such as monopole approximation were used to make the theory calculable. 6,17 Without approximations, exactly solvable models are limited to a few cases; for example, diffusion-influenced reaction with multiple point-like static traps (sinks), 18 and diffusion-limited reaction with two identical spherical reactants. 7 These exact solutions are needed to evaluate how accurate the approximations used in the theories are, or to develop another theory based on them (e.g., many-body theory based on an exact pair correlation function).
The two-sphere model has been widely used for testing theories; 7,17,19,20 due to its simple geometry, it is applicable to many situations with the consideration of the excluded volume and size effects of particles. Beyond serving as a test model, the original two-sphere model 7 by Samson and Deutch was generalized to incorporate other features such as finite reactivity and different size, and to include multiple spheres to consider multiple reactants, for application to various situations. 9, 10, 21-24 However, in most analytical works, only two-sphere cases were mainly considered in the calculations of physical quantities such as rate constant. This is because of difficulty in solving differential equations with more than three spheres (or with many local boundary conditions). Consideration of multiple spheres is needed for understanding the surrounding effect on kinetics, and probing the influence of the microscopic environment (e.g., molecular structure surrounding an active site in an enzyme).
The goal of our research is to understand how the local environment surrounding a reactant influences the kinetics of the reactant using multi-body models, in which the reactant and its local environment are represented by spheres. Precisely speaking, it is known that in the absence of any neighboring particles, the diffusion-limited rate coefficient is simply k = 4πσ 1 D (σ 1 : effective radius of reactant, D: relative diffusion constant), the well-known Smoluchowski rate coefficient, 1 and we focus on how this rate is changed in the presence of other neighboring particles around the reactant. Specifically, we consider varying numbers of spheres at a range of distances from the reactant sphere, as well as changing the reactivity of the surrounding spheres. These effects are investigated by numerical solution of the differential equation.
The outline of our paper is as follows. In Sec. II, we explain our reaction model, the associated mathematical expressions, and our numerical calculation method. In Sec. III, we compare our numerical calculations with the known analytical solutions in the two-sphere models, to evaluate how accurate our calculations are. Then in the two-sphere models, we further investigate the excluded volume effect as well as the size and the surface reactivity effects. Beyond the two-sphere model, we study the three-, four-, and five-sphere models to understand how the surrounding spheres influence the kinetics of the reactant sphere. Additionally, these models are used to address a local many-body effect associated with the angular distribution and local concentration of the surrounding spheres. In Sec. IV, we summarize our findings with the concluding remarks.

II. MODELS, THEORY, AND NUMERICAL CALCULATION
In this work, we study bimolecular reactions in the presence of neighboring particles centered around one reactant, to understand their influence on reaction kinetics. In particular, we consider diffusion-limited reactions of the following form: Here, the products can be one or many species, depending on the reaction type. In general chemical reactions, the reactants A and B can produce many different chemical species, and in protein-protein aggregation, the reactants produce only one product (the bound complex). For enzyme reactions, the reactant A or B may appear as products. These scenarios are all dependent on the initial reaction event between A and B, which is the subject of this study. To simplify the diffusion model, we make several assumptions: (1) the reactant B moves much slower than the reactant A and thus, we fix the position of reactant B in our coordinate system (target problem); (2) the concentration of the reactant B is very low, so that we consider only one reactant B; (3) the A reactants are non-interacting to each other; (4) the neighboring particles N around the B reactant are at fixed locations; and (5) the particles and reactants are represented as spheres. These assumptions permit tractable evaluation of a variety of factors that influence the kinetics on a representative system in Fig. 1. For convenience, we use a mathematically equivalent model with the effective radii, σ 1 (= R A + R B ) and σ 2 (= R A + R N ), where the reactant A can be considered as a point particle.
In the model systems, the diffusion-limited reaction can be described by a simple diffusion equation with absorbing boundary condition at the surface of the target. In this initial study, we only consider the time-independent problem or the steady-state case. Since we are interested in steady state, the time-independent equation we have to solve is given below. In our target problem, a reactant B is fixed at a certain point in the coordinate system and the A reactants are diffusing around the reactant B, where D is diffusion constant of the reactant A and c A ( r) is the concentration of the reactant A at the position of r. This differential equation has two types of boundary conditions; one is when the reactant A is far away from the reactant B (outer boundary condition), and the other one is when it is in the local environment around the reactant B (inner boundary condition). Particularly, the inner boundary condition is used for describing reactions with the reactant B and the neighboring particles (spheres). In our model, the outer boundary is lim | r|→∞ c A ( r) = c ∞ but the inner boundary condition is different depending on the local environment provided by the neighboring spheres around the reactant B. In the absence of neighboring spheres, the inner boundary condition is simply expressed by one equation, c A ( r)| | r|=σ 1 = 0 representing full reactivity (fully absorbing boundary condition). However, in the presence of neighboring spheres, it is difficult to express the inner boundary condition with a single equation, as is in the absence case. The only exceptional case is the two-sphere model with bispherical coordinates. 7,25 Thus, in general, to represent inner boundary conditions, it is convenient to introduce local coordinates whose origin is the center of each sphere. In this case, similarly to the case above, the boundary condition for each sphere i is represented by either for finite reactivity (partially reflecting boundary condition), wheren i ( r i ) is the normal vector from the reactive surface at r i and κ i is the inherent rate constant characterizing the reactivity of sphere i. Once Eq. (1) is solved for c A ( r) in a given set of boundary conditions, we can calculate the rate constant R for the reactant B through R = D∇c A ( r) · d S B , which is the integral over the surface area of the reactant B. However, for our convenience, we use the normalized rate constant or rate coefficient k, defined by k = R/c ∞ . In the absence of any neighboring particles, k becomes k 0 = 4πσ 1 D, the Smoluchowski rate coefficient. 1 Along with the rate coefficient, to get scale-independent results, we report our results using dimensionless variables u A and λ i (inverse reactivity), 9 as shown in the below, and dimensionless units such as d/σ (d: separation distance, σ : radius of sphere): In our study, we use D = 780 μm 2 s −1 from Ca 2+ diffusion in solution, 26 although the actual value is arbitrary since we consider steady state and dimensionless units. All other parameters are discussed in the text. Although this is a well-defined mathematical problem, in practice, it is very challenging to solve the differential equations analytically in the presence of the neighboring particles. The presence of particles makes it difficult to describe their boundary conditions in a single global coordinate frame. The only exceptional case is that there exists only one neighboring sphere around the reactant sphere. In this case, the bispherical coordinates can be employed to simply describe the boundary condition and then c A ( r) can be obtained in terms of bispherical coordinates. 7 However, in general cases, where more than two neighboring spheres exist, there are no such global coordinates for simply expressing multiple isolated local boundary conditions as one equation. Thus, instead of using a single coordinate frame, one may use as many sets of local coordinates as the boundary conditions to describe each boundary condition. Then, using mathematical relationships between the local coordinates, in principle one can solve the differential equation. This approach has been applied to twosphere models to study the size and reactivity effects, 9,22 and the differential equations were solved by using a shift formula or an addition theorem of Legendre functions 27 to relate two boundary conditions expressed in different local coordinates. However, as more spheres are introduced, more local coordinates and more relations between the coordinates are needed to solve the differential equations, which makes the problem insolvable.
To avoid the mathematical difficulty encountered in analytical approach, we solve the three-dimensional diffusion equation via the finite element method, 28 as was done previously for a single enzyme binding description. 26 We extend the python-based program Smolfin 26 to numerically solve the differential equations of the problems described above. The molecular mesh for defining spatial domain was constructed by GAMER plug-in Ref. 29 for varying number of neighbors, radii, and particle distances, and according to the boundary conditions, the boundary regions in the meshes are marked.
Smolfin reads the marked finite-element mesh and defines appropriate boundary conditions according to the marks. Then Smolfin solves the diffusion equation via the finite element method, based on the weak form representation of the partial differential equation (PDE) given by ∇u A ( r) · ∇v( r)d r = 0, where v is a test function in a piecewise linear function space over the mesh. The solution, u A , to the weak form of the PDE is obtained by using the default direct linear solver in FEniCS solver (http://fenicsproject.org). 30 Smolfin estimates the reaction rate R = D∇c A ( r) · d S B . Additional details and references on the Smolfin solver are given in Sec. A 3 of Ref. 26.

III. RESULTS AND DISCUSSION
Our primary goal is to understand how neighboring particles localized around a reactant B influence the kinetics of the reaction of the reactant B with the reactant A. To achieve this goal, we have used sphere models that are simple yet capture the essential physics. In the sphere models, the neighboring particles (spheres) are introduced and they are characterized by geometrical and the chemical properties; geometrical properties include the neighboring particles' number, size, and distribution about the reactant B, while the chemical property describes the reactivity of the neighboring particles. These models can be applicable for general situations. For instance, if the neighboring particles have reactivity with the reactant A, they can compete with the reactant B for the same reactant A, but if they have zero reactivity or they are neutral for the reaction, they can be simply considered as crowders by increasing the excluded volume, or the volume which is not available for reactant A.

A. Comparison with the analytically solvable models
Before we systematically study the effect of the neighbors with the aforementioned factors on the kinetics, we first examine the accuracy of our numerical solutions by comparing with the known analytical solutions. For this comparison, we use the two identical sphere model studied by Samson and Deutch 7 and a more generalized two-sphere model by Mc-Donald and Strieder. 9 In the Samson-Deutch model, the reactant B has one neighbor that is identical to the reactant B in all aspects, and both of them are described by fully absorbing spheres. Samson and Deutch obtained the exact mathematical expression of the rate coefficient k using bispherical coordinates, where Here, k 0 is the rate coefficient when neighboring particles are absent, d is distance between two spheres and σ = σ 1 = σ 2 is the radius of the spheres.
With the analytical solution, we perform the comparison with our numerical result in small and large separation regions. The comparison result is shown in Fig. 2. At small separations, Fig. 2(a) shows excellent agreement between the exact and the numerical solutions. The overall error is less than 1%. As is expected, the normalized rate coefficient increases with the separation between the two spheres. This clearly can be explained by the competition effect induced by a neighbor. When the neighbor is close to the target reactant B, the concentration of the reactant A around the reactant B is reduced because of the reaction at the neighbor and this leads to the decrease in the rate coefficient for the reactant B. However, when the neighbor is moving away from the reactant B, the competition effect is decreasing and this is why the rate coefficient increases. We also calculate the rate coefficient at large separations; as shown in Fig. 2(b), the normalized rate coefficient approaches 1.0, the Smoluchowski limit. However, the coefficient is slightly overestimated, compared to the analytical solution, but by less than 2% which could be attributed to numerical or meshing error.
We also consider a generalized two-sphere model proposed by McDonald and Strieder, in which the ratio between the radii of two spheres and the reactivity of each sphere are given as model variables. They obtained the analytical expressions for the concentration of the reactant A and the rate coefficient in recursive forms, which we use to validate our calculations. The exact expressions can be found in Ref. 9. Figure 3 shows the comparison of the rate coefficient for the reactant B in the three cases discussed in the work of McDonald and Strieder. The cases in Fig. 3(a) are that the reactant B is 10 times smaller than the neighbor, and the reactivity of the reactant B is changed from large to small (λ 1 = 0.02 − 50), while the reactivity of the neighbor is fixed as a relatively large constant (λ 2 = 0.02). Note that small value of λ implies large reactivity. The cases in Fig. 3(b) are the same with the ones in Fig. 3(a) except that the reactivity of the neighbor is now relatively small (λ 2 = 50), whereas the cases in Fig. 3(c) are the same with the ones in Fig. 3(a) except that the size ratio between two spheres is reversed (B:neighbor = 10:1).
For the three cases, Fig. 3 demonstrates the comparison between the analytical and the numerical results. In Figs. 3(a)  and 3(b), the results are in good agreement, except at contact and very close separations (d/(σ 1 + σ 2 ) = 1.0-1.1), where obtaining reasonable finite element meshes is difficult. The averaged relative error to the exact value over d/(σ 1 + σ 2 ) = 1.2-6.0 is also calculated. In Fig. 3(a), the errors are 9%, 8%, 4%, 1%, 3%, 5%, and 6% for λ 1 = 0.02, 0.1, 0.5, 1, 5, 10, and 50, respectively. Likewise, in the same order of λ 1 , the errors in Fig. 3(b) are < 1%, 1%, 5%, 7%, 12%, 13% and 14% and the errors in Fig. 3(c) are < 1%, 1%, 6%, 10%, 16%, 17%, and 18%. As we can see here, there is a general trend that as the reactivity decreases (λ 1 increases), the error becomes larger. Since we consider diffusion-limited reaction with λ 1 = 0 for the reactant B, we expect that the maximum error in non-contacting regions (d/(σ 1 + σ 2 ) ≥ 1.2) would be ∼10% from the errors in Fig. 3(a). However, when two spheres are very close, the error is significant in Fig. 3(b) as well as Fig. 3(a). The error is probably from the difficulty in creating a good quality of volumetric mesh around two asymmetric spheres at contact or near contact and the numerical errors associated with the discretization for largely changing variables. In fact, the cases of Figs. 3(a) and 3(b) are extreme for our study, in that the neighbor is much larger than the reactant B. In our study, we mainly use the same size of sphere for the reactant B and neighbors, except the study for the size effect in Sec. III B. Indeed, when the size of the neighbor is comparable to, or less than the size of the reactant B, the error is 1% or so in the small separation region including the contact region, as is seen in Fig. 3(c) as well as Fig. 2. Thus, the worst cases that we have to be careful about in numerical calculation are when relatively large neighbors are very close to the reactant B. However, even in such cases, our result is qualitatively correct, in that the rate coefficients with different λ 1 s are arranged in the same order in the both analytical and numerical results. Thus, our numerical method is useful to understand general behaviors in the influence of neighbors.

B. Effects of reactivity and relative size in the two-sphere model
The McDonald-Strieder model suggests that the normalized rate coefficient (k/k 0 ) for the reactant B is dependent on the reactivity and the relative size of a neighbor. For the reactivity effect, the two rate coefficients shown in the top of Figs. 3(a) and 3(b) are quite different in that the one in Fig. 3(a) is slowly increasing with the separation, while the one in Fig. 3(b) is rapidly increasing, although two systems are exactly the same except the reactivity of the neighbor (λ 2 = 0.02 vs. λ 2 = 50). It is apparent that this difference comes from the relative strength of competitive reaction determined by the reactivity. Similarly, the comparison between the rate coefficients in Figs. 3(a) and 3(c) can explain the size effect due to the neighbor because the two systems are the same except the relative size (σ 1 :σ 2 = 1:10 vs. σ 1 :σ 2 = 1:0.1). When the neighbor is much smaller than the reactant B, the influence by the neighbor on the rate coefficient is insignificant, hence the coefficient is nearly constant in Fig. 3(c).
To provide greater insight into general trends in the size and reactivity effects, and the quantitative dependences of the rate coefficient on the relative size and the reactivity, we systematically change the reactivity and the size of the neighbor. In the study of reactivity, we simply employ the Samson-Deutch model and investigate how the rate coefficient is modified by changing λ 2 . Here, we include the limiting cases with λ 2 = 0 and λ 2 = ∞, to clarify the influence of the neighbor's reactivity. λ 2 = 0 and ∞ correspond to fully absorbing and fully reflecting boundary conditions, respectively, while the finite value of λ 2 implies partially reflecting boundary condition. For a given reactivity, the rate coefficient for the reactant B is calculated as a function of separation between the reactant B and the neighbor. The results are shown in Fig. 4(a). As the reactivity increases (λ 2 decreases), the rate coefficient decreases because of the competition reaction (see the change along the red line in Fig. 4(a)). Interestingly, when λ 2 = ∞, the rate coefficient is slightly less than 1, which implies that the reflecting sphere does not greatly influence the rate coefficient. 24 Nevertheless, in all the cases, the rate is less than 1, and the presence of a nearby spherical neighbor always decreases the rate coefficient for the reactant B, compared to the case of an isolated reactant B. This result clearly demonstrates that the reactivity plays an important role in the modification of the rate coefficient; depending on the strength of reactivity, it decreases the coefficient by up to 30%. Also, this result implies that the reactivity may be more influential than   FIG. 4. (a) Normalized rate coefficient as a function of dimensionless distance with the inverse reactivity of the neighboring sphere (λ 2 ). (b) Normalized rate coefficient as a function of the inverse reactivity of the neighboring sphere (λ 2 ) for distances (d/σ = 2.1, 2.3, 2.5 (shown in (a)), 2.7, 2.8). The result is obtained by taking the average over three data sets (σ = 2.5, 5, 7.5) and the error bar represents the standard deviation. Insets clockwise from the left: concentration contours about reactant A at d/σ = 2.5 for λ 2 = 0.01, λ 2 = 2, and λ 2 = 1, respectively. the separation distance; for example, compare the two cases (d = 2σ with λ 2 = ∞ and d = 10σ with λ 2 = 0). Additionally, from this result, we can obtain the relation between the rate coefficient and the inverse reactivity (λ 2 ). In Fig. 4(b), we plot the rate coefficient as a function of the inverse reactivity (λ 2 ). As we can see, the curve is nonlinear. The slope of the rate coefficient is relatively large in the interval of λ 2 from 0 to 1, but beyond 1, it is small. Thus, roughly speaking, when λ 2 > ∼1, the boundary condition of the neighbor is primarily reflective; when λ 2 < ∼1, the boundary condition has more absorbing character. Particularly, from Fig. 4(a), we notice that the curve with λ 2 = 0.01 is very close to the one with λ 2 = 0. Thus, when λ 2 < ∼0.01, the boundary condition is essentially fully absorbing. Finally, to understand the physical meaning of the change with λ 2 , we draw the contours of concentration profile of reactant A around two spheres, as shown in the insets of Fig. 4(b). When λ 2 = 0.01, the contours are symmetric but when λ 2 ≥ ∼1, the contours are asymmetric. We understand that in the interval of λ 2 between 0 and 1, the transition between symmetric and asymmetric contours should take place. Also, here, as we can expect, we can notice that the contours with λ 2 = 1 and λ 2 = 2 are very similar.
In addition to the straightforward study on the reactive effect, this study serves as another purpose for understanding the excluded volume effect of a neighbor, when λ 2 = ∞. In general, the nearby particles affect the kinetics of the reactant B through two ways, providing the excluded volume blocking the reactant A to access the reactant B through them, or consuming the reactant A at the particles by the competition reaction, which is discussed in Fig. 4. Since the former is dependent on the volume while the latter depends on the surface area with the reactivity, these give rise to two separate influences on the reaction rate: the volume effect and the surface effect, respectively.
To understand the volume effect, we use fully reflecting boundary condition (λ 2 = ∞) for a neighbor, which eliminates the surface reactivity effect from the neighbor. In this study, we consider the size of the neighbor as 10 times larger or 10 times smaller than the size of the reactant B. First, in the cases of large neighboring spheres, the Fig. 5(a) shows that as the neighbor is larger, the rate coefficient is significantly decreasing. However, for small neighboring spheres, the effect is negligible as the size decreases (see the inset of Fig.  5(a)). Physically, the neighbor increases the excluded volume to the reactant A. Thus, compared to the case without any diffusion barrier, this will reduce the concentration of the reactant A around the reactant B, and hence the flux into the reactant B. Consequently, it will lead to the decrease of the rate coefficient. This effect is more apparent for decreasing distances between the reactant B and the neighbor. Considering that the rate coefficient hardly varies as a function of distance when the reactant B and the reflecting neighbor are of the same size (see black in Fig. 5(a)), it is interesting that this trend changes significantly as the neighbor size increases. Although the results may not be quantitatively accurate in the contact region as we discussed in Sec. III A, it is a quite significant drop, even in the non-contact region (e.g., d/(σ 1 + σ 2 ) = 1.2). Hence, the excluded volume effect is especially prevalent when the neighbor is large. When the excluded volume and the surface effects are considered in tandem, the influence of the neighbor is much larger, as is shown in Fig. 5(b). To study the limiting case of the reactive neighbor, we use fully absorbing boundary condition (λ 2 = 0). Compared to the cases in Fig. 5(a), the coefficient is much smaller at the given distance, due to the competition reaction. However, as the size of neighbor decreases, even with full reactivity, the decrease of the coefficient could be negligible, as is shown in the case of σ 1 :σ 2 = 1:0.1 (see the inset of Fig. 5(b)).
From the comparison between Figs. 5(a) and 5(b), it is apparent that the excluded volume effect alone is relatively short-ranged. Beyond d/(σ 1 + σ 2 ) = ∼2, the coefficient is nearly unity and practically independent of separation. However, when the surface reactivity is considered, the radius of influence of the neighboring particle is strongly extended. The critical distance in the transition between short-and longranged behavior is around ∼2(σ 1 + σ 2 ).

C. Dependence of the angular distribution of the surrounding particles in the three-sphere model
Next, we extend the two-sphere model to multiple-sphere models, where the reactant B is surrounded by more than one sphere. When we consider multiple neighbors surrounding the reactant B, the kinetics are dependent on their spatial distribution and the local concentration around the reactant B. We study these aspects in this section and in Sec III D. Specifically, we investigate two limiting cases for the boundary condition of the surrounding neighbors: fully absorbing and fully reflecting boundary conditions. The result is obtained by taking the average over three data sets (σ = 2.5, 5, 7.5) and the error bar represents the standard deviation.
As the simplest case for multiple-particle models, we use a three-sphere model. Our three-sphere model represents the scenarios in which two neighbors are equally distant from the reactant B but the distance between the two neighbors varies. This case probes how the rate coefficient is dependent on the spatial distribution of the surrounding particles around the reactant B. Figure 6 shows the rate coefficient change as a function of separation in the given angular distribution of the neighbors around the reactant B. As we can see in Fig. 6(a), with reflecting boundary condition of the neighbors, the coefficient does not depend on the angle. Within the error bars, the values of the coefficient with different angles are not distinguishable, which means the relative angular distribution is not important to the kinetics. However, with fully absorbing boundary condition, we clearly see a trend that at a given separation, the rate coefficient increases as the angle becomes smaller. When the angle is 180 • , or the two surrounding particles are maximally separated to each other, the rate coefficient is lowest. That is, at the largest angle, the competition reaction due to the neighbors is maximized. However, when the two surrounding neighbors are in contact to each other, the rate coefficient is highest. This can be explained by the weakening of the competition effect; as the surrounding neighbors approach one another, they are competing for the reactant A, which lessens their competition with reaction B (thus yielding a larger effective rate). In the case where two neighbors are directly adjacent, the two neighbors can be considered one neighbor with a smaller effective reactive surface area. To summarize, for chemically neutral particles, their angular distribution around the reactant B is less important in kinetics but for chemically active particles, it is important.

D. Many-body effect in the multiple-sphere models
In a physical system, reactant B may be surrounded by more than one particle and the kinetics may be highly dependent on the local concentration of the surrounding particles. To study this subject, we add more spheres to the two-sphere model used in Sec. III B and investigate how the rate coefficient changes with the number of spheres. To simplify the models, all the spheres are at the same distance from the reactant B and arranged in such a way to maximize their separation and hence their influence on the kinetics of reactant B as discussed in Sec. III C. That is, we use line geometry (180 • ) for two neighbors, triangle for three neighbors, and tetrahedron for four neighbors, where the center of the surrounding neighbors coincides with the center of reactant B. We consider two limiting cases for the reactivity in Fig. 7: fully absorbing and fully reflecting boundary conditions. For both cases, the rate coefficient decreases with the number of surrounding neighbors; consistent with the previous results from the two-sphere model, for the fully reflecting condition in Fig. 7(a), the decreasing effect is relatively small, yet is large for the fully absorbing boundary condition in Fig. 7(b), due to the surface reactive effect. We also observe that the effect becomes longer ranged as the number of the surrounding neighbors increases. This suggests that with more neighbors, the concentration profile of the reactant A is more perturbed and the influence can affect the profile around the reactant B at larger distances. Our next question is whether this apparent many-body effect on the reaction kinetics due to the surrounding neighbors is cooperative or synergetic in nature. For this, we consider a general perturbation formalism (similar to the density expansion of equation of state 31 ), described by k n (r) /k 0 = 1 − α (r) · n + β (r) · n 2 − γ (r) · n 3 + · · · , (5) where n is the number of the surrounding neighbors and r is the distance between the surrounding neighbors and the reactant B, and α(r), β(r), γ (r), . . . are the coefficients. Here, the coefficient of n i term is a function of r and it is related to the interaction among i neighbors. Thus, if there is no cooperative effect between the surrounding neighbors, all the coefficients except α(r) would be zero and Eq. (5) is simply reduced to the following equation: In this case, the normalized rate coefficient is linear with the number of neighbors. Moreover, this linear effect would allow reduction of the many-body problem into a two-body problem (reactant B and one neighbor). Precisely speaking, once α(r) in k 1 (r)/k 0 = 1 − α(r) is determined, k n (r)/k 0 can be calculated for any n-neighbor case. To distinguish this linear case in Eq. (6) from the general case in Eq. (5), we can define k linear n (r) /k 0 :≡ 1 − α (r) · n and we calculate it when n = 1, 2, 3, 4 for the cases in Fig. 7.
Since we have four data sets from the two-, three-, four-, and five-sphere models, α(r) can also be determined in the four ways under the linear relationship given by Eq. (6), α i (r) = (1 − k i (r)/k 0 )/i for i = 1, 2, 3, 4. In the actual calculation, we use the average value, α(r) = (α 1 (r) + α 2 (r) + α 3 (r) + α 4 (r))/4 to calculate k linear n (r) /k 0 for n = 1, 2, 3, 4 and we plot k linear n (r) /k 0 as a function of r in Fig. 8. Interestingly, the cases with reflecting boundary condition in Fig. 8(a) demonstrate excellent agreement between k n (r)/k 0 and k linear n (r) /k 0 , especially in small separation region. It means that the decrease in the rate coefficient is simply proportional to the number of the reflecting neighbors. In fact, this excellent agreement is also reflected in the calculation of α(r) . As we can see in the inset of Fig. 8(a), practically there is no difference between α i (r)s obtained from the different cases. It seems that α i (r) does not depend on the number of surrounding neighbors. Thus, the standard deviation (error bars in the figure) is relatively small, and so any choice of α i (r) for the calculation of the rate coefficient would give more or less the same result.
We also consider the cases with absorbing boundary condition. The results are shown in Fig. 8(b). In contrast to the cases with reflecting boundary condition, the inset of Fig. 8(b) exhibits that α i (r) shows noticeable changes with the number of the absorbing neighbors. As i increases, α i (r) decreases and thus, the average value α(r) is positioned in the middle between α 1 (r) and α 4 (r). For a reference, we also plot α DFS (r :≡ d/σ ) = r/(1 + r) = σ /(σ + d) predicted from the simple approximation for the two absorbing sphere model in the Deutch-Felderhof-Saxton theory, k 2 (r :≡ d/σ )/k 0 ≈ (1 + 1/r) −1 = (1 + σ /d) −1 . 6, 7 α DFS (r) seems to be an upper bound for all α i (r). Because of the dependence on i, the pre-  Fig. 7(a)). (b) Comparison of the normalized rate coefficients from the numerical calculation (solid line) and from the linear approximation (dotted line) for the cases of multiple absorbing spheres ( Fig. 7(b)). The inset shows α i for all the cases and their average, which is used for the approximation. The colors of the curves in (b) have the same meaning in the ones of (a). dicted rate coefficients show discrepancy from the numerically calculated ones, as is shown in Fig. 8(b). This implies that in the cases of many absorbing neighbors, the rate is not simply calculated from a linear term α(r) and it is necessary to consider higher-order terms such as β(r) and γ (r). In this case, the influence is so significant that not only do the surrounding neighbors influence the concentration profile of reactant A around the reactant B, but they also influence the profile around themselves (competition effect between themselves). However, at large separations (d = 9 ∼ 10σ ), the linear approximation works well, as we can expect from the convergence of α i (r) in the large separation region. This is because the distance between the neighbors is so large that their interference can be ignored.
Why does the linear approximation work well for reflecting boundary condition while it does not for absorbing boundary condition? As we discussed in Sec. III B, when a neighbor of the same size with the reactant B is fully reflective, it has only the excluded volume effect and the effect is not significant. It decreases the coefficient maximally at the contact by ∼5%. In fact, this small reduction explains why the linear approximation shows good agreement. That is, the perturbation of the concentration field by the neighbor is so weak that the effect is additive. Thus, the overall influence by the reflecting neighbors is the sum of the influences of individual neighbors and we can ignore higher-order correlation terms in the consideration. This also explains why there is no angle dependency in the three-sphere model in Sec. III C. The perturbation is again so weak that practically there is no interference between the surrounding neighbors. The insets in Fig. 7(a) evidence this weak perturbation, as the concentration contours (projected on the plane defined by the reactant and neighbors) are nearly circular and resemble the contours in the neighbor-free reaction. In contrast, for the absorbing boundary condition, the perturbation is large and thus we cannot ignore the higher order terms representing the interaction between neighbors, which gives rise to an angle dependence in the three-sphere model, as well as nonlinearity in multiplesphere models.
Specifically, as a result of neighbor-neighbor interactions, the concentration profile between absorbing neighbors is significantly perturbed and large depletion zones arise. For example, the insets in Fig. 7(b) show a dumbbell-shaped concentration profile for the two-neighbor case and a triangleshaped profile for the three-neighbor case. The depletion zone due to the many-body effect is relatively long-ranged and influences the kinetics in a non-additive way (e.g., cases at d = 5σ in Fig. 8(b)), even when the neighbors are largely separated.

IV. SUMMARY AND CONCLUSION
Using a spherical model representative of globular proteins and small molecules, we have studied the kinetics of bimolecular diffusion-limited reaction in the presence of reactive neighboring particles. The influence of the neighbors on the reaction is dependent on their distribution about the reactant, their relative size, their reactivity, and their local density. We systematically investigated these factors separately to isolate their contribution to the rate coefficient. In the two-sphere model, we studied the excluded volume effect and the reactive surface effect (competition effect), separately. If the volume of neighboring particle is comparable to or less than the volume of the reactant, the excluded volume effect is negligible; however, the reactivity effect can substantially influence the kinetics, particularly for large reactivity values (when the inverse reactivity λ 2 is less than 1) and as such, high reactivity could be more important than close proximity in the kinetics. This is due to the fact that the reflecting boundary condition (non-reactive case) minimally perturbs the concentration of the diffusing reactant adjacent to the target reactant, in contrast to the long-range influence of the absorbing boundary condition (reactive case). Nevertheless, the excluded volume effect could be significant when the volume of neighbor is much larger than the volume of reactant.
Although the excluded volume effect is relatively shortranged for a single reflecting neighbor case, it could be longer ranged when more reflecting neighbors are added to the space near the reactant. This is because as the concentration increases, its influence on kinetics spans larger distances, even with reflecting neighbors. Interestingly, this effect is linear with the number of reflecting neighbors, because the perturbation is relatively weak. By the same rationale, it is independent of the neighboring particles' angular distribution about the reactant in the three-sphere model.
When the neighbors are fully absorbing, the perturbation of the concentration profile around the reactant is quite significant. In this case, we observed the angle dependence of the kinetics in the three-sphere model stemming from the separation between neighbors. The neighbors' influence is minimized when they are nearly in contact, while it is maximized when they are maximally separated (straddling the target reactant). At small inter-neighbor separations, the neighbors compete with each other and relieve their competition with the target reactant; at maximal neighbor separations, this competition between neighbors is minimized. The dependence of the reaction rate on both the neighbor-target and the neighborneighbor distances gives rise to a many-body relationship governing the chemical kinetics, for which linear approximations are not well-suited.
The simplicity of our model study permits a broad survey of factors influencing diffusion-limited reactions, which may have substantial implications for a large class of chemical reactions. Such reactions typically occur in heterogeneous environments, wherein the primary chemical species are surrounded by particles secondary to the reaction. These particles are suitably described as crowders when they do not directly participate in the primary reaction, and are best described as competitors when they exhibit significant reactivity with the reactants. Our study thus provides important insight into the dependence of chemical kinetics on the nature of surrounding particles in a heterogeneous system. Our numerical model furthermore constitutes a step toward examining the role of substrate charge in electrostatically steered reactions, [32][33][34] and is the subject of a follow-up study.