Integral equation study of effective attraction between like-charged particles mediated by cations: Comparison between IPY2 and HNC closures

Effective interactions between like-charged particles immersed in an electrolyte solution were calculated using two integral equation theories, hypernetted-chain (HNC)-Ornstein--Zernike (OZ) and ionic Percus--Yevick 2 (IPY2)-OZ. When the HNC-OZ theory was adopted, the electrolyte concentration dependence of the effective interaction showed a reentrant behavior. By contrast, the IPY2-OZ theory did not indicate the behavior. Monte Carlo simulations were performed for one of the model systems, and the results agreed qualitatively with those calculated using the HNC-OZ theory.


Introduction
Acidic proteins in an electrolyte solution attract each other under certain conditions [1][2][3][4][5]. The condensation behavior depends on the electrolyte concentration. When the electrolyte concentration is low, the proteins repel each other; but if the electrolyte concentration increases to some mM level, condensation appears, and with further concentration increase, the protein aggregations redissolve. This reentrant condensation behavior has been reported by Zhang et al. when the electrolyte solution has multivalent cations, such as Y 3+ [1][2][3][4][5].
To study effective interaction, theoretical approaches such as the integral equation theories (IETs) were examined . Some results indicated an effective attraction between like-charged particles. Despite the high costs due to the strong direct attraction between anions and cations, some researchers performed simulation studies. The results also showed the existence of effective attraction [10][11][12][13][14][15][16][17]. By contrast, the reentrant behaviors were not shown theoretically until 2011. A study using the hypernetted-chain (HNC)-Ornstein-Zernike (OZ) theory showed the reentrant behaviors in effective interaction [36,37], and some studies followed to describe the reentrant behavior in effective interaction [38][39][40]. Finally, a thermodynamic perturbation approach with effective interactions showed the reentrant condensations [41].
The results obtained from IET are affected by the closure approximation employed. Under each condition, the HNC approximation was shown to be valid for monovalent ionic systems such as the 1:1 electrolyte solution [43][44][45]. However, it was also shown that for highly associating ionic systems such as 2:2 electrolytes, the HNC results become unsatisfactory qualitatively. To overcome such obstacles, the ionic Percus-Yevick 2 (IPY2) closure relation was proposed [42]. The results obtained under the IPY2 approximation reproduced the Monte Carlo (MC) simulation data for 2:2 electrolyte solutions well [42]. In protein condensation, which is the main target of this study, a mixture of multivalent and monovalent ions is typically involved as an experimental condition. For such a case, it is still unclear which theory is more appropriate between the HNC and IPY2.
To study the existence of an effective attraction and reentrant behavior, we calculated the effective interaction between like-charged particles immersed in an electrolyte solution. In the present study, the results calculated using the IPY2-OZ and HNC-OZ theories are compared. If the results in the discussion of the adequate closure are contradictory, MC simulation will indicate which one is crucial.

Direct Interaction between Particles
A charged hard-sphere model was adopted in the present study. The direct interaction between particles is shown as follows: where u ij (r) is the interaction between particles i and j at the distance r, and Q i , ε 0 , ε r , and σ ij are the charge of particle i, the permittivity of vacuum, the relative permittivity of water (78.5), and the contact distance between i and j, respectively. The contact distance σ ij is (σ ii + σ jj )/2, where σ ii is the diameter of particle i.

Model A
The first model, Model A, was prepared to calculate the effective interaction between macroanions immersed in an electrolyte solution. This model was similar to that used in Ref. [37]. We deal with this model to compare our previous results, obtained using the HNC closure in Ref. [37], with those obtained using the IPY2 closure. In Model A, we adopted an explicit solvent model considering the solvent granularity. Thus, the solvent consists of uncharged hard spheres. The solution also has a 1:1 electrolyte. A pair of macroanions are immersed in the solution. "A", "C", "M", and "V" denote the anion, cation, macroanion, and solvent, respectively. The anions, cations, and solvents have the same size. The parameters are shown in Table 1. In this model, the total packing fraction is kept at 0.383, which is the packing fraction of water. This total packing fraction is maintained during the electrolyte concentration change. The temperature is 298 K. This is a simplified model of an aqueous electrolyte solution. We calculate the effective interactions between macroanions using the OZ equation with two closures, namely the HNC and the IPY2 closures.

Model B
We prepared the second model, Model B, to compare the radial distribution functions g(r) obtained using the IET and the MC simulations. Because the calculation cost for "Model A", that was mentioned above, was too expensive, the solvent hard spheres and the macroanions were removed from it. The results calculated using the HNC closure showed that the solvent hard spheres intensify the effective attraction between like-charged particles [39,40]. However, the explicit solvent model was not required, and the implicit solvent model also showed the effective attraction and the reentrant behavior in our studies [39,40]. Therefore, this model was selected to reduce the computational costs. In this model, the g AA (r)s values are compared. In the relationship of g(r) = exp[−βW (r)], g(r) and W (r) are the radial distribution function and the effective interaction, respectively. Thus, the calculated effective interactions can be compared. By contrast, the differences between the radial distribution functions were more apparent than those between the effective interactions in the important r region. Therefore, we decided to use g(r). We constructed this model referring to the previous study [39]. In this model, σ AA is 3.62Å and σ CC is 1.70Å. These values correspond to Cl − and Mg 2+ , whereas Q C is 2.4 e. Table 1 shows the parameters. When these parameters are adopted, the radial distribution function between anions can be expected to show the reentrant behavior [39].
where h ij (r), c ij (r), and ρ i are the total correlation function, the direct correlation function, and the number density of i, respectively. Here, ρ M goes to zero because we are studying the dilution limit of the macroanion. At the dilution limit, we obtain three OZ equations. One is the OZ equation for the bulk electrolyte solution as follows: and another is the OZ equation for the macroanion and electrolyte solution, as follows: At first, the OZ equation for the bulk system is solved with a closure relation. The calculated correlation functions become the input functions for the second OZ equation. Eq. (4) with a closure provides the correlation functions h M j (r), c M j (r). The effective interaction W M M (r) can be calculated using the correlation functions h M j (r), c M j (r), and the third OZ equation with a closure. The Fourier transformation is as follows: where , and c lM (r), respectively; and k is the wave number. When we obtain T M M (k) (or t M M (r)), we can calculate W M M (r) using a closure relation: Before the explanation of W M M (r), we introduce the closure relations. The closures are expected as follows: In the present study, we examined two closure relations: HNC and IPY2.
In the HNC, b(r) = 0 in Eq. (6). It has been used to study electrolyte solutions [26][27][28][29][30][31][32][33][34][36][37][38][39][40]. The equation is as follows: In the IPY2 closure equation, the bridge function b ij (r) for a like-charge pair is defined as follows [42]: whereas that for an unlike-charge pair is defined as follows: The definition of τ ij (r) is τ long ij (r) in Eq. (10) is given by In this equation, κ is Debye screening constant, which is given as This approximation was applied to some molten salt systems [42]. The form of W M M (r) = −k B T ln(h M M (r) + 1) was prepared as follows, Therefore, W M M (r) can be obtained without the iteration process using the abovementioned equation and t M M (r), which is obtained from the third OZ equation (Eq. (5)).

The IETs and MC Simulation to Calculate Model B
The radial distribution function g AA (r) = h AA (r) + 1 is analyzed in Model B. We calculate g AA (r)s using the IETs. In this case, we adopt the OZ equation, namely Eq. (3) for the bulk. Moreover, Model B does not have explicit solvent particles.
The OZ equation is solved with the HNC and the IPY2 closures. MC simulations were performed only for Model B. Because of the Coulomb system, we adopted the Ewald method in the potential energy calculation. The sampling of the Coulomb system was difficult. Therefore, the thermal replica-exchange method was adopted to avoid the trapping at the local minimum. The temperatures of replicas were 298 K, 320 K, 340 K, and 360 K. The replicas m and n were exchanged using the transition probability [46,47], P acc mn : where U m is the potential energy of replica m.
To discuss the electrolyte concentration dependence, the molar concentrations C A = 0.01 M, 0.1 M, 0.3 M, 0.4305 M, and 1 M were examined. According to Ref. [39], the reentrant behavior can be expected in the concentration region for this model. The charge neutrality condition was maintained. Thus, the anion and cation numbers in the MC simulation (N A and N C ) were 60 and 25, respectively (Q A = − e, Q C = 2.4 e). The basic cell sizes were set according to the concentration.

Results and Discusstions
The effective interactions between two macroanions are shown in Fig. 1. Figure 1 (a) presents the effective interactions obtained by the HNC approximation reported in [37], and Fig. 1(b) presents those obtained by the IPY2 approximation. The results will be briefly discussed, as per Ref. [37]. When the electrolyte concentration was 10 interaction changed from repulsive to attractive, and the attraction disappeared as the electrolyte concentration increased. Thus, under the HNC, a clear reentrant association was observed between two macroanions. By contrast, under the IPY2 closure, we obtained the result shown in Fig.  1 (b). When the electrolyte concentration was 10 −5 M, a significant repulsion was observed between the two macroanions. As the electrolyte concentration increased, the repulsion gradually disappeared, whereas a clear attraction was not observed in the concentration considered in this study. The results under the IPY2 approximation were qualitatively different from those under the HNC approximation.
Here, we define the dimerization free energy of macroanions as the contact value of the effective interaction between macroanions. The electrolyteconcentration dependence of the dimerization free energy is shown in Fig. 2. Figures 2 (a) and (b) correspond to the results obtained by HNC and IPY2, respectively. Under the HNC approximation ( Fig. 2 (a)), the dependence of the dimerization free energy on the concentration was not monotonic: the dimerization free energy decreased with increasing electrolyte concentration, to take a local minimum in the vicinity of 10 −3 M, followed by a significant increase. In the low concentration region (i.e., from 10 −5 M to 10 −4 M), the dimerization free energy was positive, whereas in the middle concentration region (i.e., in the vicinity of 10 −3 M), it took a negative value, which was more than 20 times the thermal energy, k B T . When the concentration became high (i.e., higher than 10 −1 M), the dimerization free energy was nearly equal to zero. The negative dimerization free energy in the vicinity of the local minimum observed in Fig. 2 (a) indicated that the effective interaction between macroanions was attractive, as mentioned above. The curve drawn in Fig. 2 (a) clearly shows the reentrant behavior of the effective interaction under the HNC. By contrast, the dimerization free energy of the effective interaction under the IPY2 approximation decreased monotonically with concentration ( Fig. 2 (b)). In the high concentration region, such as the vicinity of 10 0 M, the dimerization free energy was negative (at best twice the thermal energy, k B T ). Such a clear attraction as observed under the HNC approximation did not appear when the IPY2 closure was adopted. Moreover, no reentrant behavior was perceived under the IPY2 approximation for the effective interaction between macroanions. These results were in contrast to those of the HNC closure.
The closures showed qualitatively different results: the HNC-OZ theory presented effective attractions between macroanions and the reentrant behavior, whereas the IPY2-OZ theory showed no such behaviors. The effective repulsion between macroanions decreased monotonically with the increasing electrolyte concentration, indicating only screening by the electrolyte in the IPY2-OZ theory. In the present study, we adopted molecular simulations to study this difference. We surveyed models and thermodynamic conditions for reentrant behavior using the HNC-OZ theory under conditions that could be calculated by molecular simulation [39]. We obtained the correlation functions for Model B using the HNC-OZ and IPY2-OZ theories and MC simulation. To discuss the validity of the IETs, we compared the radial The radial distribution functions between Cl − anions g AA (r) are shown in Fig. 3. Figures 3 (a), (b), and (c) show g AA (r) obtained using the HNC-OZ theory, the IPY2-OZ theory, and the MC simulation, respectively. When the electrolyte concentration was lower than 1.0 × 10 −3 M, all g AA (r)s were smaller than 1, meaning that the effective interactions between anions were repulsive. As the concentration increased, g AA (r) became larger. When the concentration was 1.0×10 −2 M, only the term g AA (r) for the HNC-OZ theory had values that were larger than 1. Others, namely (b) and (c), followed this behavior, and all g AA (r)s had values that were larger than 1 for 1.0 × 10 −1 M in the vicinity of the contact distance, 5.5Å. This behavior, namely the disappearance of the direct repulsion caused by the increase in electrolyte, is common to all calculation methods when the concentration is lower than 5.0 × 10 −2 M. Although we found a common behavior under the lower concentration condition, we also found that the HNC approximation emphasized effective attraction in g AA (r). The concentration dependence of the "first peak" value is shown in Fig.4. (The definition of the "first peak" is given later.) The first peak for the HNC approximation was the highest in g AA (r) at 5.0 × 10 −2 M. This behavior seems to be due to the well-known short-ranged overcorrelation yielded for many model systems for the HNC approximation [45,[48][49][50]. In particular, the first peak of the spatial distribution function is emphasized when using the approximation. By contrast, the values for the IPY2 approximation were closer to those for the MC simulation at 5.0 × 10 −2 M. The values for IPY2 were smaller than those of MC. The deviations caused by the IPY2 and HNC approximations were opposite, and the absolute values of the deviations were at the same level. Therefore, there was no superiority between the two closures in this concentration region. The dependence on the electrolyte concentration qualitatively depended on the approximation when the concentration became higher than 5.0 × 10 −2 M. g AA (r)s for the HNC-OZ theory and the MC simulation indicated the reentrant behavior. Their first peaks decreased as the concentration increased to certain levels. In Fig. 3, we also found the highest first peaks for the HNC-OZ theory and MC simulation at 5.0×10 −2 M and 3.0×10 −1 M, respectively. Under the high concentration condition, g AA (r) for the HNC-OZ theory became closer to that for the MC simulation. For example, g AA (r) for the HNC-OZ theory was almost the same as that for the MC simulation when the concentration was 3.0×10 −1 M. By contrast, the value of g AA (r) for the IPY2-OZ theory did not describe reentrant behavior. The concentration dependence was monotonic. Although the concentration dependence for the IPY2-OZ theory had a different qualitative behavior from that for the MC simulation, g AA (r)s also got closer to each other at 4.305 × 10 −1 M.
The previous paragraph concluded on the similarity and the difference of g AA (r)s obtained using three methods with Model B. They are clearly shown in Fig. 4. However, we must explain how the plots were created because the first peaks of g AA (r) are unclear under the dilute condition. For example, the first peaks are unclear when the concentration is 1.0×10 −3 M. We determined the separation distances between two anions for the first peak using the most attractive g AA (r)s for each method. The separations of the first peak for the HNC, IPY2, and MC simulation were 5.39Å, 5.61Å, and 5.60Å, respectively. Based on this definition, we obtained the plots in Fig. 4. The dots (the IPY2-OZ theory) did not show the reentrant behavior, although others (the HNC-OZ theory (the solid line) and the MC simulation (the squares)) showed it. At the beginning of this section, we discussed the results for Model A. The difference between the HNC and the IPY2 approximations is similar to that between them in Model B. Our MC simulations suggest that the HNC-OZ theory gives qualitatively reasonable results even when the effective attraction between anions appears. Under this condition, the cations mediate the effective attraction [37]. Therefore, it seems that the mediate mechanism must be included in the HNC approximation, whereas the mediation effect is excessively emphasized in the approximation.
According to Ref. [42], the IPY2-OZ theory gave reasonable correlation functions for highly associating ionic systems such as the 2:2 electrolytes, whereas for monovalent ionic systems the HNC-OZ theory was known to be reasonable [40,43,44]. The model adopted in this study was a mixture of monovalent and multivalent ions, which is an intermediate system between the abovementioned references [40,[42][43][44]. Moreover, this study has dealt with the effective attraction between like-charged particles, which was not considered in the previous studies [40,[42][43][44]. Based on the abovementioned discussion, we can expect the following conclusion. The HNC is better for the mixture of monovalent and multivalent ions than the IPY2. In particular, it seems that the IPY2 approximation includes an insufficient mediation effect. As a result, the approximation did not show the effective attraction and the reentrant behavior.

Conclusion
We calculated the effective interaction between like-charged particles immersed in an electrolyte solution using two IETs, namely HNC-OZ and IPY2-OZ. Two models were adopted, Models A and B. In Model A, two highly charged macroanions (−10 e) were immersed in a 1:1 electrolyte solution [37]. The explicit solvent model was adopted in this model. By contrast, in Model B, the system had small monovalent anions and small multivalent cations. The implicit solvent model was adopted in the latter.
The electrolyte concentration dependences of the effective interaction were examined. In Model A, the results of the HNC-OZ theory indicated a reentrant behavior. As the concentration increases, the effective interaction turns from repulsive to attractive and returns to repulsive in the reentrant behavior (see Figs. 1 (a) and 2 (a)). By contrast, we did not find the reentrant behavior in the results of the IPY2-OZ theory (see Figs. 1 (b) and 2 (b)). This closure dependence of the reentrant behavior was the same in Model B (Figs. 3 and 4).
We also performed a MC simulation for Model B. The simulation results qualitatively agreed with those of the HNC-OZ theory. We conclude that the HNC-OZ theory provides adequate results qualitatively. Our previous studies show that the counter-charged ions cause an effective attraction between likecharged particles at the mediation site [36][37][38][39][40]. The simulation results suggest that the IPY2 approximation included an insufficient mediation effect. By contrast, the HNC approximation emphasizes the effect [45,[48][49][50].