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

Abstract

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 the competitions between reactive patches arising in real applications.

1. 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 curvature-dependent 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, such 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.

2. Reaction Models and Computational Methods

2.1. Reaction Models

In this investigation, we use two types of reaction models to study a diffusion-limited (or diffusion-controlled) bimolecular reaction between reactants A and B. Here, reactant A is the one we are most interested in because it is a square reactant with a reactive patch distribution, while reactant B is simply a point-like particle. One type is for the situation where reactant A is localized on a certain chemically-inert plane. The other type is for the situation where reactant A has a large surface area or is enlarged for detained analysis; thus, its reactive patches are 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 $D$. 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 diffusion-controlled) 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, such that the concentration of B or $c_B$ 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:

$$D{\nabla }^2{c}_B\left(\overrightarrow{r}\right)=0,$$

(1)

where $D$ is the diffusion coefficient of reactant B and $c_B\left(\overrightarrow{r}\right)$ is the concentration of reactant B at the position $\overrightarrow{r}$.

To solve Equation (1), we set up the boundary conditions as follows: (a) on the reactive patches of A, $c_B$ is zero (absorbing boundary condition), or ${c_B\left(\overrightarrow{r}\right)|}_{reactivepatches}=0$, where $\overrightarrow{r}$ refers to the reactive patches; (b) on the remaining non-reactive surface area of A, the flux of B onto A is zero (reflecting boundary condition), or ${\widehat{n}\left(\overrightarrow{r}\right)·D\nabla c_B\left(\overrightarrow{r}\right)|}_{non-reactivearea}=0$, where $\widehat{n}\left(\overrightarrow{r}\right)$ is the normal vector from the reactant surface at $\overrightarrow{r}$; (c) on the area of the x-y plane outside reactant A, the flux of B onto the area is also zero (reflecting boundary condition), or ${\widehat{n}\left(\overrightarrow{r}\right)·D\nabla c_B\left(\overrightarrow{r}\right)|}_{x-yplaneoutsidereactantA}=0$; (d) at locations far from A, $c_B$ 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 $R_{outer}$ represents an outer boundary condition, which is ${c_B\left(\overrightarrow{r}\right)|}_{\left|\overrightarrow{r}\right|=R_{outer}andz>0}=1$. Therefore, we solve the Smoluchowski equation for a bounded domain confined by the x-y plane containing reactant A at z = 0 and a hemisphere of radius $R_{outer}$ with z > 0 (see the inner space created by the blue lines in Figure 1). In this work, we set a dimensionless value of $R_{outer}$ to 10, as was used in our previous work [25], while a dimensionless length of each side of square reactant A is 1.8.

By solving Equation (1), we obtain $c_B$, and from $c_B$, we calculate the rate constant $k$ by calculating the fluxes of B at the reactive patches of A, which is given by

$$k={\displaystyle \sum }_i{{\displaystyle \int }}^{}D\nabla c_B\left(\overrightarrow{r}\right)·d{\overrightarrow{S}}_{A\left(i\right)},$$

(2)

where $d{\overrightarrow{S}}_{A\left(i\right)}$ is the surface area element of patch $i$.

Per usual [7], we finally calculate the rate constant normalized by the rate constant $k_{entire}$ for the case in which the entire reactant surface is reactive. Specifically, this normalized rate constant is simply calculated by $k/k_{entire}$ and was denoted by $\overline{k}$. Here, we note that the maximum value of $\overline{k}$ is 1.

In the first type of model, reactant A is localized on a plane, and the reactant has a finite surface area. However, when reactant A has a very large surface area or when a particular reactant area of interest is enlarged, we need 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 is an infinite surface. In this case, we can consider the second type of model. In the second type, we remove the 4th boundary condition (d) and instead introduce the new boundary condition (d)’, i.e., periodic boundary conditions (PBCs) for the x and y coordinates, which are typically used for molecular dynamic simulations [26], and an outer boundary condition for the z coordinate, which is ${c_B\left(\overrightarrow{r}\right)|}_{z=R_{outer}}=1$. Because of the PBCs, a reactive patch structure is periodically repeated in the x and y dimensions. From the boundary condition (d)’, we, therefore, need 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 =$R_{outer}$, 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 $R_{outer}$ to 10, while a dimensionless length of each side of the square reactant is 1.8.

Figure 2 shows a schematic diagram of this second type of model. As with the first type of model, we solve the Smoluchowski equation or Equation (1) to obtain $c_B$ and calculate the rate constant $k$ by Equation (2) and eventually calculate the normalized rate constant $\overline{k}$ or $k/k_{entire}$., where $k_{entire}$ is the rate constant when the reactant surface is fully reactive.

2.2. Computational Methods

To solve the Smoluchowski equations for our models with N × N square reactants, we use the finite element method (FEM), a popular numerical method for solving partial differential equations with complicated boundary conditions. Specifically, for the FEM calculations, we use 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 use 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 is “extremely fine” in the predefined setting of COMSOL or is finer than the “extremely fine” setting, while the mesh for the remaining regions is “fine” or has a more refined mesh setting. With this mesh setting, the calculations with COMSOL yield high accuracy. For example, in reaction models where the exact solutions are known, the relative error between the exact and numerical solutions is less than or approximately 1% [25]. The FEM calculation is typically finished within a few minutes on a standard desktop computer, and even for demanding cases, it is finished at most within a few hours. COMSOL is also used to visualize the results shown in the figures with xmgrace (http://plasma-gate.weizmann.ac.il/Grace/), which is used for plotting and fitting curves to the data.

3. Results and Discussion

3.1. 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 consider various N × N patch configurations generated by the model. Since the N × N square model has N × N cells and each cell can be a reactive patch or not, the total number of configurations for the N × N square model is $2^{N\times N}$. For instance, for N = 2 and 3, the total number is $2^4$ and $2^9$, respectively. In fact, in this work, we completely analyze 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 consider only symmetrically representative configurations with the number of equivalent configurations $n_c$.

In Figure 4, we present all the symmetrically representative configurations of patches for N = 2 and 3 and the normalized rate constants $\overline{k}$ calculated from the solutions of the Smoluchowski equations. Since the rate constant is dependent on the distances between the reactive patches (cells), we calculate the average of the distances between all pairs of the reactive cells ${\overline{d}}_2$. 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 find a general trend that as ${\overline{d}}_2$ increases, the normalized rate constant $\overline{k}$ increases. Therefore, for a given number of reactive cells, the normalized rate constant $\overline{k}$ 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 are mainly the configurations that have a symmetric arrangement of the reactive patches (configurations 8, 23, 22, 16, and 3 for 2, 4, 5, 6, and 8 reactive cells, respectively).

Additionally, we calculate the ensemble-averaged value of the normalized rate constants $\langle \overline{k}\rangle$ for a given number of reactive cells, since we can determine the probability of each representative configuration from $n_c$. From the calculation, we plot the ensemble-averaged rate constant against the surface coverage ratio of reactive cells over the entire reactant area $\sigma$. $\sigma$ is easily calculated by counting the number of reactive cells; for example, in the 2 × 2 model, if only one cell is reactive, $\sigma$ is 0.25 because it is 1 over 4. The plot is shown in Figure 5a. As one can expect, as $\sigma$ increases, the normalized rate constant $\overline{k}$ increases, but as $\sigma$ approaches 1, the slope of the rate constant against $\sigma$ decreases, which is 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 $\mathsf{\sigma }$. Interestingly, as is shown in the inset of Figure 5a, the variation in the rate constant $∆\overline{k}$, i.e., the difference between the maximum and minimum values increases at low values of $\sigma$, reaches its maximum, and then decreases 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 $\sigma$, the ensemble-averaged normalized rate constant for the 3 × 3 model is larger than that for the 2 × 2 model. This outcome is consistent with the observation in a spherical reactant that a smaller patch gives a larger rate constant [7,17]. Basically, this occurrs because in the case of a higher value of N, more configurations exist in which the reactive patches are separated, which reduces the competition between patches. Therefore, one can expect that as N increases in the N × N model, one can obtain a higher ensemble-averaged rate constant.

In this work, we plot the normalized rate constant $\overline{k}$ as a function of the dimensionless distance ${\overline{d}}_2/a$, where $a$ is the length of the unit cell. The result is shown in Figure 5b. Since the cases of the 2 × 2 square model are trivial, we present only the cases of the 3 × 3 square model. From the plot, we find a general trend that as ${\overline{d}}_2$ increases, the normalized rate constant $\overline{k}$ increases. Specifically, from linear regressions, the result shows that a strong correlation exists between the normalized rate constant $\overline{k}$ and ${\overline{d}}_2$, but the linear lines are not monotonously increasing, which means that ${\overline{d}}_2$ is a good parameter but is not a perfect parameter for evaluating the competition between reactive patches in a given configuration.

3.2. Second Model Type: N × N Square Reactant Model

The main difference between the first and the second model types comes from whether PBCs are applied to the x and y dimensions. In this section, we discuss the second model type with PBCs. When we apply PBCs to the x-y plane, the number of symmetrically representative configurations that we have to consider is significantly reduced, as shown in Figure 6. As in Figure 4, we calculate ${\overline{d}}_2$, $n_c$, and $\overline{k}$ for each representative configuration for the 2 × 2 and 3 × 3 square models. For a given number of reactive cells, we arrange the configurations in increasing order of ${\overline{d}}_2$. As expected, the normalized rate constant $\overline{k}$ generally increases with ${\overline{d}}_2$.

The ensemble-averaged normalized rate constant $\langle \overline{k}\rangle$ is shown in Figure 7a. We also observe that as $\sigma$ increases, the rate constant increases, which is the same trend observed in Figure 5a. However, differences exist. Notably, the normalized rate constant $\overline{k}$ is very high even at low values of $\sigma$. 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 $\sigma$ in Section 3.3. Another noticeable property from Figure 7a is that the variation in the rate constant is 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 is topologically similar to the surface over a sphere in that they are 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 $\overline{k}$ and the dimensionless average pairwise distance ${\overline{d}}_2\left(\mathrm{PBC}\right)/a$ for symmetrically representative configurations in the 3 × 3 square reactant model. The linear regression lines show a strong correlation between them.

3.3. Dependence of Kinetics on the Fraction of the Total Reactive Patch Area

The results of the N × N square reactant models in Section 3.1 and Section 3.2 indicate that the normalized rate constant $\overline{k}$ is apparently dependent on the fraction of the total area of the reactive patches $\sigma$. 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 ($\sigma =1/9$), Case 1 of the 2 × 2 model with one reactive cell ($\sigma =1/4$), Case 1 of the 3 × 3 model with four reactive cells ($\sigma =4/9$), and Case 1 of the 2 × 2 model with four reactive cells ($\sigma =1$), we understand the dependence of kinetics on $\sigma$ in that the rate constant for a single square increases in a nonlinear way as $\mathsf{\sigma }$ increases; $\overline{k}$= 0 ($\sigma =0$), 0.323 ($\sigma =1/9$), 0.486 ($\sigma =1/4$), 0.654 ($\sigma =4/9$), and 1 ($\sigma =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 ($\sigma =1/9$), Case 1 of the 2 × 2 model with one reactive cell ($\sigma =1/4$), Case 2 of the 3 × 3 model with four reactive cells ($\sigma =4/9$), and Case 1 of the 2 × 2 model with four reactive cells ($\sigma =1$), we see that the rate constant $\overline{k}$ increases with $\sigma$ in another nonlinear way; $\overline{k}$= 0 ($\sigma =0$), 0.892 ($\sigma =1/9$), 0.952 ($\sigma =1/4$), 0.982 ($\sigma =4/9$), and 1 ($\sigma =4/4$) (see the filled circles in Figure 8).

However, because of the limitations in exploring various cases with different values of $\sigma$ from the 2 × 2 and 3 × 3 square reactant models, we consider another square reactant model in which one reactive square cell is placed at the center of the square reactant. Here, we change the size of the reactive cell, such that we have different values of $\sigma$. The result is shown in Figure 8. We see that this result from the square model with one reactive cell is consistent with the results from the N × N square model from Figure 4 and Figure 6 (filled symbols). The main feature is that as $\sigma$ increases, the normalized rate constant $\overline{k}$ increases, and the effect is much more significant with PBCs. That is, for a given value of $\sigma$, the rate constant in the presence of PBCs is larger than the corresponding rate constant without PBCs. In the presence of PBCs, even with $\sigma$= 0.01, the normalized rate constant $\overline{k}$ is 61.7% of the maximum rate constant with $\sigma$= 1, while the corresponding rate constant in the first model type without PBCs is 9.78%. Additionally, we notice that in the inset of Figure 8, the concentration color at $\sigma$= 0.1 is almost blue (concentration is near zero) for the case with PBCs; in this case, the normalized rate constant $\overline{k}$ is 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 consider the division-separation procedure of reactive cells in Section 3.5 and Section 3.6 to see how the procedure enhances the rate constant.

Additionally, from the agreement of $\overline{k}$ 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 are useful in obtaining a general trend of the dependence of the kinetics on $\sigma$, although the results do not cover the full range of $\sigma$. This result confirms the validity and usefulness of the simple N × N square reactant models.

3.4. Effect of the Reactive Patch Distribution on the Kinetics

The results of the N × N models in Section 3.1 and Section 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 $∆\overline{k}$ at a fixed value of $\sigma$ in the insets of Figure 5a and Figure 7a. The general trend is that when the patches (cells) are aggregated, the normalized rate constant $\overline{k}$ is lower than when the patches are largely separated. Here, we notice that the variation in the rate constant $∆\overline{k}$ is dependent on the fraction of the area of the unit cell ${\sigma }_1$ in the N × N model. That is, as N increases, ${\sigma }_1$ 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, ${\sigma }_1$ is 1/4 and in the 3 × 3 model, ${\sigma }_1$ is 1/9.

Since the 2 × 2 and 3 × 3 square reactant models do not provide a variety of patch arrangements, we consider the spatial arrangements of patches with $\sigma$= 1/9 in the 9 × 9 square model with ${\sigma }_1$= 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 $\sigma$= 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 ₈₁C₉$\approx 2.61\times {10}^{11}$; 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 consider 10 representative cases from the one with mostly aggregated patches to the one with largely separated patches, as shown in Figure 9a. We also use 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 calculate the normalized rate constant $\overline{k}$ for the first and second model types. We display the results as a plot of the normalized rate constant $\overline{k}$ versus the dimensionless average pairwise distance ${\overline{d}}_2/a$ or ${\overline{d}}_2\left(\mathrm{PBC}\right)/a$, as shown in Figure 9b.

In Figure 9b, as one can expect from the competition between patches, as the reactive cells move away from each other or ${\overline{d}}_2$ increases, the normal rate constant $\overline{k}$ increases. This effect is more apparent in the first type of model without PBCs. Interestingly, for a given value of $\sigma$= 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 $\overline{k}$ by more than two times, from 0.323 to 0.671; the normalized rate constant $\overline{k}$ in Case 10 ($\sigma$ = 1/9) is comparable to the ensemble-averaged rate constant $\langle \overline{k}\rangle$ for three reactive cells ($\sigma$= 1/3) in the 3 × 3 model (see Figure 4b; 0.653). Therefore, in addition to increasing the $\sigma$ 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 observe the trend that the rate constant $\overline{k}$ increases with ${\overline{d}}_2\left(\mathrm{PBC}\right)$, as is indicated by the linear regression line. However, with PBCs, except for Case 1 (the most aggregated case), the normalized rate constant $\overline{k}$ does not strongly depend on the details of the patch arrangements in that the values of $\overline{k}$ are more or less the same. Another interesting observation is that the case giving the highest rate constant is Case 6 (the most symmetrical arrangement), where the reactive cells are 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 is related to the solutions of the Thomson problem or the Tammes problem [7]. Specifically, our competition-minimization problem with n reactive patches is 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.

3.5. Symmetric Reactive Patch Distributions

Based on the effect of the patch distribution in Section 3.4, we further study the dependence of normalized rate constant $\overline{k}$ on the number of reactive patches $n_p$ for a given value of $\sigma$. For this study, we employ the square reactant models of $\sigma$= 1/9 with symmetric patch distributions, in which the patches are uniformly distributed over the reactant surface (see Case 6 of Figure 9a). In the patch distributions, we consider $n_p$= 1², 2², 3², 4², 5², 6², 7², and 8² and arrange the patches symmetrically (see Figure 10a). Note that the case with $n_p$= 3² is exactly the same as Case 6 of Figure 9a. For each case, we calculate the normalized rate constants $\overline{k}$ for the first and second model types, the results of which we present in Figure 10b. The general trend is that the normalized rate constant $\overline{k}$ increases with $n_p$. As a result, for the first type of model, the case with $n_p=$8² and $\sigma$= 1/9 has a high rate constant (0.807), comparable to the rate constant (0.827) of the case with $n_p=$1 and $\sigma$= 0.7 (see Figure 8). However, for the second type of model, even for $n_p=$1, the normalized rate constant $\overline{k}$ is very high (0.891); thus, the enhancement of $\overline{k}$ due to the division and separation of patches is not as dramatic as that in the first type of model. Another feature in the plots of Figure 10b is that as $n_p$ increased, the effect of the enhancement of $\overline{k}$ is reduced; the enhancement is largest when $n_p$ changed from 1 to 4.

3.6. 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 $\overline{k}$ and the number of reactive patches $n_p$ for the cases with a fixed value of $\sigma$. For the discussion, we consider sequential construction procedures for systematically generating the reactants with symmetric patch distributions of smaller patches. Precisely speaking, the procedures consist of N steps of division and symmetric arrangement of patches, and the procedures are sequential in that the reactant at the Mth step is obtained from the reactant at the previous M-1th step.

To explain the procedure in detail, we start 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 divide the square reactive patch into four equal smaller patches. In Figure 11a, the division of the patch is represented by the dotted lines. We also perform the same type of division for the entire square reactant, resulting in four smaller square blocks. For each smaller square block, we place 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 corresponds to the symmetric patch distribution with $n_p$= $2^2$ = $4^1$ in Figure 10a. Note that each block with a smaller reactive patch has 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 enlarge its area by four times, the enlarge 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 is $n_p$= $2^2\times 2^2$=$4^2$; 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 $\overline{k}$ 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 denote the rate constant before the procedure as $k_0$, the rate constant after the procedure as $k_1$, and the overall increasing factor of the rate constant due to the reduction in the competition by the first procedure as ${\alpha }_1$. Then, the overall rate constant $k_1$ is $k_1=4\frac{k_0}{4}{\alpha }_1=k_0{\alpha }_1$, where ${\mathsf{\alpha }}_1>1$. Note that if the rate constant of each patch is simply proportional to its area, the overall rate constant $k_1$ is $k_1=4\frac{k_0}{4}=k_0$, since there are four patches with a rate constant of $\frac{k_0}{4}$. Therefore, the reduction in competition is taken into account only through the factor ${\alpha }_1$, which may be different from system to system.

Similarly, from the second division-separation procedure, the competitions among the patches are further reduced. After the second procedure, we obtain the third patch distribution, in which 16 identical blocks are generated. Each block has a reactive square patch at its center. Because of the way we construct the reactant, a similarity symmetry exist between a four-block square in the third patch distribution and a four-block square (entire reactant) in the second patch distribution (see the green square lines in Figure 11a). Because of the similarity symmetry appearing in the sequential procedures, we may assume that the increasing factor ${\alpha }_2$ due to the reduction in the competition by the second procedure is the same as the increasing factor ${\alpha }_1$ by the first procedure. Under the assumption that ${\alpha }_1={\alpha }_2=\cdots =\alpha$, the overall rate constant after two divisions, or $k_2$, can be written as $k_2=4\frac{k_1}{4}{\alpha }_2=k_1{\alpha }_2=k_1\alpha$. Together with $k_1=k_0{\alpha }_1=k_0\alpha$, we finally obtain $k_2=k_1\alpha =\left(k_0\alpha \right)\alpha =k_0{\alpha }^2$.

In this way, we can repeat the division-separation procedure under the similarity symmetry, the competition is continuously reduced, and as a result, the overall rate constant increases. Moreover, with the assumption that $\alpha$ is constant, after an M-fold symmetric division and arrangement, we finally obtain the overall rate constant $k_M=k_0{\alpha }^M$. That is, we obtain the power–law relationship between the rate constant $k_M$ and the number of division-separation procedures $M$. Subsequently, we obtain the normalized rate constant ${\overline{k}}_M$, which is calculated by dividing $k_M$ by the rate constant for $\sigma =1$ ($k\left(\sigma =1\right)$ or $k_{entire}$), and the same power law ${\overline{k}}_M={\overline{k}}_0{\alpha }^M$, where ${\overline{k}}_0=k_0/k\left(\sigma =1\right)$. However, note that this power–law scaling behavior has limitations. That is, we consider only the nearest neighbor competitions through the similarity symmetry, but if we take into account the competitions arising from all possible pairs of patches, the behavior of the rate constant will deviate from the power–law behavior as M or $n_p$ increases. Moreover, since the upper value of ${\overline{k}}_M$ is 1 (normalized rate constant for $\sigma =1$), we expect that as the number of procedures M increases, ${\overline{k}}_M$ will eventually become saturated, which means that $\alpha$ decreases; otherwise ${\overline{k}}_M$ diverges as M increases.

In the case above of the division of patches by four, since the number of patches $n_p$ and the number of procedures $M$ are related by $n_p=4^M$, we can represent ${\overline{k}}_M$ in terms of $n_p$ as ${\overline{k}}_M={\overline{k}}_0{\alpha }^M={\overline{k}}_0{\alpha }^{{\log }_{10}{n}_p/{\log }_{10}4}$. That is, ${\log }_{10}{\overline{k}}_M={\log }_{10}{\overline{k}}_0+\frac{{\log }_{10}\alpha }{{\log }_{10}4}{\log }_{10}{n}_p$, which means that ${\log }_{10}{\overline{k}}_M$ is proportional to ${\log }_{10}{n}_p$. Additionally, note that from a logarithmic property, we can rewrite ${\overline{k}}_M$ as ${\overline{k}}_M={\overline{k}}_0{n}_p{}^{{\log }_{10}\alpha /{\log }_{10}4}$. Therefore, the simple scaling argument provides a power law of $n_p$ for the rate constant ${\overline{k}}_M$, generally denoted by $\overline{k}$.

3.6.1. Analysis for the Case of $\sigma =$ 1/9

For the reactants with symmetric distributions in Figure 10a, to see what the plot of ${\log }_{10}\overline{k}$ versus ${\log }_{10}{n}_p$ looks like according to the above division-by-four procedure, we calculate ${\log }_{10}\overline{k}$ for the rate constant $\overline{k}$ shown in Figure 10b and display the log-log plot in Figure 11b. In the plot, the y-intercept and the slope correspond to ${\log }_{10}{\overline{k}}_0$ and $\frac{{\log }_{10}\alpha }{{\log }_{10}4}$, respectively. Again, the physical interpretation of slope ($s$) as $\frac{{\log }_{10}\alpha }{{\log }_{10}4}$ means that if we increase $n_p$ by $4^M$ times (e.g., $n_p$= 1, 4, 16, 64, …), the increasing factor $\alpha$ is given by ${10}^{s{\log }_{10}4}$. However, the trend deviates from linearity as ${\log }_{10}{n}_p$ increases. In fact, to find the slope, we fit the data with a linear equation for ${\log }_{10}\overline{k}$ versus ${\log }_{10}{n}_p$ using the xmgrace program. The correlation coefficients for the first and the second model types in Figure 10 are 0.985 and 0.960, respectively. These high values explain why the relation between $\overline{k}$ and $n_p$ can be described by power laws ($\overline{k}~{\alpha }^{{\log }_{10}n/{\log }_{10}4}$ or $\overline{k}~n_p{}^{{\log }_{10}\alpha /{\log }_{10}4}$). From the linear equations for the regressions, we can easily determine $\alpha$ for the division-by-four procedure; for a determined value of $s$, $\alpha$ is given by ${10}^{s{\log }_{10}4}$. Note again that here, the factor of ${\log }_{10}4$ is due to the division by four carried out in the procedure. In general, if we break a patch into $n$ smaller patches at every step of the division, $\mathrm{then}\alpha ={10}^{s{\log }_{10}n}$. In this case, when each patch was divided into $n$ patches in the procedures, $\alpha$ indicate how the normalized rate constant $\overline{k}$ scales with $n_p$ (=$n^M$).

In fact, from Figure 11b, we find power laws for $\overline{k}$ by determining the values of $\alpha$ and ${\overline{k}}_0$; the slopes are 0.2158 and 0.02348 for the first and second model types, respectively, and accordingly, the $\alpha$ values are 1.349 and 1.033, respectively. In fact, in the inset of Figure 11b, we plot the normalized rate constants $\overline{k}$ calculated by the power laws of $\overline{k}={\overline{k}}_0{\alpha }^{{\log }_{10}{n}_p/{\log }_{10}4}={\overline{k}}_0{n}_p{}^{{\log }_{10}\alpha /{\log }_{10}4}$ (dashed lines), along with the original ones shown in Figure 10b (solid lines). The two curves are generally in good agreement, but the rate constant from the power–law deviates from the original rate constant as $n_p$ increases.

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, $\alpha$ decreases since the competition decreases less; for example, in a division-by-two procedure, the $\alpha$ values ($\alpha ={10}^{s{\log }_{10}2}$) for the first and second model types are 1.161 and 1.016, respectively. In contrast, if we break the patch into a larger number of patches, $\alpha$ will increase since the competition decreases more. For example, in a division-by-eight procedure, the $\alpha$ values ($\alpha ={10}^{s{\log }_{10}8}$) for the first and second model types are 1.566 and 1.050, respectively.

We now consider the general properties of the power laws observed in Figure 11b. First, for all cases, $\alpha$ is greater than 1, which means that $\overline{k}$ always increases when the division and separation take place. Second, as $n_p$ increases, we expect that the slope ($d\overline{k}/d{n}_p$) decreases and eventually becomes zero ($\alpha =1$) since $\overline{k}$ has an upper bound of 1; otherwise, $\overline{k}$ would be larger than 1 as $n_p$ increases.

3.6.2. Analysis for the Cases of $\sigma =$ 0.01 and 0.001

In a reactant with symmetrically distributed reactive patches, we expect that as $\sigma$ increases, the interference between a reactive patch and its other patches beyond the nearest neighboring patches is more significant. As a result, the increasing factor $\alpha$ is more dependent on the number of patches, which means that the assumption that $\alpha$ is constant, used to derive the power law of $\overline{k}$, becomes questionable.

In contrast, if $\sigma$ decreases, the interference is less significant and the assumption that the increasing factor $\alpha$ is constant makes more sense. To see if $\overline{k}$ shows more of a power–law behavior for small values of $\sigma$, we consider the case of $\sigma$= 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 consider the case of $\sigma$= 0.001, but it is only for the second type of model; moreover, we calculate the normalized rate constant $\overline{k}$ only up to $n_p$= 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 $\overline{k}$ as a function of $n_p$ and in Figure 12b, we show the log-log plot of $\overline{k}$ against $n_p$. Remarkably, the correlation coefficient for the case of the first model type with $\sigma$= 0.01 is very high (0.997). In fact, the correlation coefficients for other cases are also high; 0.970 and 0.988 for the second model type with $\sigma$= 0.01 and 0.001, respectively. From the slopes, we also determine the $\alpha$ values for the division-by-four procedure, i.e., 1.723 (first type, $\sigma$= 0.01), 1.1494 (second type,$\sigma$= 0.01), and 1.448 (second type,$\sigma$= 0.001).

4. 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 affect 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 $\sigma$= 0 (non-reactive) to $\sigma$= 1 (fully reactive), where $\sigma$ 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 $\sigma$ increases, the rate constant increases, but because of large competitions between patches, at large $\sigma$ values, the increasing rate decreases. Additionally, we noticed that for a given value of $\sigma$, 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 $\sigma$, 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 $\sigma$= 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 $\overline{k}$ 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 $\sigma$ 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 $\overline{k}$ and the number of patches $n_p$, and the relationship is more accurate when $\sigma$ 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.