Nonpolar solutes enhance water structure within hydration shells while reducing interactions between them

MD Simulations. MD simulations were computed by using the gromacs program package as previously described (18). Each simulation consisted of a single solute, either benzene or the chair conformation of cyclohexane, in a 40-Å cubic box filled with TIP4P (four-point transferable intermolecular potential) water (31). The OPLS (optimized potentials for liquid simulation) force field (32) was used to model the solutes. Ten simulations, each 1 ns in length, were run at 300 K at constant temperature and pressure (33). The Particle Mesh Ewald summation (34) was used for the electrostatic interactions, and a cutoff of 10 Å was used for the van der Waals interactions. A 0.002-ps time step was used, and the coordinates were output every 2 ps. The first 500 ps of each simulation were not used in the subsequent analysis. A set of pure TIP4P water simulations was also run under identical conditions.

Spatial Distribution of Water Position and Orientation. Previously, we computed the SDFs of the positions of the water molecules surrounding benzene and cyclohexane (18). Unlike the radial distribution function, g(r), the SDF is resolved in Cartesian space and not rotationally averaged around the solute. The simulation box was rotated onto reference axes centered on the solute, and the positions of the water O atom centers were stored in 3D tables (20-Å cube composed of 0.25-Å cube volume elements). The O density, ρ_O, was then normalized relative to the density of a water-only simulation (ρ = 1.0). In the current study, the average orientation of these water molecules as a function of position was also measured; we refer to these functions as SOFs. The orientation of a water molecule is described by its orientation vector, \begin{equation*}{\vec{{\mu}}}\end{equation*} , which is the vector of unit length originating at the O atom center and bisecting the two H atoms (see Fig. 2 A). The individual components of \begin{equation*}{\vec{{\mu}}}\end{equation*} were summed and stored in a 4D table, where the first three dimensions correspond to the coordinates of the water O position as in the SDF, and the fourth dimension was used to sum the x, y, and z components of \begin{equation*}{\vec{{\mu}}}\end{equation*} . After all of the MD frames were analyzed, the sums of the individual components for each SOF voxel were divided by the total number of water O counts in the corresponding SDF voxel. The resulting components define the average water orientation vector, \begin{equation*}{\langle}{\vec{{\mu}}}{\rangle}\end{equation*} , for each volume element in the SOF.

Radial/Angular Distribution Functions. Solute centroid to water O radial distribution functions, g_Sol-O(R, cos ω), were computed from the MD simulations, but were expanded in a second dimension that corresponds to the tilt angle, ω. The tilt angle was defined as the angle between the water orientation vector, \begin{equation*}{\vec{{\mu}}}\end{equation*} , and the radial vector, \begin{equation*}{\vec{R}}\end{equation*} , defined as the vector between the water O center and the solute centroid (see Fig. 2 A). The value of cos ω was used because this quantity is expected to be uniform for an isotropic distribution of orientations. The distance bin width was 0.05 Å, and the angle bin width was 1°. The density in the radial (R) dimension was normalized in the usual fashion. Because of the oblate shapes of both benzene and cyclohexane, g_Sol-O(R,cos ω) was computed in two volume regions, the axial region and the equatorial region (ref. 18; see Fig. 2B).

Hydration Shells Defined by Distance Cutoff. To study the H-bonding structure within and between the hydration shells surrounding the solute molecules, water molecules corresponding to each shell were identified in an explicit manner. For the DC method, water molecules were classified as S1 if their O centers were <5.4 Å away from the nearest solute C atom. S2 water molecules were those that fell within 5.4–8.9 Å of the nearest solute C atom. All other water molecules were considered to be bulk (B). The 5.4-Å cutoff value was taken from the first minimum in g_Sol-O(r) from a methane–water simulation (data not shown), and the 8.9-Å cutoff was set to yield the same number of water atoms in the second shell as obtained by the VP method. For the pure TIP4P water control, cutoffs of 4.9 Å for S1 and 8.4 Å for S2 were used to match the numbers of S1 and S2 waters determined by VP.

Hydration Shells Defined by VP. The VP surrounding the heavy atoms (O and C) in the simulation box were computed by using the contact area program implemented by Patrice Koehl (35). A contact is defined as two atoms that share a polyhedron face. Based on the VP-generated contact lists, a water molecule was labeled as S1 if its O contacted at least one solute C atom. A water molecule was labeled as S2 if it contacted at least one S1 water molecule and did not contact the solute. All other water molecules were considered to be B. This geometric method differs notably from the DC method in that it depends on the instantaneous arrangement of all of the atoms in the simulation. For instance, two water molecules positioned the same distance from the solute can have differing shell labels, depending on the proximity of their neighbors. With the VP definition, S1 and B waters will never be in contact.

Measurements of Water–Water Contact Distance and Angle. Rather than define a strict set of geometric criteria for defining a H-bond, we use the method of Sharp (17, 24) and compute distributions of water O···O contact length and O···O–H angle within and between adjacent shells defined by either the DC or the VP method. The O···O–H angle, θ, is defined as the smallest of the four possible angles for each water pair (see Fig. 4B). This procedure is equivalent to finding the deviation from linearity of the best H-bond between the two water molecules. The distance and angle bin widths were 0.05 Å and 1°, respectively. These g_OO(r, θ) distributions were computed for five classes of water contacts: S1–S1, S1–S2, S2–S2, S2–B, and B–B. For the VP-based shells, the distributions were computed only for those water pairs found to be in contact. In this case, the distributions are reported as probabilities of observing particular geometries within the classes of interactions, ρ_OO(r, θ).

Normalization of DC-Based Water–Water Contact Distance and Angle Distributions. For the DC-based shells, the distribution of distances and angles for all labeled waters was computed. The distribution was normalized in the contact length dimension to give g_OO(r). For molecules that are spherically symmetric, g(r) is normalized by assuming uniform bulk density in each spherical shell. For an asymmetric solute such as benzene or cyclohexane, such an analytical approach cannot be used. Thus, we generated an expected distance distribution from a random sampling (36) of points distributed within (and between) the shell volumes observed in the MD simulations (see Computational Methods, which is published as supporting information on the PNAS web site).

SOFs of Benzene and Cyclohexane. The SDFs of water molecules surrounding benzene and cyclohexane were previously published (18). To complement that study, we computed the SOFs for the water molecules surrounding benzene and cyclohexane. Fig. 1 shows slices through the SOFs for both solutes. Projections of the longest \begin{equation*}{\langle}{\vec{{\mu}}}{\rangle}\end{equation*} vectors into the respective planes are drawn as arrows. The most striking feature of the benzene SOF (Fig. 1 A) is the strong orientation of the water molecules located along the z-axis of the benzene; they have \begin{equation*}{\langle}{\vec{{\mu}}}{\rangle}\end{equation*} vectors that point directly toward the benzene centroid. This finding is consistent with our previous SDF results and the results of other simulations (19–23) and experiments (37, 38) that have shown a preference for H-bonding to the faces of the benzene ring. The axial ordering extends beyond the first hydration shell, with longer arrows seen further out along the z-axis, corresponding to second shell water molecules. In the equatorial region around benzene (Fig. 1C), the ordering is less pronounced. The \begin{equation*}{\langle}{\vec{{\mu}}}{\rangle}\end{equation*} vectors of the waters closest to the edge of benzene tend to point away from the solute, influenced by the localized positive charge in that region.

The SOF of the water molecules around cyclohexane is quite different from that of benzene. In Fig. 1 B and D, the water molecules in the first hydration shell show less ordered \begin{equation*}{\langle}{\vec{{\mu}}}{\rangle}\end{equation*} vectors than those surrounding benzene. The water molecules closest to the cyclohexane show a slight tendency to point outwards, but the majority of the volume occupied by the first hydration shell is devoid of any long \begin{equation*}{\langle}{\vec{{\mu}}}{\rangle}\end{equation*} . This finding indicates that ordering in the first hydration shell around a hydrophobic solute does not occur by forcing the waters into a fixed orientation; instead water molecules that occupy the same location in space sample a variety of orientations.

Radial/Angular Distribution Functions for Benzene and Cyclohexane. To further explore the orientation preferences in the first hydration shell, we measured a combined solute–water radial/angular distribution function, g_Sol-O(R, cos ω). Fig. 2 shows the g_Sol-O(R, cos ω) for the axial and equatorial volume regions surrounding benzene and cyclohexane. Fig. 2C confirms that the water molecules axial to benzene occur at high density and with a strong inward orientation. The maximum of cos ω in this region is 0.66, corresponding to ω = 48.7°. This value is very close to the expected value of 52.26° for orienting an O–H bond vector directly along \begin{equation*}{\vec{R}}\end{equation*} (the H–O–H bond angle for TIP4P is 104.52°) (Fig. 2G). For the first hydration shell waters that are either equatorial to benzene (Fig. 2D) or axial or equatorial to cyclohexane (Fig. 2 E and F), the orientation preference is less strong, but the maximum of cos ω occurs between –0.46 to –0.50, corresponding to an ω angle of 117–120°. This angle positions an O–H bond vector roughly tangential to the surface of the solute, allowing for formation of a H-bond network within the first hydration shell (Fig. 2H).

Identification of Hydration Shells by Voronoi Contacts and Cutoff Distances. To test the assumption that the tangential orientation preference of the water molecules in the first hydration shell surrounding a hydrophobic solute allows for enhanced H-bonding, we decided to rigorously quantify the first and second hydration shells surrounding benzene and measure the H-bonding structure within and between the shells. For the DC method, the average number of water molecules in S1 for benzene was 31, and in S2 it was 98. The average number of hydration shell waters determined by VP was similar to DC, with 33 waters in S1 and 99 in S2 for benzene. An illustration of the shell assignments for a single MD frame is shown in Fig. 3. While there is significant overlap in the shell assignments between the VP and DC methods, some differences in labeling are observed.

Water–Water Contact Distance and Angle Distributions for Interactions Within and Between Hydration Shells. The distributions of water pair O···O distances (r) and minimum O···O–H angles (θ) were measured for all inter- and intrashell contacts. For the DC method, all water molecules positioned within 6 Å of each other were considered in contact. The g_OO(r, θ) distributions for the five classes of water interactions found by using the DC method of shell labeling on the benzene simulations are shown in Fig. 4 A–E. Each graph has a characteristic shape, with a tall peak at short r and low θ (maximum at 2.7 Å and ≈9°), and a shorter, wide peak at longer r and a broader distribution of angles (maximum at 3.9–4.3 Å and 42.5–49.5°). The first peak corresponds to H-bonded water pairs, and the second peak corresponds to contacting waters that are not H-bonded to one another but are instead both H-bonded to a third water molecule (Fig. 5). The largest differences are seen between the S1–S1 distribution and the other classes of contacts. Within S1, the population of water molecules that fall within the H-bonding geometry is larger than that within the other classes. The H-bonding peak is lowest in the S1–S2 distribution, but the remaining three classes, S2–S2, S2–B, and B–B, are nearly identical to each other.

The distributions of water O···O distances and minimum O···O–H angles, ρ(r, θ), for the VP method applied to the benzene simulations are shown in Fig. 4 F–J. The VP distributions have a shape similar to that of the DC distributions. The primary peak has a maximum at r = 2.8 Å and θ = 8.5–10.5°, and the secondary peak is at r = ≈4.2 Å and θ = 42.5–48.5°. The relative height of the second peak is more pronounced because it is not flattened by normalization as is the case in a typical g(r). The shapes of the boundaries of the distributions at longer contact length (r > 5.4 Å) are determined by the implementation of the Voronoi construction used (35), which will not allow any two water molecules separated by a distance that would accommodate an intervening water molecule to be in contact.

The water structure distributions for the VP-based shells show structural differences between the five classes of interactions. The S1–S1 interactions show the most H-bonding structure, with the largest difference between the first and second peak heights. Likewise, the S2–S2 interactions show enhanced H-bonding, relative to the B–B interactions. Interestingly, the interactions between S1 and S2 show reduced H-bonding, with the majority of the population falling into the second peak. A similar pattern is seen in S2–B. These results show enhanced H-bonding structure within the first and second hydration shells, accompanied by decreased structure between them. While the DC method shows that the interactions within S1 are enhanced, the reduced structure in S1–S2 is observable, but more subtle. The results obtained for the cyclohexane simulations (not shown) were practically identical to those of benzene.

Fig. 6 shows for benzene the 2D water structure distributions from Fig. 4 projected into either the r or θ dimension. The results for cyclohexane (not shown) are indistinguishable from those of benzene. While both the VP and DC methods show an increase in H-bonds within S1 coupled with a reduction of H-bonds between S1 and S2, the changes are much more marked for the VP method (Fig. 6 A and B). In addition, the VP method shows an increase in H-bonding within S2 coupled with a reduction of H-bonds between S2 and B. Similar trends are seen for the angular distributions, with the height of the first peak measuring the extent of H-bonding (Fig. 6 G and H). The assignment of shell labels by using either the DC or VP method selects for water molecules that fall within geometrically defined volumes. The results shown in Fig. 4, therefore, may be influenced by these geometric restraints. To better understand these effects, we applied our analysis methods to a control simulation of pure TIP4P water molecules, specifying a single water molecule as the “solute.” The resulting measures of shell interaction geometries (Fig. 6 C and D) show that there are structural differences within and between the hydration shells of a single water molecule defined by both methods. Ratios of first and second peak heights for both solutes are given in Table 1, which is published as supporting information on the PNAS web site. Both methods show an apparent increase in the tangential H bonds within S1. In addition, the VP method shows similar results within S2. Both methods also show an apparent decrease in radial H-bonds (between S1 and S2). At present, we do not understand the cause of this geometric effect. We can, however normalize for this bias by dividing the benzene results by those for pure water (Fig. 6 E and F). These normalized distributions still show that S1···S1 hydrogen bonds are increased, whereas the S1···S2 hydrogen bonds are decreased.

The origins of the hydrophobic effect are intimately linked with the structural changes induced by a nonpolar solute in water molecules that surround it. In this paper, we have measured the effects of the hydrophobic solutes benzene and cyclohexane on the orientation preferences of the hydration shell water molecules. These results not only highlight the obvious ordering of water molecules H-bonded to benzene but also reveal subtle biases in orientation in the first hydration shell associated with the hydrophobic effect. A key observation of our work is enhanced H-bonding structure within the first hydration shells combined with reduced hydrogen bonding with the second shell.

The SOFs show that the \begin{equation*}{\langle}{\vec{{\mu}}}{\rangle}\end{equation*} of the water molecules in the axial regions of benzene are strongly oriented inward, providing confirmation of H-bonding with the faces of the ring. The two H-bonded waters are only a small subset of the ≈30 water molecules found in the first hydration shell. The water molecules closest to the edge of benzene have a weaker, yet detectable, bias for an outward orientation. This result makes sense in the context of the localized partial positive charge on the edge of the ring. This ordering is not seen evenly throughout the majority of the volume of S1, however, which has relatively few long \begin{equation*}{\langle}{\vec{{\mu}}}{\rangle}\end{equation*} vectors.

The SOF results for cyclohexane are interesting, showing less order than benzene, and a S1 that is not well populated with long \begin{equation*}{\langle}{\vec{{\mu}}}{\rangle}\end{equation*} vectors. There is no strong absolute orientation bias induced in S1 by cyclohexane, and the water molecules in this region maintain significant rotational degrees of freedom.

Whereas the SOF shows relatively little orientation structure, g_Sol-O(R,cos ω) is somewhat more sensitive to the structural biases of water near nonpolar solutes. In the regions around cyclohexane and equatorial to benzene, the water molecules sample a broad range of orientations, with the peak of the distribution corresponding to ω = 120°. This tilt angle would allow one O–H bond to orient roughly tangential to the solute surface, accommodating H-bonding with other S1 neighbors. Water molecules in the first hydration shell surrounding nonpolar surfaces, therefore, show a subtle but significant outward orientation preference.

Measuring H-bonding interactions within and between the hydration shells probes the orientation structure of water. Both the DC and VP methods show clear enhancement of the H-bonded population of neighboring waters in S1 and a reduction between S1 and S2. Thus, the strong H-bond network in S1 does not link the first and second shells.

The use of VP to identify both the shell labels and contacts of neighboring water molecules enhances the structural differences within and between the shells. The controls on pure water show this is a geometrical effect, which also occurs in a pure water control calculation. Our results emphasize the crucial role of controls in simulations and when using complicated geometrical procedures.

Author contributions: T.M.R. and M.L. designed research; T.M.R. performed research; T.M.R. and M.L. analyzed data; and T.M.R. and M.L. wrote the paper.

This paper was submitted directly (Track II) to the PNAS office.

Abbreviations: MD, molecular dynamics; DC, distance cutoff; VP, Voronoi polyhedra; SDF, spatial distribution function; SOF, spatial orientation function; S1, shell 1; S2, shell 2; B, bulk.

We thank Patrice Koehl for providing his Bestfit and Contact Area subroutines and Rachel Kolodny for helpful discussions. This work was supported by National Institutes of Health Grant GM-41455. M.L. was supported by the Fondation de l'Ecole Normale Supérieure Chaires de Blaise Pascal.