Osmotic second virial coefficients for hydrophobic interactions as a function of solute size Free

Hydrophobic interactions are the effective interactions between nonpolar molecules, or between nonpolar groups in molecules, in aqueous solutions.¹ The self-assembly of amphiphilic molecules, such as surfactants, lipids, and biological molecules, in an aqueous solution illustrates that hydrophobic interactions between nonpolar groups are one of the driving forces for the formation of micelles, bilayers, and other multimolecular aggregates in aqueous environments.^2–4 

Quantification of the effective interactions between hydrophobic molecules of various sizes and shapes in liquid solutions is, by far, incomplete compared to that of the molecular interactions in gas phases. Direct force measurements are limited to those between macroscopic hydrophobic surfaces;⁵ measurements of osmotic second virial coefficients B are available for proteins and other amphiphilic molecules in water; however, experimental data of B for simple, nonpolar species in water, which will be the bases for understanding a variety of hydrophobic interactions, are still sparse and so theoretical and molecular simulation studies play significant roles. The recent progress in understanding the nature of hydrophobic interactions is, for example, in reviews by Berne et al., ⁶ by Pratt et al.,⁷ and by Ben-Amotz.^8,9

The hydrophobic interactions between small molecules at ambient conditions are not as strong as the dispersion forces in a vacuum. For example, the osmotic second virial coefficient B of methane in water was evaluated as a function of temperature, pressure, and salt concentration based on molecular simulation of realistic model systems.^10–12 A conclusion derived from those studies is that around and below room temperature, at atmospheric pressure, and in pure water, the hydrophobic interactions are less attractive than the direct interactions. Recent theoretical studies also showed that the water-induced parts of the effective interactions between hydrophobic groups in proteins act as repulsive forces.^13,14

Hydrophobic pair interactions are precisely described by the potential w(r) of mean force between two solute molecules distant r apart. The overall magnitudes of pair interactions are quantified by the osmotic second virial coefficient B, as it is related to w(r) by

with k being Boltzmann’s constant, T being the temperature, ρ being the number density of the solute, and dτ being the volume element. The integral is over the whole space. This expression was originally obtained by McMillan and Mayer¹⁵ and is now known as one of the Kirkwood–Buff formulas.¹⁶ 

Size dependences of the effective pair potential w(r) between particles in liquids have been an important subject of theoretical studies of liquids,^17–23 most of which examine hard-sphere mixtures whose structure and phase behavior are purely determined by entropic effects. In 1999, Lum, Chandler, and Weeks (LCW) developed a theory of solvation of solutes in water and showed that there is a crossover in the size dependence of the solvation free energy of hydrophobic solutes between small (angstroms) and large (nanometers) length scales, thus leading to large attractive forces between hydrophobic surfaces on the large length scale.²⁴ Since around that time, the question of how the hydration free energy and the hydrophobic interaction change with the molecular size of hydrophobes has been a topic of great interest,^3,25–37 and substantial progress has been made in understanding the hydration free energy.^{25,27–29,31,32,38} In contrast, a number of questions regarding the hydrophobic interactions remain unanswered.

Here, we calculate the osmotic second virial coefficients B of hydrophobic solutes in water to answer the question of how B depends on molecular diameter σ of the solutes. We know the answer if the solvent is “vacuum,” i.e., the system is a dilute gas of the solutes. In that case, B is just the second virial coefficient B_gas of that gas, which is given by B_gas = −(1/2)∫[e^−ϕ(r)/kT − 1]dτ, the same form as Eq. (1) with the effective pair potential w(r) now replaced by the direct pair potential ϕ(r). For typical pair potential functions ϕ(r), B_gas is proportional to σ³ if the well depth of ϕ(r) is fixed. Some illustrations of B_gas are given in the supplementary material. If water or any other solvent exists around solute molecules, there appears the solvent-induced pair potential w(r) − ϕ(r), which depends in some nontrivial way on σ as well as temperature, pressure, and composition of the solvent. It is this part of w(r) that would presumably make solute-size (σ) dependence of B qualitatively different from that of B_gas.

To compute w(r), and ultimately B, molecular dynamics simulations are performed using GROMACS 2018.1³⁹ for the model systems consisting of 4000 TIP4P/2005⁴⁰ water molecules and 10, 20, or 40 Lennard-Jones (LJ) solute particles. The smallest solute has the LJ parameters for methane:⁴¹ σ_CH4 = 0.373 nm and ɛ_CH4 = 1.23 kJ/mol. The largest one has the diameter σ = 2σ_CH4. The energy parameters ɛ for all the solutes are set to ɛ_CH4. We adopt three key strategies of evaluating B accurately: first, replacing the solute–solute LJ pair potential ϕ(r) in MD simulations by some repulsive pair potentials, e.g., the repulsive part of Weeks–Chandler–Andersen (WCA) potential,⁴² to suppress aggregation of solute particles during simulation and later constructing the desired g(r); second, correcting the finite size effect on B employing the method proposed by Krüger et al.;^43,44 third, correcting the finite solute-concentration effect on B using an approximation proposed previously.^45,46 Computational details are described in the supplementary material.

Figure 1(a) shows the radial distribution functions g(r) for pairs of hydrophobic particles with different diameters in water at 300 K and at 1 bar. The corresponding potentials w(r) of mean force are plotted in the inset. Remember that the solute–solute LJ energy parameters ɛ are set to the common value ɛ_CH4 and so do the solute–solvent energy parameters. Nevertheless, the first peak of g(r) grows rapidly with the peak position shifting to large r as the solute size σ increases and correspondingly the minimum of the potential curve w(r) decreases (becomes more negative). Since the well depth ɛ of the solute–solute LJ pair potential is being kept fixed when the solute size σ increases, those changes in g(r) and in w(r) are entirely due to the solute-size effect on the water-mediated interaction.

The osmotic second virial coefficients B resulting from those g(r) or w(r) are plotted in Figs. 1(b) and 1(c). The red solid curves in the figures are the second virial coefficients B_gas of the LJ gases at the same temperature. As shown in Fig. 1(b), B and B_gas are both negative and decreasing functions of the solute diameter σ. However, the rate of change in B with σ is greater than that in B_gas. That is, as the molecular size of the solute increases, the hydrophobic pair interactions become attractive more rapidly than the direct pair interactions. More precise information on the solute-size dependence of B is found from the log–log plot in Fig. 1(c). The relation between B and σ on the log–log graph is near-linear with the slope close to 6. This suggests

The constant A is negative for the LJ solutes in water. The behavior of the osmotic second virial coefficient is significantly different from that of the second virial coefficient B_gas for the LJ pair potential

where A_LJ is a function of reduced temperature kT/ɛ alone and is negative below the Boyle temperature and positive otherwise. The analytical expression of A_LJ is available.⁴⁷ 

We also calculated the osmotic second virial coefficients at higher temperatures 330 and 360 K in addition to 300 K. Figure 2(a) shows log–log plots of B vs σ. Power law behaviors of B are found at the higher temperatures, too. At 330 K, B ∼ Aσ⁶ fits very well the numerical data while at 360 K, the fifth or sixth power of σ is the best fit. We note that the effect of temperature on B at any given σ agrees with the results reported earlier:^10,11 B becomes more negative with increasing T, i.e., the hydrophobic interaction is enhanced at higher temperatures.

We also considered another set of solutes, which have stronger attraction with each other and with water. To be specific, they are LJ particles interacting with each other with the energy parameter ɛ = 2ɛ_CH4 = 2.46 kJ/mol and interacting with water molecules with the energy parameter greater by a factor of $2$ than before. Figure 2(b) shows the log–log plot of B vs σ for the two sets of solutes. Unlike B for the hydrophobic solutes (blue circles in the figure), B for the solutes with stronger attraction (red squares) is not fit by a single power law. Instead, there seems to be a crossover on the solute sizes of σ/σ_CH4 between 1.5 and 1.6: the cubic power law fits B for small length scales, while the sixth power law seems to follow for larger ones. Further calculations are required to confirm the sixth power law for larger solutes. The cubic power law observed for smaller solutes is the same as the size dependence of B_gas of the LJ gases, i.e., in the absence of the solvent-mediated interaction. Thus, it is conjectured that the cubic power law appears when the direct pair interaction between solute molecules is sufficiently stronger than the solvent-mediated interaction. To verify the conjecture, we examined the effect of the direct pair interaction on the behavior of B. We then found that the exponent α in B ∼ Aσ^α changes from 6 to 3 as the solute–solute LJ energy parameter increases from ɛ_CH4 to 3ɛ_CH4. See details and the plots of B vs σ for three sets of solutes with ɛ/ɛ_CH4 = 1, 2, and 3 in the supplementary material: effects of the solute–solute direct interaction on B.

Our main goal is to examine the solute-size dependence on B, a measure of the overall strength of the hydrophobic interaction, and now we find that B for the hydrophobic solutes in water obeys B ∼ Aσ^α with α = 6 for T = 300 and 330 K and α = 5 or 6 for 360 K. Then, naturally, the following question would arise: Is the power law with α ≃ 6 characteristic of the hydrophobic interaction between solute molecules in water? As the first attempt to answer the above question, we calculated B for the LJ solutes in a nonpolar solvent. The solvent molecules are now LJ particles whose size is similar to that of water.^48,49 The potential parameters are given in the supplementary material. The number density of the nonpolar solvent is fixed to that of water at 300 K and at 1 bar. Figure 2(c) displays the plots of B for the LJ particles in the nonpolar solvent and in water. First, we find that B > 0 in the nonpolar solvent, indicating that the effective solute–solute pair interactions are overall repulsive. This means that the LJ solutes are not solvophobic with respect to the nonpolar solvent examined here. It should be noted, however, that the sign of B is sensitively dependent on the combination of the solute–solute, solute–solvent, and solvent–solvent interactions as well as the thermodynamic state. Second, as shown in the log–log plot in the inset, B for the nonpolar solutes in the nonpolar solvent obeys the power law with the same exponent α = 6 as that found for those solutes in water. The above results suggest that whether B is positive or negative, it follows the power law with α ≃ 6 as long as the contribution of the solvent-induced pair interaction is dominant over that of the direct pair interaction.

The osmotic second virial coefficient given by Eq. (1) is the correlation function integral over the whole space: $B=−(1/2)∫0∞h(r)4πr2dr$⁠, where h(r) = g(r) − 1 = e^−w(r)/kT − 1. Then, one can ask whether the power law dependence [Eq. (2)] arises from h(r) at short distances r or over long distances. To answer this, we calculate the fraction B^R of B that comes from the integral of h(r) from r = 0 to a finite distance R. When R is chosen to be the minimal distance such that h(r) = 0 or equivalently w(r) = 0, which is close to σ, then B^R ∝ σ³. This is because B^R here is essentially proportional to the excluded volume of each solute molecule. Now, if R is the distance of the first local minimum of h(r) or equivalently the first local maximum of w(r), then B^R obeys nearly the same power law as B does. The plots of B^R are given in the supplementary material, Fig. S3. Thus, the characteristic solute-size dependence of B is already anticipated from the short-range solute–solute correlation.

Further study is needed to examine both how robust the power law [Eq. (2)] is in aqueous and non-aqueous solutions and how the power α ≃ 6 is rationalized by theory of liquids. Now, we just note the following thermodynamic identity:⁵⁰ 

where B″ is the coefficient in the expansion of the osmotic pressure Π = ρkT(1 + B″ρ + ⋯) just like the definition of B, but now the expansion is at a fixed solvent density instead of a fixed solvent chemical potential; v is the partial molecular volume of the solute, and χ is the isothermal compressibility of the pure solvent. It was shown earlier⁵¹ that B″ = −(1/2)∫c(r)dτ with c(r) being the solute–solute direct correlation function. The second term on the right-hand side of Eq. (4), which is negative, is asymptotically proportional to σ⁶ for large σ. This is because v ≫ kTχ for typical liquids and v ∝ σ³. Cerdeiriña and Widom have shown using equations of state that B″ (>0) largely cancels the second term in Eq. (4) to give B, which is much smaller in magnitude than them.⁵² Values of B and −(v − kTχ)²/2kTχ obtained from our MD simulations are consistent with their findings. It is then concluded that B″, too, must be asymptotically proportional to σ⁶. Just to note an exceptional case, as the critical point of the solvent is approached, the second term would be diverging while B″ remains finite, and so B would be essentially equal to the second term.

We have seen that values of B for the solutes in the nonpolar solvent are positive [Fig. 2(c)]. Here, we note that the magnitudes of the first and second terms in Eq. (4) are also different from those for water. The compressibility χ of the nonpolar solvent is larger than that of water, and, as a consequence, (v − kTχ)²/2kTχ is only 38% of that of water. From the numerical data of B and (v − kTχ)²/2kTχ, it is concluded that the first term B″ (>0) in Eq. (4), too, is smaller in the nonpolar solvent than in water, but the difference is so small that B becomes positive in the nonpolar solvent.

In summary, for model systems of hydrophobic solutes in water and in a simple liquid, we have obtained numerically accurate values of the osmotic second virial coefficients B as a function of solute molecular diameter σ using MD simulations with corrections to finite-size and finite-concentration effects. It is well known that the second virial coefficient B_gas for LJ potentials is proportional to σ³ with the proportionality constant being a function of ɛ/kT. In contrast, the numerical data of B (<0) for LJ particles, whose σ ranges from σ_CH4 to 2σ_CH4 and whose ɛ is fixed to ɛ_CH4, in water are best fit by Aσ^α with α ≃ 6. The same power law behavior has been found for the LJ particles in an LJ solvent in which the effective pair interaction is overall repulsive (B > 0). It remains to be examined by simulation or molecular theory of liquids whether or not the sixth power law holds beyond the range of σ examined (1.2 ≤ σ/σ_v ≤ 2.4 for TIP4P/2005 water and 1.4 ≤ σ/σ_v ≤ 2.9 for the LJ solvent, with σ_v being the LJ diameter of solvent molecules). We note that the second term in Eq. (4), which is proportional to σ⁶ for large solutes, is characteristic of the osmotic second virial coefficient; if the solvent turns to an infinitely dilute gas, that term vanishes and B = B″ = B_gas. Therefore, what one may conclude from Eq. (4) is that if B is asymptotically a power function of σ, the exponent α must be equal to or greater than 6. This conjecture will be tested or refined by studying different models ranging from hard-sphere mixtures to solutions of simple solutes with varying solute–solvent interactions. We hope to examine the robustness of the power law behavior or the possibility of even stronger solute-size dependences of B in future work.

See the supplementary material for illustrations of B_gas, computational details, the effect of the solute–solute direct interaction on B, and solute-size dependence arising from short-range correlation.

This work was supported by the JSPS KAKENHI (Grant Nos. 18KK0151, 18K03562, and 20H02696). Part of the computation was carried out using the Research Center for Computational Science, Okazaki, Japan (Project No. 21-IMS-C125).