Origin of hydrophobicity and enhanced water hydrogen bond strength near purely hydrophobic solutes

Significance Hydrophobicity governs a wide range of fundamental physicochemical processes, but its physical origin is unclear. The classical explanation of hydrophobicity is that tiny “icebergs” are formed near such solutes; however, no experimental proof has been advanced for their existence. Here, we used four small purely hydrophobic solutes (methane, ethane, krypton, and xenon) in water to study hydrophobicity at the most fundamental level. We present unequivocal experimental proof for strengthened water hydrogen bonds near purely hydrophobic solutes, matching those in ice and clathrates. The water molecules involved in the enhanced hydrogen bonds display extensive structural ordering resembling that in clathrates, thus indicating the fundamental interconnection between electrostatic screening (shielding) and the hydrophobic effect. Hydrophobicity plays an important role in numerous physicochemical processes from the process of dissolution in water to protein folding, but its origin at the fundamental level is still unclear. The classical view of hydrophobic hydration is that, in the presence of a hydrophobic solute, water forms transient microscopic “icebergs” arising from strengthened water hydrogen bonding, but there is no experimental evidence for enhanced hydrogen bonding and/or icebergs in such solutions. Here, we have used the redshifts and line shapes of the isotopically decoupled IR oxygen–deuterium (O-D) stretching mode of HDO water near small purely hydrophobic solutes (methane, ethane, krypton, and xenon) to study hydrophobicity at the most fundamental level. We present unequivocal and model-free experimental proof for the presence of strengthened water hydrogen bonds near four hydrophobic solutes, matching those in ice and clathrates. The water molecules involved in the enhanced hydrogen bonds display extensive structural ordering resembling that in clathrates. The number of ice-like hydrogen bonds is 10–15 per methane molecule. Ab initio molecular dynamics simulations have confirmed that water molecules in the vicinity of methane form stronger, more numerous, and more tetrahedrally oriented hydrogen bonds than those in bulk water and that their mobility is restricted. We show the absence of intercalating water molecules that cause the electrostatic screening (shielding) of hydrogen bonds in bulk water as the critical element for the enhanced hydrogen bonding around a hydrophobic solute. Our results confirm the classical view of hydrophobic hydration.

Hydrophobicity plays an important role in numerous physicochemical processes from the process of dissolution in water to protein folding, but its origin at the fundamental level is still unclear. The classical view of hydrophobic hydration is that, in the presence of a hydrophobic solute, water forms transient microscopic "icebergs" arising from strengthened water hydrogen bonding, but there is no experimental evidence for enhanced hydrogen bonding and/or icebergs in such solutions. Here, we have used the redshifts and line shapes of the isotopically decoupled IR oxygen-deuterium (O-D) stretching mode of HDO water near small purely hydrophobic solutes (methane, ethane, krypton, and xenon) to study hydrophobicity at the most fundamental level. We present unequivocal and model-free experimental proof for the presence of strengthened water hydrogen bonds near four hydrophobic solutes, matching those in ice and clathrates. The water molecules involved in the enhanced hydrogen bonds display extensive structural ordering resembling that in clathrates. The number of ice-like hydrogen bonds is 10-15 per methane molecule. Ab initio molecular dynamics simulations have confirmed that water molecules in the vicinity of methane form stronger, more numerous, and more tetrahedrally oriented hydrogen bonds than those in bulk water and that their mobility is restricted. We show the absence of intercalating water molecules that cause the electrostatic screening (shielding) of hydrogen bonds in bulk water as the critical element for the enhanced hydrogen bonding around a hydrophobic solute. Our results confirm the classical view of hydrophobic hydration.
hydrophobic hydration | hydrogen bonding | IR spectroscopy | electrostatic screening | ab initio molecular dynamics D espite its great importance in numerous phenomena, the origin of hydrophobicity remains one of the most disputed topics in science (1)(2)(3)(4)(5). Experimental studies have shown that small purely hydrophobic solutes (alkanes and noble gases) in water increase the order (6, 7) and restrict the mobility (8) of neighboring water molecules. There are several opposing views on how to explain these data. The classical view is that a solute modifies water structure by forming transient, semiordered clathrate-like clusters ("icebergs"; used here only as a loose term) around it, arising from enhanced water hydrogen bonding (H bonding) (6,7). This enhancement is brought about by either strengthening (9) or increasing the number of water to water H bonds (10). The classical view explains the characteristic changes in the thermodynamic variables of hydrophobic hydration (positive ΔG, ΔC p , negative ΔS, and ΔH) and the restricted mobility of water molecules observed by NMR (8). Neutron diffraction (11,12) and extended X-ray absorption fine structure (EXAFS) (13) studies, however, show that the water molecules around small purely hydrophobic molecules do not differ significantly from those in pure liquid water. According to the dynamic view, the hydrophobic solute causes slowdown of the dynamics of the nearby water molecules by obstructing the jump mechanism of rotational relaxation, whereas water structure and H-bonding strength remain basically unchanged (5,14). The main problem of these two views is that there is essentially no direct experimental support for enhanced H bonding and resulting formation of icebergs in the hydration shell of small purely hydrophobic solutes (15) or a jump mechanism of rotational relaxation of the water molecules (5). An alternative view, arguing that the characteristic changes in thermodynamic variables of hydrophobic hydration can be explained by the difficulty of a solute molecule to be accommodated in an appropriate cavity in water because of small size of the water molecules, was proposed by Lee (16). Baldwin (17) recently suggested that extensive ordering of water molecules in the hydration shell of hydrophobic molecules is caused by van der Waals attraction between hydrophobic C and water O atoms.
Experimental proof for strengthened water H bonds and icebergs near purely hydrophobic solutes (alkanes and noble gases) would confirm unequivocally the classical view. Vibrational spectroscopy of the oxygen-hydrogen (O-H) stretching mode (ν OH ) is commonly used as the most reliable, sensitive, and model-free approach for determining relative strengths of H bonds (18)(19)(20)(21)(22)(23)(24). Any frequency downshift (redshift) of the ν OH mode, detected by comparing spectra of bulk water with those of the water molecules perturbed by solute, would indicate enhanced H-bond strength. Moreover, the spectral line width of ν OH is commonly used to assess structural order: a narrower line width indicates greater structural order (23,25). Measuring redshifts and line widths of the water molecules perturbed by hydrophobic solutes is very difficult because of the latter's very low solubility and further complicated by the formation of solid clathrates at the high pressures and low temperatures used to increase solubility. Vibrational spectra of the water molecules perturbed by small purely hydrophobic solutes have not been measured. However, the redshift of water ν OH has been observed in the hydration shell of several amphipathic solutes at high concentrations (20)(21)(22)24). These measurements are, however, not convincing proof of enhanced water H-bond strength near hydrophobic solutes, because polar parts of the molecules and concentration effects could affect the results considerably (14,26). Here, we used the redshifts and line shapes of the isotopically Significance Hydrophobicity governs a wide range of fundamental physicochemical processes, but its physical origin is unclear. The classical explanation of hydrophobicity is that tiny "icebergs" are formed near such solutes; however, no experimental proof has been advanced for their existence. Here, we used four small purely hydrophobic solutes (methane, ethane, krypton, and xenon) in water to study hydrophobicity at the most fundamental level. We present unequivocal experimental proof for strengthened water hydrogen bonds near purely hydrophobic solutes, matching those in ice and clathrates. The water molecules involved in the enhanced hydrogen bonds display extensive structural ordering resembling that in clathrates, thus indicating the fundamental interconnection between electrostatic screening (shielding) and the hydrophobic effect.
decoupled IR oxygen-deuterium (O-D) stretching mode (ν OD ) of four purely hydrophobic solutes in water together with ab initio molecular dynamics (MD) simulations of methane in water to study the origin of hydrophobicity at the most fundamental level.

Results and Discussion
IR Spectroscopy of Water Molecules in the Neighborhood of Purely Hydrophobic Solutes. The redshifts and line widths of the O-D stretching mode of the water molecules were measured when perturbed by small purely hydrophobic solutes: methane, ethane, krypton, and xenon. The very low solubility of these solutes in water was increased by using high pressure (14-56 bar) and low temperature (285-299 K). The schematic setup of the IR experiment is shown in SI Appendix, Fig. S1.
We focused on the decoupled O-D stretching mode (ν OD ) centered at ∼2,500 cm −1 (Fig. 1A and SI Appendix, Fig. S2 Because of the large differences between the frequencies of ν OH and ν OD , the intra-and intermolecular couplings are completely removed. The spectrum of ν OD of the water molecules perturbed by solute ( Fig. 1 B and C and SI Appendix, Fig. S4) was obtained by subtracting the spectra of pure H 2 O and pure, unperturbed HDO from the spectrum of the same solution containing a small amount of hydrophobic solute. This double-subtraction procedure, introduced by Lindgren and coworkers (28), was modified to account for variable cell thickness and tested thoroughly using pure solvent (2.8% HDO in H 2 O) and solutions of NaCl and methanol. The effects of temperature and pressure on the difference spectra have been investigated (SI Appendix). The spectra of H 2 O and HDO diluted in H 2 O used for subtractions were obtained at the same temperatures and pressures as the spectra of solutions of alkanes and noble gases. The same high-pressure transmission cell was used to obtain all of these spectra. To preclude formation of solid clathrates, we used pressures and temperatures positioned on the right-hand side of the equilibrium line in the pressure-temperature (PT)-phase diagram ( Fig. 1D and SI Appendix).
The resulting spectra of the water molecules that are perturbed by the solutes are shown in Fig. 1 B and C and summarized in Table  1. All four solutes display a redshift of ν OD (Δν OD ) of ∼60 cm −1 , the same as those of HDO ice and HDO clathrates ( Fig. 1C and Table  1). This observation shows that the strengths of the H bonds near purely hydrophobic solutes are enhanced to the level observed in ice and clathrates. Neither the temperature nor the pressure influence Δν OD significantly. The number of enhanced or ice-like H bonds per methane molecule is estimated to be between 10 and 15 (SI Appendix). The number of water molecules in hypothetically frozen structures is between five and eight per methane molecule. This amount is only a fraction of the number of waters in the first solvation layer of methane (i.e., 16) measured by neutron scattering (11). Our estimated number is in good agreement with the number obtained by Kauzmann (7). He estimated that, if icebergs around a methane molecule in solution are highly crystalline, they must contain less than half a dozen water molecules (7).
The line widths H w of ν OD of the water molecules perturbed by hydrophobic solutes are only slightly larger than those of clathrates and much smaller than those of the liquid solvent ( Fig. 1C and Table 1). Temperature and pressure do not influence H w significantly (by <6 cm −1 ) (  (32)(33)(34)(35)(36). Ab initio MD simulations currently provide the most accurate approach to addressing the structural and dynamic aspects of water-water H bonding in a condensed phase under ambient conditions, because in these methods, the interatomic forces are derived on the fly at the quantum mechanical level. Such simulations can explicitly capture electronic polarization effects and account for the proper response of water molecules to the local environment. We performed ab initio MD simulations of methane in D 2 O at three temperatures (283, 293, and 300 K) using the dispersion-corrected density functional theory (DFT) approach as implemented in the program package VASP, version 5.3 (37) (Materials and Methods and SI Appendix). The most appropriate choice of exchangecorrelation functional was found to be revPBE (38,39) in combination with the Grimme D3 correction (40), which was shown in a series of test simulations of bulk water according to the ability to best reproduce experimental radial distribution functions, diffusion constants, line shape of ν OD , and solvation enthalpy of liquid water at ambient conditions (SI Appendix). Time-dependent O-D stretching frequencies of individual water molecules were derived using the wavelet transform of the corresponding O-D distance time series (Materials and Methods).
Our analysis was focused on comparing the structural and energetic properties of the H-bond network of water molecules in the immediate vicinity of methane and relative to those of bulk water (Fig. 2). The ν OD frequencies as a function of distance to methane ( Fig. 2A) clearly show notable depressions in the region of the first hydration shell. The redshifts range from 14 to 7 cm −1 with increasing temperature. Note, however, that the redshifts obtained by the ab initio MD simulations are not directly comparable with the experimental redshift (∼60 cm −1 ) (SI Appendix). Furthermore, we have analyzed the electrostatic interaction energies of water molecule pairs located in either the methane hydration layer or the bulk water. The electrostatic interaction potentials of water molecule pairs located in the methane hydration layer are consistently lower than those in the bulk water, indicating higher H-bond strength (Fig. 2B). The difference between hydration and bulk water pairs (Fig. 2B) shows clearly that the first four nearest water molecules surrounding a water molecule located in the first hydration shell around the methane interact more strongly with the tagged water molecule (by 0.1-0.15 kcal/mol) than molecules in bulk water, whereas the difference for the more distant neighbors (fifth onward) becomes negligible. This result shows that the structure of water is tetrahedral (four neighbors). The number of H bonds per water molecule and the tetrahedrality (Fig. 2 C and D and SI Appendix), both as a function of distance to the methane, indicate that water H bonds in the hydration shell are more numerous and more Temperature T (Kelvin), pressure p (bar), frequency at maximum intensity ν OD (centimeters −1 ), frequency redshift Δν OD (centimeters −1 ), and half-width H w (centimeters −1 ) are shown. H w is defined as the width of a peak at onehalf height. Errors in measuring T and p are estimated at 0.5 K and 0.5 bars, respectively. All samples contain a solvent mixture of 2.8% HDO in H 2 O. *Solvent mixture (2.8% HDO in H 2 O) under four different conditions. † Spectral parameters of the water molecules perturbed by solutes (Fig. 1B). ‡ Solution of 0.03 M NaCl used for verification of the subtraction procedure (SI Appendix, Figs. S4 and S6). Note that the H w of its ν OD is equal to that of liquid solvent. § Solid state (SI Appendix, Fig. S3).  tetrahedrally oriented than those in bulk water. Other structural and energy parameters that characterize H-bond strength are consistent with the results presented above (SI Appendix, Fig. S12). The correlation times, calculated by the ab initio MD simulations, show restricted orientational mobility of the water molecules near a methane solute (SI Appendix, Fig. S13), in accord with the NMR data from the work by Haselmeier et al. (8). The O-D vibrational redshift, the tetrahedrality, and the number of H bonds are larger at lower temperatures. The simulations performed in this work clearly show that water molecules in the vicinity of methane have stronger, more numerous, and more tetrahedrally oriented H bonds than those in bulk water and that their mobility is restricted, which is consistent with the experimental results presented above.
Origin of Strengthened Water H Bonds near Hydrophobic Solutes. The results presented above show that the water H bonds in the neighborhood of purely hydrophobic solutes are as strong as those in ice or clathrates. These ice-like H bonds cause greater water ordering because of formation of icebergs as postulated by Kauzmann (7). Determining the physical origin of ice-like H bonds near hydrophobic solutes is, therefore, crucial for understanding hydrophobicity at the most fundamental level. Theoretical simulations show that water H bonds straddle small hydrophobic solutes in a way similar to that of H bonds in clathrates to maximize the number of H bonds (33). Such a constraint imposed on the H-bonding network can cause ordering of water molecules per se; however, it does not necessarily make H bonds stronger. The absence of a clear connection between steric constraint and H-bond strength led us to focus on electrostatic interactions, because water H bonds are predominantly electrostatic in nature (41). H bonds in water are known to be strongly cooperative (23,(42)(43)(44); their strength, therefore, depends on electrostatic interactions of an H-bonded pair with neighboring water molecules.
We, therefore, analyzed the effects of neighboring water molecules on the strength of the H bond between a pair of waters (Fig.  3A, yellow) during ab initio MD simulations. It has been shown that the instantaneous frequency, ν OH , is proportional to the projection of an electric field on an H atom along the bond vector O-H (E O-H ) (23,25,45). The electric field on an H atom depends on the positions, orientation, and point atomic charges of water molecules in the vicinity. To verify this relationship, we calculated the electric field on D atoms using the Hirshfeld's point atomic charges (46) and the corresponding frequencies of the O-D stretching mode ν OD for short segments of the ab initio MD trajectory using the continuous wavelet transformation method (47). A remarkably simple linear relationship is seen to exist between the ν OD frequency (i. The ν OD of an H-bonded water pair is affected predominantly by water molecules in the first solvation layer (cutoff < 5 Å) (Fig. 3B).
Three distinct classes of water molecule (Fig. 3B) were identified in the neighborhood of a water H-bond pair in liquid water. The first two classes comprise the water molecules that are H-bonded to the H-bonded pair; their contributions to the value of E D. . .O do not differ significantly from the contributions of the corresponding waters in ice. In ice, a pair of H-bonded water molecules is H-bonded to six water molecules located at the vertices of two fused tetrahedrons (Fig. 3C). In liquid water, the locations of H-bonded waters resemble fused tetrahedrons (Fig. 3D); however, some vertices are empty or overpopulated. The contributions of the first two classes of water molecules to the value of E D. . .O are determined primarily by the point atomic charges of atoms that are closest to the atom D of the H-bonded pair controlled by donor or acceptor characters of the corresponding H bonds. The first class is the H bond-enhancing water molecules, which increase the value of E D. . .O (Fig. 3B, red). The second class is the H bond-weakening waters, which decrease E D. . .O (Fig. 3B, blue). The average numbers of these two classes per H bond are 3.4 and 1.7; therefore, the net effect is to strengthen the H bond. In ice, the corresponding numbers are four and two.
The most interesting is the third class of intercalating water molecules. They are able to occupy transiently the space closest to the atom D of the H-bonded pair, which is empty in ice (Fig. 3B,  green). To accommodate the intercalating water molecule, the angle O-D . . . O of the H-bonded pair must bend (Fig. 3D). Only ∼6% of the H-bonded pairs in the bulk water (at cutoff of 2.65 Å) have intercalating waters. The intercalating waters represent a subgroup of the so-called fifth, interstitial, nontetrahedral, or mismatched water molecules (29,(48)(49)(50) that interact with the tagged water molecule by, predominantly, van der Waals forces. These water molecules are specific for liquid water and absent in ice. The intercalating water molecules are not directly H-bonded to the H-bonded pair; therefore, they have less negative solvation energy than the bulk water. The orientation and thus, also, the contribution of these water molecules to the value of E D. . .O are determined mainly by the electric field generated by the first two classes. Fig.  3D shows that the average dipole moments of the intercalating waters are oriented in the direction of the average local electric field (Fig. 3E). The intercalating waters screen (shield) this electric field. Consequently, the contributions of the intercalating water molecules to E D. . .O are generally negative (Fig. 3B, green) and primarily responsible for the lower strength of H bonds in liquid water than those in ice.
In the transfer of the H-bonded pair from the bulk to the first solvation layer of methane (Fig. 3E), the fraction of intercalating water molecules per H-bonded pair is reduced from 6 to 2%, because methane and intercalating water molecules tend to occupy the same space lateral to the direction of H bonds (Fig. 3 D and E). Methane prefers this position, because water H bonds tend to straddle small hydrophobic solutes in a way similar to that in clathrates. By removing the intercalating water molecules, the value of E D. . .O increases. A higher value of the projection of electric field E D. . .O means that a larger force is pushing a positively charged D atom toward the acceptor O atom, which strengthens the H bond. This mechanism explains why water H bonds near methane are stronger than those in bulk water.
The explanation described above for the strengthening of water H bonds near hydrophobic solutes is in accord with the general concept that electrostatic interactions are stronger near hydrophobic groups. Berry and coworkers (51) have shown that water molecules in the neighborhood of a hydrophobic solute have a smaller dielectric susceptibility or in microscopic terms, less effective electrostatic screening (52). It has been shown that electrostatic screening is the main reason for the distinct backbone conformational preferences of amino acid residues in peptides and proteins (53), the nearest neighbor effect (54), and the formation of transient β-strands in unfolded proteins (55). It has also been shown, using amide to ester mutations, that H bonds in the protein interior are stronger than those exposed to solvent (56)(57)(58).

Conclusion
By applying methods of superior accuracy for measuring and calculating subtle effects of H-bonding in water, we have obtained evidence that supports, unequivocally, the iceberg view of hydrophobicity proposed by Frank and Evans (6) and Kauzmann (7). Our experimental results show that water H bonds near purely hydrophobic solutes are strengthened to the level observed in ice and clathrates. We have proposed a physical origin for hydrophobicity based on electrostatic screening that couples hydrophobic with electrostatic interactions. This coupling may be crucial in understanding protein folding and other complex phenomena.
Instruments and Sample Preparation. IR spectra were recorded on Bruker Vertex 80 and Tensor 27 Spectrometers and collected in the transmission mode with a nominal resolution of 4 cm −1 . Typically, 128 interferograms were averaged and apodized using the Gapp-Henzel function. A deuterated triglycine sulphate (DTGS) detector was used throughout. The aperture of the IR beam was set to values between 2 and 4 mm. A Harrick (HPL-TC-13-3) high-pressure and hightemperature sample cell was used for high-pressure measurements (SI Appendix, Fig. S1). The errors in measuring temperature and pressure are estimated at 0.5 K and 0.5 bar, respectively.
OPUS software (Ver. 6.5) was used for subtraction of spectra. Elimination of CO 2 bands has been performed manually.
Solutions of 2.8% (by volume) HDO in H 2 O were used for all IR measurements. They were prepared by mixing 1.4% D 2 O and 98.6% H 2 O. HDO/H 2 O solution was used to eliminate inter-and intramolecular vibrational coupling and prevent saturation of the O-D stretching mode by significantly reducing its intensity. The small amount of HDO in H 2 O ensures that intensities of the O-D stretching peak are below saturation level. The solution of 2.8% HDO in H 2 O is a compromise between factors that govern band saturation (i.e., dynamic range of the fully decoupled O-D stretching mode and subtraction of the H 2 O combination band). Spectra of pure H 2 O were recorded at the same temperatures and pressures as those for the solutions.
Ab Initio MD Simulations. The ab initio MD simulations were performed using the Vienna Ab Initio Simulation Package (VASP, version 5.3) (37). This code, one of the fastest ab initio computer packages available, combines plane-wave implementation of the DFT in combination with the projector-augmented wave (PAW) pseudopotential (59). In general, the PAW pseudopotential is more accurate than the ultrasoft pseudopotential, because it provides an exact valence wave function in the core region of the electron orbital.
The molecular systems were constructed at three different temperatures (283, 293, and 300 K) by inserting one methane molecule per cubic box containing 91 D 2 O molecules. Sizes of simulation cells were set according to the experimental densities of liquid water at three corresponding temperatures, assuming that the volume of one methane molecule occupies approximately the volume of two water molecules. Accordingly, three cubic boxes of sizes 14.06605, 14.07299, and 14.08090 Å were used for the simulations at 283, 293, and 300 K, respectively. The initial molecular setup was achieved by running short classical MD simulations in the microcanonical (NVE) ensembles using the TIP4P/2005 water model (60) with methane parameters introduced by Docherty et al. (61). This step was followed by extensive ab initio MD simulations, in which each system was first equilibrated for 100 ps and simulated further for 250 ps at 283 and 293 K and 100 ps at 300 K to produce corresponding canonical (NVT) ensembles for analysis. Shorter (100 ps) control MD runs were performed at slightly increased and decreased densities by varying the size of the simulation box by ±1% to explore the effect of pressure on the results.
Simulated trajectories were saved every 1 fs and used for structural and O-D stretching vibrational analyses. The O-D stretching vibrational spectrum was calculated based on continuous wavelet transformation of the trajectories as introduced by Vela-Arevalo and Wiggins (47). The vibrational frequency of an individual O-D bond at a given moment is obtained through a time series analysis of the O-D distance. This approach enabled position-and timedependent O-D frequencies to be obtained, which could be correlated with the structural and dynamic properties of the H-bonded network. The point atomic charges used to calculate electric fields were calculated from the instant electron density of the system according to the Hirshfeld definition (46).