Diffusion-Limited Reaction Kinetics of a Reactant with Square Reactive Patches on a Plane

: We present a simple reaction model to study the influence of the size, number, and spatial arrangement of reactive patches on a reactant placed on a plane. Specifically, we consider a reactant whose surface has an N × N square grid structure, with each square cell (or patch) being chemically reactive or inert for partner reactant molecules approaching the cell via diffusion. We calculate the rate constant for various cases with different reactive N × N square patterns using the finite element method. For N = 2, 3, we determine the reaction kinetics of all possible reactive patterns in the absence and presence of periodic boundary conditions, and from the analysis, we find that the dependences of the kinetics on the size, number, and spatial arrangement are similar to those observed in reactive patches on a sphere. Furthermore, using square reactant models, we present a method to significantly increase the rate constant by sequentially breaking the patches into smaller patches and arranging them symmetrically. Interestingly, we find that a reactant with a symmetric patch distribution has a power–law relation between the rate constant and the number of reactive patches and show that this works well when the total reactive area is much less than the total surface area of the reactant. Since our N × N discrete models enable us to examine all possible reactive cases completely, they provide a solid understanding of the surface reaction kinetics, which would be helpful for understanding the fundamental aspects of competition between reactive patches arising in real applications.


Introduction
In a reaction between A and B, two reactant A molecules can compete with each other for a partner reactant B molecule, the outcome of which is strongly dependent on how close the competing A molecules are to each other [1,2].A similar competition also occurs between two reactive sites on a reactant surface [3][4][5] and can be generalized to cases of multiple sites [6,7].One specific example of the latter is a cell surface with multiple receptors competing for ligands [8].This case has been intensively studied through a simple model, in which the cell surface is modeled as a sphere with multiple circular patches [8][9][10][11][12][13][14][15][16][17][18][19][20][21], using analytical theories, such as the effective medium treatment [11], the homogenization theory [9,12,14,15], the matched asymptotic analysis [16], and the curvaturedependent kinetic theory [19], and numerical methods, such as Brownian dynamics simulations [10,12,13,17], the spectral boundary element method [18], and the finite element method [7].Moreover, recently, other related models have been proposed and studied using the matched asymptotic analysis [22,23].Despite these advances in the study of multiple reactive patches on a sphere, the study of the effect of the spatial patch arrangement on kinetics is mostly limited to a random, symmetric, or certain type of arrangement of patches [7].A study involving the complete examination of all possible arrangements has not been performed.
In a reactant with multiple reactive sites, any two reactive sites can kinetically compete with each other, depending on their distance, and, therefore, their overall arrangement on the surface can significantly affect the overall reaction kinetics [7].However, for continuum models such as reactive sites on infinite planes, since the sites can be placed at any location as long as the sites do not overlap, an infinite number of site arrangements exist, which increases the difficulty in investigating all possible arrangements.Even for a spherical reactant with a finite surface area, because of the space continuum, an infinite number of arrangements of reactive sites (patches) also exist, except for the cases in which the sites are maximally occupied on the surface such that no room exists for one to rearrange the sites, meaning that only one arrangement is available.Therefore, in this work, we focus on the development of a model that can generate a finite number of site arrangements, and we study the kinetics of that model.
To generate a finite number of arrangements, we consider the discretization of the surface of the reactant.One intuitive way to accomplish this is to create an N × N grid on the surface, with each cell defined by the grid lines, being either reactive or inert.When a cell is reactive, we can consider it a reactive patch.Furthermore, for our convenience, we assume that the shape of the reactant is a square, and consequently, a square with N × N grid generates N × N cells.Depending on which cell is reactive, the model can generate a variety of reactants; for example, when all cells are reactive, the reactant is fully reactive.We also assume that the reactant is on a plane, thus that the N × N cells are identical in shape.If the reactant is on a curved surface, the surface curvature must be considered, and the N × N cells will have different shapes.Because of the N × N square shape of the reactant, we call this model the N × N square reactant model.Another discrete model based on a triangular lattice or square lattice has been proposed, but that model study was focused on clusters of reactive N-disk patches [24]; in our model, we consider the cases of both separate patches on a reactant and clustered ones.
In this work, we first study the simplest cases of N = 2, 3 to see how the square reactant model works.As we did previously for a spherical reactant with multiple reactive patches [7], we study how the unit cell size, number, and arrangement of reactive patches affect the reaction kinetics.However, the variation in patch arrangement for N = 2, 3 is limited because of the small number of cells on the surface (4 and 9 cells for N = 2 and 3, respectively).To overcome this limitation and to obtain a general idea regarding the effect of the patch arrangements, we also study the extended model of N = 9 with some representative arrangements of reactive cells.Based on this understanding, we discuss how one can increase the overall rate constant for a given total reactive patch area.For the discussion, we propose a systematic way of increasing the rate constant by properly arranging reactive cells in a symmetric way.
This paper is organized as follows.In Section 2, we explain in detail our N × N square reactant models and computational methodologies.In Section 3, we perform a complete kinetic study of the 2 × 2 and 3 × 3 square models for two types of models (one with normal boundary conditions and one with periodic boundary conditions) and discuss the significance of the results.Using the square models, we also investigate the dependence of kinetics on the fraction of the total reactive patch area and the effect of the patch distribution.Based on these investigations, we propose a systematic way to increase the rate constant by breaking and symmetrically arranging reactive patches for a fixed surface coverage ratio of reactive patches.Furthermore, we discuss a very interesting power-law relationship between the rate constant and the number of patches.In Section 4, we conclude by summarizing our findings from this work.

Reaction Models
In this investigation, we used 2 types of reaction models to study a diffusion-limited (or diffusion-controlled) bimolecular reaction between reactants A and B. Here, reactant A was the one we were most interested in because it was a square reactant with a reactive patch distribution, while reactant B was simply a point-like particle.One type was for the situation where reactant A was localized on a certain chemically-inert plane.The other type was for the situation where reactant A had a large surface area or was enlarged for detained analysis; thus, its reactive patches were assumed to exist on the entire surface of the model.
The first type of model was used in our previous study [7,25].A schematic diagram of the reaction model is shown in Figure 1.Reactant A is placed at the origin of the coordinate system and is on a nonreactive x-y plane at z = 0.The surface of reactant A has a square grid with N × N cells, where N is the number of cells on each side.Each cell is reactive or inert for the reaction with the reactant B; when the cell is reactive, it is considered a reactive patch.Around reactant A (black square lines on the x-y plane in Figure 1), we assume that many point-like molecules of reactant B exist that do not interact with each other and freely move with a diffusion coefficient .When a molecule of reactant B contacts one of the reactive patches of reactant A, the reaction occurs immediately.Specifically, the reaction rate is solely determined by the transport rate at which B molecules diffuse to the reactive patches of A; this is why we called this reaction a diffusion-limited (or diffusioncontrolled) reaction.The reaction process above can be described by the Smoluchowski equation for the concentration of reactant B, which is basically a diffusion equation with boundary conditions describing reactions.Here, to further simplify, we normalize the concentration of reactant B by the bulk-solution concentration of B, which is the concentration far away from A, so that the concentration of B or is a dimensionless variable.We also assume that our model system is in a steady state.Consequently, the Smoluchowski equation governing the reaction kinetics is described by a time-independent diffusion equation as follows: where is the diffusion coefficient of reactant B and ( ⃗) is the concentration of reactant B at the position ⃗. is 1 (bulk-solution concentration), meaning that no reactions occur when B is far from A. For boundary (d), we can define an outer boundary; for example, in Figure 1, the blue hemisphere of radius represents an outer boundary condition, which is . Therefore, we solved the Smoluchowski equation for a bounded domain confined by the x-y plane containing reactant A at z = 0 and a hemisphere of radius with z > 0 (see the inner space created by the blue lines in Figure 1).In this work, we set a dimensionless value of to 10, as was used in our previous work [25], while a dimensionless length of each side of square reactant A was 1.8.
By solving Equation (1), we obtained , and from , we calculated the rate constant by calculating the fluxes of B at the reactive patches of A, which was given by where ⃗ ( ) is the surface area element of patch .Per usual [7], we finally calculated the rate constant normalized by the rate constant for the case in which the entire reactant surface was reactive.Specifically, this normalized rate constant was simply calculated by ⁄ and was denoted by .Here, we note that the maximum value of is 1.
In the first type of model, reactant A was localized on a plane, and the reactant had a finite surface area.However, when reactant A had a very large surface area or when a particular reactant area of interest was enlarged, we needed to consider only the reactant surface without other parts of the x-y plane; thus, we may assume on that length scale that the reactant surface was an infinite surface.In this case, we can consider the second type of model.In the second type, we removed the 4th boundary condition (d) and instead introduced the new boundary condition (d)', i.e., periodic boundary conditions (PBCs) for the x and y coordinates, which were typically used for molecular dynamic simulations [26], and an outer boundary condition for the z coordinate, which was ( ⃗)| = 1.Because of the PBCs, a reactive patch structure was periodically repeated in the x and y dimensions.From the boundary condition (d)', we, therefore, needed to solve the Smoluchowski equation for a bounded domain along the z-axis from the reactant surface at z = 0 to the outer boundary at z = , with the periodicity along the x and y axes (see the inner space of the rectangular prism at the center in Figure 2).In this work, we set a dimensionless value of to 10, while a dimensionless length of each side of the square reactant was 1.8. Figure 2 shows a schematic diagram of this second type of model.As with the first type of model, we solved the Smoluchowski equation or Equation (1) to obtain and calculated the rate constant by Equation (2) and eventually calculated the normalized rate constant or ⁄ ., where was the rate constant when the reactant surface was fully reactive.

Computational Methods
To solve the Smoluchowski equations for our models with N × N square reactants, we used the finite element method (FEM), a popular numerical method for solving partial differential equations with complicated boundary conditions.Specifically, for the FEM calculations, we used the software COMSOL [COMSOL Multiphysics ® v.5.4,COMSOL AB, Stockholm, Sweden], an FEM software package.In our models, to achieve high accuracy in the FEM calculations, we used very fine meshes in the regions of and near the reactive patches; in the mesh construction, the mesh for the region of the reactive patches was "extremely fine" in the predefined setting of COMSOL or was finer than the "extremely fine" setting, while the mesh for the remaining regions was "fine" or had a more refined mesh setting.With this mesh setting, the calculations with COMSOL yielded high accuracy.For example, in reaction models where the exact solutions were known, the relative error between the exact and numerical solutions was less than or approximately 1% [25].The FEM calculation was typically finished within a few minutes on a standard desktop computer, and even for demanding cases, it was finished at most within a few hours.COMSOL was also used to visualize the results shown in the figures with xmgrace (http://plasma-gate.weizmann.ac.il/Grace/), which was used for plotting and fitting curves to the data.

First Model Type: N× N SQUARE Reactant Model
In the N × N square reactant model of the first type, to study the effects of the size, number, and arrangement of reactive patches on the kinetics, we first considered various N × N patch configurations generated by the model.Since the N × N square model had N × N cells and each cell can be a reactive patch or not, the total number of configurations for the N × N square model was 2 × .For instance, for N = 2 and 3, the total number was 2 and 2 , respectively.In fact, in this work, we completely analyzed all the configurations for N = 2, 3, applying symmetric operations to the configurations.Note that a configuration of reactive patches can generate eight kinetically equivalent configurations by rotations and reflections through planes, as shown in Figure 3.However, when the configuration has symmetries in the arrangement of patches, the number of configurations generated is reduced because some configurations are duplicated.Therefore, we considered only symmetrically representative configurations with the number of equivalent configurations .In Figure 4, we present all the symmetrically representative configurations of patches for N = 2 and 3 and the normalized rate constants calculated from the solutions of the Smoluchowski equations.Since the rate constant was dependent on the distances between the reactive patches (cells), we calculated the average of the distances between all pairs of the reactive cells ̅ .Note that a similar study based on the geometric mean pair distance between patches has been performed [17].The average pairwise distances are also displayed in Figure 4. Here, we found a general trend that as ̅ increased, the normalized rate constant increased.Therefore, for a given number of reactive cells, the normalized rate constant could differ depending on the distances between the patches or the configurations of the patches.Notably, in the 3 × 3 square model, the configurations giving the highest rate constant were mainly the configurations that had a symmetric arrangement of the reactive patches (configurations 8, 23, 22, 16, and 3 for 2, 4, 5, 6, and 8 reactive cells, respectively).representative configurations from 0 cells (no reactive) to 9 cells (full reactive) are considered.For a given number of reactive cells, we arrange the configurations in order of increasing ̅ .Note that the unit of ̅ is the length of unit cell of the N × N square model.< > represents the ensembleaveraged value of for a given number of reactive cells.
Additionally, we calculated the ensemble-averaged value of the normalized rate constants < > for a given number of reactive cells, since we can determine the probability of each representative configuration from .From the calculation, we plotted the ensemble-averaged rate constant against the surface coverage ratio of reactive cells over the entire reactant area .was easily calculated by counting the number of reactive cells; for example, in the 2 × 2 model, if only one cell was reactive, was 0.25 because it was 1 over 4. The plot is shown in Figure 5a.As one can expect, as increased, the normalized rate constant increased, but as approached 1, the slope of the rate constant against decreased, which was due to the increase in the competition between reactive patches.We also show the maximum and the minimum values of the normalized rate constant for a given value of σ.Interestingly, as is shown in the inset of Figure 5a, the variation in the rate constant ∆ , i.e., the difference between the maximum and minimum values increased at low values of , reached its maximum, and then decreased at large values.This characteristic behavior was also observed in the reactive patches of a spherical reactant [7,17].Another feature from Figure 5a is that for a given value of , the ensemble-averaged normalized rate constant for the 3 × 3 model was larger than that for the 2 × 2 model.This outcome was consistent with the observation in a spherical reactant that a smaller patch gave a larger rate constant [7,17].Basically, this occurred because in the case of a higher value of N, more configurations existed in which the reactive patches were separated, which reduced the competition between patches.Therefore, one can expect that as N increases in the N × N model, one can obtain a higher ensembleaveraged rate constant.
In this work, we plotted the normalized rate constant as a function of the dimensionless distance ̅ / , where is the length of the unit cell.The result is shown in Figure 5b.Since the cases of the 2 × 2 square model were trivial, we presented only the cases of the 3 × 3 square model.From the plot, we found a general trend that as ̅ increased, the normalized rate constant increased.Specifically, from linear regressions, the result showed that a strong correlation existed between the normalized rate constant and ̅ , but the linear lines were not monotonously increasing, which means that ̅ was a good parameter but was not a perfect parameter for evaluating the competition between reactive patches in a given configuration.

Second Model Type: N × N Square Reactant Model
The main difference between the first and the second model types came from whether PBCs were applied to the x and y dimensions.In this section, we discuss the second model type with PBCs.When we applied PBCs to the x-y plane, the number of symmetrically representative configurations that we had to consider was significantly reduced, as shown in Figure 6.As in Figure 4, we calculated ̅ , , and for each representative configuration for the 2 × 2 and 3 × 3 square models.For a given number of reactive cells, we arranged the configurations in increasing order of ̅ .As expected, the normalized rate constant generally increased with ̅ .PBCs.All representative configurations from 0 cells (no reactive) to 9 cells (full reactive) are considered.For a given number of reactive cells, we arrange the configurations in order of increasing ̅ .Note that the unit of ̅ is the length of unit cell of the N × N square model.< > represents the ensemble-averaged value of for a given number of reactive cells.
The ensemble-averaged normalized rate constant < > is shown in Figure 7a.We also observed that as increased, the rate constant increased, which was the same trend observed in Figure 5a.However, differences exist.Notably, the normalized rate constant was very high even at low values of .That is, even a small surface coverage of reactive patches can produce a high rate constant.This significant increase due to the boundary condition is discussed further with the dependence of the kinetics on in Section 3.3.Another noticeable property from Figure 7a is that the variation in the rate constant was small, as indicated in the inset of Figure 7a.These interesting properties with the PBCs were also found in the reaction systems of a spherical reactant with multiple reactive patches [10].For a spherical reactant, small patches can produce a large rate constant, and the distribution of rate constants is narrow, which means that the difference between the maximum and minimum values is relatively small.These common features may come from the fact that the reactant surface with PBCs in the second model type was topologically similar to the surface over a sphere in that they were closed surfaces and have no boundary; if one walks straight starting from a certain position on the surface, one will eventually return to the starting position.As in Figure 5b, we also show a scatter plot of the normalized rate constant and the dimensionless average pairwise distance ̅ (PBC)/ for symmetrically representative configurations in the 3 × 3 square reactant model.The linear regression lines show a strong correlation between them. is easily calculated by the number of reactive cells over the total number of cells.(b) For numbers of reactive cells from two to nine in the 3 × 3 square model with PBCs, the plot of the normalized rate constant ( ) versus the dimensionless average pairwise distance ( ̅ (PBC)/ ), which is the average pairwise distance ( ̅ ) over the length of the unit cell ( ) in the presence of PBCs.The dotted lines represent linear regression lines obtained from the xmgrace program.In the linear regressions, the correlation coefficients are 0.978 (3 cells), 0.918 (4 cells), 0.966 (5 cells), and 0.997 (6 cells).

Dependence of Kinetics on the Fraction of the Total Reactive Patch Area
The results of the N × N square reactant models in Sections 3.1 and 3.2 indicate that the normalized rate constant is apparently dependent on the fraction of the total area of the reactive patches .For example, let us consider the first model type without PBCs.From Figure 4, if we consider Case 1 of the 3 × 3 model with one reactive cell ( = 1/9), Case 1 of the 2 × 2 model with one reactive cell ( = 1/4), Case 1 of the 3 × 3 model with four reactive cells ( = 4/9), and Case 1 of the 2 × 2 model with four reactive cells ( = 1), we understand the dependence of kinetics on in that the rate constant for a single square increases in a nonlinear way as σ increases; = 0 ( = 0), 0.323 ( = 1/9), 0.486 ( = 1/4), 0.654 ( = 4/9), and 1 ( = 4/4) (see the filled squares in Figure 8).Similarly, for the PBC cases, from Figure 6, if we consider Case 1 of the 3 × 3 model with one reactive cell ( = 1/9), Case 1 of the 2 × 2 model with one reactive cell ( = 1/4), Case 2 of the 3 × 3 model with four reactive cells ( = 4/9), and Case 1 of the 2 × 2 model with four reactive cells ( = 1), we see that the rate constant increases with in another nonlinear way; = 0 ( = 0), 0.892 ( = 1/9), 0.952 ( = 1/4), 0.982 ( = 4/9), and 1 ( = 4/4) (see the filled circles in Figure 8).However, because of the limitations in exploring various cases with different values of from the 2 × 2 and 3 × 3 square reactant models, we considered another square reactant model in which one reactive square cell was placed at the center of the square reactant.Here, we changed the size of the reactive cell, thus that we had different values of .The result is shown in Figure 8.We see that this result from the square model with one reactive cell was consistent with the results from the N × N square model from Figures 4 and 6 (filled symbols).The main feature is that as increased, the normalized rate constant increased, and the effect was much more significant with PBCs.That is, for a given value of , the rate constant in the presence of PBCs was larger than the corresponding rate constant without PBCs.In the presence of PBCs, even with = 0.01, the normalized rate constant was 61.7% of the maximum rate constant with = 1, while the corresponding rate constant in the first model type without PBCs was 9.78%.Additionally, we notice that in the inset of Figure 8, the concentration color at = 0.1 is almost blue (concentration is near zero) for the case with PBCs; in this case, the normalized rate constant was 88.2% of the maximum rate constant.
The high rate constants obtained with PBCs may suggest that in the first type of model, even in a single small patch in a reactant, if one breaks it into many smaller patches and arranges them such that any patch looks like the one with PBCs in Figure 8, the overall rate constant can be significantly increased.To check this idea, we considered the division-separation procedure of reactive cells in Sections 3.5 and 3.6 to see how the procedure enhances the rate constant.
Additionally, from the agreement of between the square model of a single reactive cell and the 2 × 2 and 3 × 3 square models in Figure 8, we note that the 2 × 2 and 3 × 3 square models were useful in obtaining a general trend of the dependence of the kinetics on , although the results did not cover the full range of .This result confirms the validity and usefulness of the simple N × N square reactant models.

Effect of the Reactive Patch Distribution on the Kinetics
The results of the N × N models in Sections 3.1 and 3.2 show that the rate constant is dependent on the distribution of the reactive patches, as demonstrated by the variation in the rate constant ∆ at a fixed value of in the insets of Figure 5a and Figure 7a.The general trend is that when the patches (cells) are aggregated, the normalized rate constant is lower than when the patches are largely separated.Here, we notice that the variation in the rate constant ∆ is dependent on the fraction of the area of the unit cell in the N × N model.That is, as N increases, decreases, which allows more diverse arrangements of patches on a reactant and increases the variation of the rate constant.Note that in the 2 × 2 square reactant model, is 1/4 and in the 3 × 3 model, is 1/9.Since the 2 × 2 and 3 × 3 square reactant models did not provide a variety of patch arrangements, we considered the spatial arrangements of patches with = 1/9 in the 9 × 9 square model with = 1/81.In this case, the number of reactive cells is 9.Note that in the 3 × 3 square model, the number of reactive cells corresponding to = 1/9 is 1.Moreover, note that in contrast to the 2 × 2 and 3 × 3 square models, in this case, considering all the possible patch configurations is practically impossible because the total number of arrangements of patches is enormous; when we ignore symmetrically equivalent patch arrangements and PBCs, the number of patch configurations is 81C9≈ 2.61 × 10 ; in fact, when we consider the symmetry of arrangements and PBCs, the actual number of symmetrically unique patch configurations is less than this number.Therefore, instead of considering all the cases, we considered 10 representative cases from the one with mostly aggregated patches to the one with largely separated patches, as shown in Figure 9a.We also used the average value of the pairwise distance as a way to measure the degree of scattering of patches over the reactant surface.For each case, we calculated the normalized rate constant for the first and second model types.We displayed the results as a plot of the normalized rate constant versus the dimensionless average pairwise distance ̅ / or ̅ (PBC)/ , as shown in Figure 9b.In Figure 9b, as one can expect from the competition between patches, as the reactive cells moved away from each other or ̅ increased, the normal rate constant increased.This effect is more apparent in the first type of model without PBCs.Interestingly, for a given value of = 1/9, by simply redistributing patches from the most aggregated patches (Case 1) to the maximally separated patches (Case 10), we can increase the normalized rate constant by more than two times, from 0.323 to 0.671; the normalized rate constant in Case 10 ( = 1/9) is comparable to the ensemble-averaged rate constant < > for three reactive cells ( = 1/3) in the 3 × 3 model (see Figure 4b; 0.653).Therefore, in addition to increasing the value, another way to increase the rate constant is to arrange the patches such that they are more separated from each other.
In the second type of model with PBCs, we also observed the trend that the rate constant increased with ̅ (PBC), as is indicated by the linear regression line.However, with PBCs, except for Case 1 (the most aggregated case), the normalized rate constant does not strongly depend on the details of the patch arrangements in that the values of are more or less the same.Another interesting observation was that the case giving the highest rate constant was Case 6 (the most symmetrical arrangement), where the reactive cells were uniformly distributed in the presence of PBCs with equal distance between neighboring cells.This arrangement of reactive cells corresponds to the symmetric distribution minimizing their competitions in a sphere, which was related to the solutions of the Thomson problem or the Tammes problem [7].Specifically, our competitionminimization problem with n reactive patches was similar to the Thomson problem aimed to find the spatial arrangement of n electrons on a sphere giving the minimum potential energy and the Tammes problem aimed to find the spatial arrangement of circles maximizing the minimal distance between circles in that their solutions are generally uniformly distributed arrangements.

Symmetric Reactive Patch Distributions
Based on the effect of the patch distribution in Section 3.4, we further studied the dependence of normalized rate constant on the number of reactive patches for a given value of .For this study, we employed the square reactant models of = 1/9 with symmetric patch distributions, in which the patches were uniformly distributed over the reactant surface (see Case 6 of Figure 9a).In the patch distributions, we considered = 1 2 , 2 , and 8 2 and arranged the patches symmetrically (see Figure 10a).Note that the case with = 3 2 was exactly the same as Case 6 of Figure 9a.For each case, we calculated the normalized rate constants for the first and second model types, the results of which we present in Figure 10b.The general trend was that the normalized rate constant increased with .As a result, for the first type of model, the case with = 8 2 and = 1/9 had a high rate constant (0.807), comparable to the rate constant (0.827) of the case with = 1 and = 0.7 (see Figure 8).However, for the second type of model, even for = 1, the normalized rate constant was very high (0.891); thus, the enhancement of due to the division and separation of patches was not as dramatic as that in the first type of model.Another feature in the plots of Figure 10b was that as increased, the effect of the enhancement of was reduced; the enhancement was largest when changed from 1 to 4.

Multiple Divisions of Patches and the Power-Law Relationship
For the square reactant model with symmetric patch distributions in Section 3.5, we can further discuss the mathematical relationship between the normalized rate constant and the number of reactive patches for the cases with a fixed value of .For the discussion, we considered sequential construction procedures for systematically generating the reactants with symmetric patch distributions of smaller patches.Precisely speaking, the procedures consisted of N steps of division and symmetric arrangement of patches, and the procedures were sequential in that the reactant at the Mth step was obtained from the reactant at the previous M-1th step.
To explain the procedure in detail, we started with a single square reactive patch at the center of the square reactant, as shown in the leftmost reactant in Figure 11a.In the first step, we divided the square reactive patch into four equal smaller patches.In Figure 11a, the division of the patch was represented by the dotted lines.We also performed the same type of division for the entire square reactant, resulting in four smaller square blocks.For each smaller square block, we placed one reactive patch at the center of the square block, resulting in the separation of the original reactive patch.The patch distribution obtained from this procedure is shown in the second reactant in Figure 11a, which corresponded to the symmetric patch distribution with = 2 = 4 in Figure 10a.Note that each block with a smaller reactive patch had exactly the same shape as that of the original initial reactant.That is, if we select any one of the four blocks in the second reactant and enlarged its area by four times, the enlarged block will be exactly the same as the original reactant (see the blue square lines in Figure 11a).Now, if we apply the same division-separation procedure to each of the four blocks in the second reactant in Figure 11a, we generate another four smaller blocks for each block.As a result, the total number of patches was = 2 × 2 =4 ; the overall patch distribution is shown in the third reactant in Figure 11a.We can repeatedly apply this procedure to a reactant at each step.Now, we consider how the overall normalized rate constant changes as we carry out the procedure repeatedly.For a quantitative discussion on the change in the rate constant due to the first division-separation procedure, we denoted the rate constant before the procedure as , the rate constant after the procedure as , and the overall increasing factor of the rate constant due to the reduction in the competition by the first procedure as .Then, the overall rate constant was = 4 = , where α > 1.Note that if the rate constant of each patch was simply proportional patch was divided into patches in the procedures, indicated how the normalized rate constant scales with (= ).In fact, from Figure 11b, we find power laws for by determining the values of and ; the slopes were 0.2158 and 0.02348 for the first and second model types, respectively, and accordingly, the values were 1.349 and 1.033, respectively.In fact, in the inset of Figure 11b, we plot the normalized rate constants calculated by the power laws of = / = ⁄ (dashed lines), along with the original ones shown in Figure 10b (solid lines).The two curves were generally in good agreement, but the rate constant from the power-law deviated from the original rate constant as increased.Instead of using four patches in the division-separation procedure, we can break a single patch into a smaller number of patches.In this case, decreases since the competition decreases less; for example, in a division-by-two procedure, the values ( = 10 ) for the first and second model types were 1.161 and 1.016, respectively.In contrast, if we break the patch into a larger number of patches, will increase since the competition decreases more.For example, in a division-by-eight procedure, the values ( = 10 ) for the first and second model types were 1.566 and 1.050, respectively.
We now consider the general properties of the power laws observed in Figure 11b.First, for all cases, was greater than 1, which means that always increases when the division and separation take place.Second, as increases, we expect that the slope ( ⁄ ) decreases and eventually becomes zero ( = 1) since has an upper bound of 1; otherwise, would be larger than 1 as increases.
3.6.2.Analysis for the Cases of = 0.01 and 0.001 In a reactant with symmetrically distributed reactive patches, we expected that as increased, the interference between a reactive patch and its other patches beyond the nearest neighboring patches was more significant.As a result, the increasing factor was more dependent on the number of patches, which means that the assumption that is constant, used to derive the power law of , becomes questionable.
In contrast, if decreases, the interference is less significant and the assumption that the increasing factor is constant makes more sense.To see if shows more of a power-law behavior for small values of , we considered the case of = 0.01 for a reactant with symmetrically distributed patches for the first and second model types.The result is shown in Figure 12.We also considered the case of = 0.001, but it was only for the second type of model; moreover, we calculated the normalized rate constant only up to = 25 because of the difficulty in creating a very refined mesh for the first and second model types.In Figure 12a, we display the normalized rate constant as a function of and in Figure 12b, we show the log-log plot of against .Remarkably, the correlation coefficient for the case of the first model type with = 0.01 was very high (0.997).In fact, the correlation coefficients for other cases were also high; 0.970 and 0.988 for the second model type with = 0.01 and 0.001, respectively.From the slopes, we also determined the values for the division-by-four procedure, i.e., 1.723 (first type, = 0.01), 1.1494 (second type, = 0.01), and 1.448 (second type, = 0.001).

Conclusions
In this work, we considered the diffusion-limited reactions between reactive patches on a reactant molecule fixed on a plane (x-y plane) and point-like partner reactant molecules and studied how the size, number, and spatial arrangement of reactive patches affected the overall reaction kinetics.For a systematic study, we proposed discrete N × N square reactant models, where the reactive patches are squares and are on the surface of a square reactant with an N × N grid; the patches occupy some of the N × N cells generated by the grid lines.Furthermore, we considered two types of reaction model.In the first model type, a reactant is placed on a chemically inert plane, and in the second, reactive patches are periodically placed on the entire reactant plane.Therefore, in the second type of model, we applied PBCs to the x and y coordinates.
For the 2 × 2 and the 3 × 3 square reactant models, we examined all possible cases of reactive patches from = 0 (non-reactive) to = 1 (fully reactive), where is the ratio of the total reactive patch area to the total reactant area.Through a complete study of the kinetics for the models, we found that as increases, the rate constant increases, but because of large competitions between patches, at large values, the increasing rate decreases.Additionally, we noticed that for a given value of , large variations in the rate constant could occur due to various arrangements of the reactive patches.This observation signifies that by rearranging patches, we may increase the rate constant.By comparing the ensemble-averaged rate constants in the 2 × 2 and 3 × 3 models, we see that for a given value of , the ensemble-averaged rate constant for the 3 × 3 model is larger than that for the 2 × 2 model, because, with smaller reactive patches, we can generate more patch arrangements that result in less competition between patches than with larger patches.Additionally, using a 9 × 9 model with = 1/9, we investigated the effect of the patch distribution on the kinetics by considering ten representative patch distributions.As one expects, when the patches are aggregated, the overall normalized rate constant is small, but when they are largely separated, the rate constant is large.In fact, for both the first and second model types, the patch distribution giving the lowest value of is the one where all reactive patches are maximally aggregated or compact.However, the patch distributions giving the highest value are different; the patch distribution for the first model type is the symmetric distribution where the reactive patches are maximally separated from each other (Case 10 in Figure 9a), and the patch distribution for the second model type is the symmetric distribution where the patches are uniformly distributed, and all the distances between nearest neighbors (including the PBC neighbors) are the same (Case 6 in Figure 9a).Interestingly, the latter patch distribution is similar to the uniform patch distribution in a sphere, maximizing the overall rate constant.
By further developing the idea of the symmetric patch distribution, we proposed a method to systematically increase the overall rate constant for a given value of by breaking a patch into smaller patches and separating them symmetrically based on the uniform arrangement of patches (Case 6 in Figure 9a).When we applied this division-separation procedure multiple times, we showed that the rate constant significantly increases, especially in the first type of model.Interestingly, we found that a power-law relationship exists between the normalized rate constant and the number of patches , and the relationship is more accurate when is smaller.One possible future research direction of this work is to show the usefulness of our discrete square reactant models.In this work, the ensemble used for calculating the averages is a microcanonical ensemble since we assume the probability of each configuration of patches being the same, meaning that the patches do not interact with each other.However, when the patches interact with themselves attractively or repulsively, the probability of each configuration depends on the arrangement of patches; thus, we may need to consider a canonical ensemble based on the potential energy of patches, which may be more realistic and more interesting.In that case, compared to continuous models, our discrete N × N square models are beneficial in that the number of patch configurations to be considered is finite, not infinite.Therefore, the study of the kinetics with the interactions between patches using our discrete models will be the subject of future investigations.

Figure 1 .
Figure 1.Schematic diagram of our first chemical reaction model.(Left) One reactant molecule A with an N × N cellular surface (black lines) is located at the center of the coordinate system and is on the xy plane at z = 0.The outer boundary condition is represented by the blue lines.(Right) The concentration profiles of reactant B from the Smoluchowski solutions for the N × N square models of N = 2 and 3 with certain reactive cell distributions are represented by color maps.The length of each side of reactant A is 1.8 ( : Unit length), and the distance from the center to the outer boundary (Router) is 10 .The red cells in the N × N square reactants indicate the reactive patches.

Figure 2 .
Figure 2. Schematic diagram of our second chemical reaction model.(Center) The arrangement of reactive patches on the bottom surface at z = 0 is repeated through periodic boundary conditions (PBCs) for the x-y dimensions.The blue top surface at z = corresponds to the outer boundary in the first model type in Figure 1, and the concentration on the blue surface is normalized to 1. (Left) The concentration profile of reactant B obtained from the Smoluchowski solution for the 3 × 3 square model with a certain reactive cell distribution is represented by a color map.(Right) The concentration profile of reactant B for the 2 × 2 square model with a certain reactive cell distribution is represented.In this unit cell, the length of each side of the bottom surface at z = 0 is 1.8 ( : Unit length), and the distance from the bottom surface to the top surface ( ) is 10 .The red cells in the N × N square reactants indicate the reactive patches.

Figure 3 .
Figure 3. Eight configurations are generated from the configuration of reactive patches shown on the left by symmetric operations of rotations around the z-axis out of the page and reflections through the planes represented by the dotted lines.

Figure 4 .
Figure 4. (a) Average pairwise distance between reactive patches (cells) ( ̅ ), the number of all configurations kinetically equivalent to a representative configuration ( ), and the normalized rate constant ( ) of a representative configuration in the 2 × 2 square model of the first type.All representative configurations from 0 cells (no reactive) to 4 cells (full reactive) are considered.(b) ̅ , , and of a representative configuration in the 3 × 3 square model of the first type.All

Figure 5 .
Figure 5. (a) Ensemble-averaged normalized rate constant < > as a function of the fraction of reactive cells over the entire reactant surface area ( ) for the 2 × 2 and 3 × 3 square reactant models of the first type.is easily calculated by the number of reactive cells over the total number of cells.(b) For the numbers of reactive cells from two to nine in the 3 × 3 square model, the plot of the normalized rate constant ( ) versus the dimensionless average pairwise distance ( ̅ / ), which is the average pairwise distance ( ̅ ) over the length of the unit cell ( ).The dotted lines represent linear regression lines obtained from the xmgrace program.In the linear regressions, the correlation coefficients are

Figure 6 .
Figure 6.(a) Average pairwise distance between reactive patches (cells) ( ̅ ), the number of all configurations kinetically equivalent to a representative configuration ( ), and the normalized rate

Figure 7 .
Figure 7. (a) Ensemble-averaged normalized rate constant < > as a function of the fraction of reactive cells over the entire reactant surface area ( ) for the 2 × 2 and 3 × 3 square models of the

Figure 8 .
Figure 8. Normalized rate constant as a function of the surface coverage ratio of reactive patch for a single reactive patch at the center of a reactant surface. is calculated as the reactive square area over the total reactant square area.The insets show the color maps of the concentration of reactant B on the reactant A surface, where a reactive patch lies without PBCs (first model type, upper panel) and with PBCs (second model type, lower panel).The filled squares and circles indicate the values from Figures 4 and 6, respectively.
Figure 8. Normalized rate constant as a function of the surface coverage ratio of reactive patch for a single reactive patch at the center of a reactant surface. is calculated as the reactive square area over the total reactant square area.The insets show the color maps of the concentration of reactant B on the reactant A surface, where a reactive patch lies without PBCs (first model type, upper panel) and with PBCs (second model type, lower panel).The filled squares and circles indicate the values from Figures 4 and 6, respectively.

Figure 9 .
Figure 9. (a) 10 representative arrangements of patches in the 9 × 9 square reactant model with the average pairwise distance between reactive patches (cells) ( ̅ ) for the first model type (without PBCs) and with ̅ (PBC) for the second model type (with PBCs).Here, is the length of the unit cell of the 9 × 9 square model.(b) The plot of the normalized rate constant versus the dimensionless average pairwise distance ̅ / for the first model type or ̅ (PBC)/ for the second model type, where is the length of the unit cell.Note that the cases giving the lowest and highest rate constants for the first type of model are Case 1 and Case 10.For the second type of model, Case 1 and Case 6 give the lowest and highest rate constants, respectively.The dotted lines represent linear regression lines obtained from the xmgrace program.In the linear regressions, the linear equations with correlation coefficients are = 0.0726 ( ̅ / ) + 0.217 with 0.990 for the first model type model and = 0.0303 ( ̅ (PBC)/ ) + 0.844 with 0.940 for the second model type.

Figure 10 .
Figure 10.(a) A square reactant model with symmetric reactive patch distributions based on the uniform patch distributions.In this model, once the number of patches is determined, the spatial arrangement of patches is completely determined.The arrangements of the reactive patches (red) for = 1 2 , 2 2 , 3 2 , 4 2 , 5 2 , 6 2 , 7 2 , and 8 2 are displayed.(b) The normalized rate constant as a function of for the first (red) and second (green) model types.The upper and lower insets present color maps of the concentration of reactant B for the first and second model types, respectively.

Figure 11 .
Figure 11.Sequential division-separation procedures of reactive patches on a reactant for a fixed value of (a) In the first division, a single reactive patch is broken into four smaller identical patches, and the smaller patches are separated and symmetrically arranged on the x-y plane.This divisionseparation procedure can be carried out any number of times .(b) Log-log plot of the normalized rate constant versus the number of patches for the case of = 1/9 shown in Figure 10b.From the linear regressions, we obtain the following linear equations for log ( ) versus log ( ): = 0.2158 − 0.4569 for the first model type (correlation coefficient: 0.985) and = 0.02348 − 0.04335 for the second model type (correlation coefficient: 0.960).These linear equations are represented by the dotted lines.The inset shows a comparison between from Figure 10b and from the linear regression.

Figure 12 .
Figure 12.Sequential division-separation procedures of reactive patches for the first model type with = 0.01 and the second model type (PBC) with = 0.01 and 0.001 (a) The normalized rate constant as a function of for the first (red) and second (green) model types.(b) Log-log plot of the normalized rate constant versus the number of patches .From the linear regressions, we obtain the following linear equations for log ( ) versus log ( ): = 0.3923 − 0.9914 for the first model type with = 0.01 (correlation coefficient: 0.997), = 0.1004 − 0.1869 for the second model type with = 0.01 (correlation coefficient: 0.970), and = 0.2672 − 0.4510 for the second model type with = 0.001 (correlation coefficient: 0.988).These linear equations are represented by the dotted lines.The inset shows a comparison between from (a) and from the linear regressions.