Size-Dependent Order-Disorder Crossover in Hydrophobic Hydration: Comparison between Spherical Solutes and Linear Alcohols.

Theory and computer simulation studies have predicted that water molecules around hydrophobic molecules should undergo an order-disorder transition with increasing solute size around a 1 nm length scale. Some theories predict the formation of a clathrate-like ordered structure around smaller hydrophobic solutes (<1 nm) and the formation of disordered vapor-liquid interfaces around larger solutes (>1 nm) and surfaces. Experimental validation of these predictions has often been elusive and contradictory. High-resolution Raman spectroscopy has detected that water around small hydrophobic solutes shows a signature similar to that of bulk water at lower temperature (increased ordering and a stronger hydrogen-bonded network). Similarly, water around larger solutes shows an increasing population of dangling OH bonds very similar to higher temperature bulk water. Thus, the solute size dependence of the structure and dynamics of water around hydrophobic molecules seems to have an analogy with the temperature dependence in bulk water. In this work, using atomistic classical molecular dynamics (MD) simulations, we have systematically investigated this aspect and characterized this interesting analogy. Structural order parameters including the tetrahedral order parameter (Q), hydrogen bond distribution, and vibrational power spectrum highlight this similarity. However, in contrast to the experimental observations, we do not observe any length-dependent crossover for linear hydrophobic alcohols (n-alkanols) using classical MD simulations. This is in agreement with earlier findings that linear alkane chains do not demonstrate the length-dependent order-disorder transition due to the presence of a sub-nanometer length scale along the cross section of the chain. Moreover, the collapsed state of linear hydrocarbon chains is not significantly populated for smaller chains (number of carbons below 20). In the context of our computational results, we raise several pertinent questions related to the sensitivity of various structural and dynamical parameters toward capturing these complex phenomena of hydrophobic hydration.

Introduction

Hydrophobic interactions facilitate a wide range of molecular phenomena in aqueous solutions such as protein folding, formation of self-assembled structures such as membranes and micelles, protein–ligand interactions, and so on.[1−6] Hydrophobic hydration and interaction are some of the most widely studied, debated, and yet to be clearly understood aspects of water-mediated interactions. In a series of seminal works, Chandler and co-workers have proposed a nanometer-scale order–disorder transition in the hydration shell of hydrophobic solutes.[7−10] While a small hydrophobic molecule such as methane can just about fit in the cavities of the hydrogen-bonded network of bulk water, it leads to a significant reduction in entropy. While some earlier theories associated this observation with the formation of ordered clathrate-like structures around the solute, several simulation and experimental studies could not detect the formation of such structures.[11] However, theories based on excluded volume effects could explain the reduction in entropy without invoking the formation of clathrate-like structures.[12] On the other hand, for larger solutes (and surfaces), water cannot maintain the hydrogen-bonded network and leads to a more disordered (possibly vapor-like) layer around the solute. Experimental observations of such a vapor-like layer have been elusive due to the transient nature of the density fluctuations at the interface. Moreover, the theoretical framework has been built over models with purely repulsive interactions between the solute and the solvent, whereas a little bit of attractive interaction may drastically alter the scenario for real molecules.

Over the last decade, Ben-Amotz and co-workers have performed a series of pioneering experiments using temperature-dependent Raman scattering measurements with multivariate curve resolution to capture the structure of water in the hydration shell of hydrophobic solutes with a very high degree of resolution.[13−15] They provided one of the first direct experimental pieces of evidence of the size-dependent order–disorder transition discussed above. Based on the observed shift in the OH stretch band, they concluded that the water around small solutes behaves similar to cold water (more structured) and those around larger solutes behaved as hot water (less structured). The population of the dangling OH bonds (non-hydrogen-bonded) increases as the solute size increases.[14] On the other hand, several scattering experiments were unable to detect any enhanced structural ordering around small hydrophobic solutes.[16,17]

Even earlier simulation works have reported apparently contradictory results. Ab initio molecular dynamics (MD) simulations have reported evidence of dangling OH bonds around hydrophobic solutes, but a significant increase in tetrahedral ordering has not been observed.[18,19] On the other hand, evidence of the presence of both “cold” and “hot” water around methane has been presented using computational analyses depending on which way the water is facing around methane.[20]

In recent times, several theoretical and experimental studies have explored the nature of hydrophobic hydration near biological and electrochemical interfaces.[21−24] It has been highlighted that the presence of such interfaces can significantly modulate the solvation thermodynamics for solutes present near the surface. Characterization of hydrophobicity near a biomolecular surface becomes particularly challenging due to the presence of both chemical and topographical heterogeneities. In a seminal work, Patel and co-workers have demonstrated that protein patches containing significant amounts of hydrophilic/polar content can still function as hydrophobic patches.[23] Thus, in the context of molecular recognition, it is extremely crucial to understand the nature of hydrophobic hydration in realistic molecular systems with chemical and structural heterogeneities since such systems often present us with counterintuitive emergent phenomena in contrast to simpler model systems such as spherical solutes.

Given the above considerations, we find it important to systematically investigate the nature of the structural and dynamical transition of a hydration shell of hydrophobic solutes with varying sizes. We also highlight the interesting analogy between the solute size dependence and the temperature dependence in bulk water. In addition to model spherical solutes, where the transition is more prominent, we also take up alcohol molecules with varying chain lengths n (CH2OH) since this is the model system used in most experimental investigations owing to the poor solubility of pure hydrocarbons in water. Our earlier work has demonstrated clearly that it is unlikely to observe a length-dependent crossover in linear hydrocarbon chains due to the presence of a shorter length scale along the cross section of the polymer.[11] Hence, we do not expect a size-dependent crossover at least for the short alcohol chains traditionally used in the experimental studies. Here also, we demonstrate that linear alcohols do not exhibit any dramatic length-dependent crossover unlike their spherical counterpart. We suggest the possibility that while the OH stretch can be a sensitive reporter of the local environment and interactions encountered by the water molecules, it does not necessarily correlate with the nanometer-scale crossover encountered for ideal spherical hydrophobic solutes.

Methodology

In this study, our objective is to examine the behavior of the hydration shell water around hydrophobic solutes of varying sizes/lengths in comparison to the temperature variation in bulk water. For this purpose, we have used two types of solutes: (i) spherical solutes with varying sizes and (ii) linear alcohol molecules with varying chain lengths (CH2OH). All atomistic classical MD simulations have been carried out using the GROMACS software.[25] Two types of water models have been used: (i) rigid TIP4P/2005[26] to look into the structural properties of water and (ii) flexible TIP4P/2005f[27] to study the vibrational power spectrum.

For the spherical solutes, we have employed our previously studied models of hydrophobic solutes,[11] namely, “single LJ” and “multi LJ”. The “single LJ” model consists of a single Lennard-Jones particle with a varying size parameter (σ), whereas the “multi LJ” model consists of a droplet of multiple methane-like particles into a spherical packing. Both the spherical hydrophobic solutes and alcohols have been modeled using united-atom OPLS-UA and all-atom OPLS-AA force fields, respectively.[28] The size of hydrophobic solutes is varied from 0.4 to 3 nm (single LJ from 0.4 to 1 nm and multi LJ from 1 to 3 nm). The simulations of hydrophobic solutes and alcohols in water were performed at a temperature of 300 K. For the bulk water simulation, the temperature was varied from 273 to 373 K with an increment of 10 K. The concentration of alcohols was fixed at 0.5 M for methanol, 1-propanol (as used by Davis et al.[13] and Perera et al.[14]), and 1-heptanol. In addition, we have also simulated a single alcohol molecule in water with the effective concentrations as listed in Table S1 (Supporting Information).

The energy minimization of all the systems was done using the steepest-descent algorithm. The energy-minimized structures were then subjected to NVT equilibration for 5 ns at the specified temperatures for the system under consideration using a V-rescale thermostat,[29] subsequently followed by NPT equilibration at 1 bar pressure for 5 ns using a Berendsen barostat.[30] We have used a Parrinello-Rahman barostat[31] for the production runs. Periodic boundary conditions have been applied in all directions. In all the simulations, the long-range electrostatics have been solved using the particle mesh Ewald method.[32] Newton’s equations of motion have been solved using a leapfrog integrator with integration time steps of 0.1 and 2 fs for the flexible and rigid water models, respectively. The frames of the trajectories have been saved at a frequency of 1 fs (for a 1–5 ns total run length) and 2 ps (for a 50 ns total run length) for the flexible and rigid water models, respectively.

Results and Discussion

Size Variation of Hydrophobic Solutes versus Temperature Variation in Bulk Water

The hydrophobic solutes of different sizes are simulated in water at 300 K, while bulk water is simulated at temperatures ranging from 273 to 373 K with an increment of 10 K. In order to characterize the variation in structural order with size and temperature, we have computed the probability distribution of two structural order parameters: (i) tetrahedral order parameter (Q) and (ii) number of hydrogen bonds (H-bond) formed by each water molecule. Moreover, the vibrational power spectrum has been computed using the analogous flexible water model. Here, we compare these distributions for the hydration shell of solutes with increasing size and for bulk water with increasing temperature. Interestingly, Laage and co-workers have highlighted that a single parameter is not enough to describe the complexity in structure and dynamics of water, rather a combination of several order parameters would be crucial.[33]

(i) Tetrahedral Order Parameter (Q)

Tetrahedral order parameter (Q) as defined by Errington and Debendetti[34] provides a simple yet elegant way to quantify the “tetrahedrality” of the water molecules around a central water molecule. The expression is given bywhere Q is the tetrahedral order parameter of the ith water molecule and θ is the angle subtended on the oxygen atom of that water molecule by each pair (given by the indices j and k) of four nearest-neighbor molecules. Hence, for a perfect tetrahedral arrangement of the four neighboring water molecules around the central ith water molecule, Q = 1, whereas for a random and uniform distribution of these angles, Q = 0.

We expect higher tetrahedral ordering in water at lower temperatures, which would gradually decrease upon increase in temperature, as shown in Figure a for bulk water with temperature variations. The peak at around Q ∼ 0.85 is dominant at a lower temperature (T = 273 K), whereas with increasing temperature, the peak at around Q ∼ 0.5 takes over (T = 373 K). Interestingly, a nearly bimodal distribution is visible at around T ∼ 350 K, which indicates a weakly first-order order–disorder transition around this temperature.

Figure 1

Distribution of the tetrahedral order parameter (Q) for (a) bulk water with a temperature variation from 273 to 373 K and (b) first hydration shell water around the hydrophobic molecules for single and multi LJ models with the size ranging from 0.4 to 3 nm, for a 2D surface of a hydrophobic molecule and bulk water at 300 K. The size-dependent data presented in (b) have been taken from our earlier work.[11]

Similarly, the distribution of Q for the hydration shell water around “single LJ” and “multi LJ” with a varying size is shown in Figure b. The calculation of Q of the reference water molecule in this case includes the solute, where the need of including the solute in the neighbor list of reference molecules has been explained in our earlier work.[11] We observe that Q gradually shifts to lower Q (∼0.5) values with an increase in size of the hydrophobic solute. Although there is a rapid change in the distribution above ∼1 nm, the peak at Q ∼ 0.5 does not dominate even for largest spheres with a size of 2–3 nm. Interestingly, there is a dramatic change in the nature of the distribution for the 2D planar interface, where the population of the disordered water molecules at Q ∼ 0.5 is significant. This can be attributed to the difference in effective interaction between the solute and water since the packing density of particles is different between the largest multi LJ sphere and the 2D surface. The dependence of the average tetrahedral order parameter (⟨Q⟩) on the temperature and size has been shown in Figure S6. While ⟨Q⟩ decreases gradually with increase in temperature in bulk water, there is a sharp decrease with increasing size at around ∼1 nm.

(ii) Number of Hydrogen Bonds

We have further looked at the distribution of the number of hydrogen bonds formed per water molecule. A hydrogen bond between two water molecules is identified by the following conditions: (i) the O–O distance between the donor and acceptor molecules should be less than 0.35 nm and (ii) the H–O (donor)–O (acceptor) angle should be less than 30 degrees. The distribution of the number of hydrogen bonds (NHB) formed by the water molecules is presented in Figure ; (a) is for the bulk water with variations in temperature and (b) is of the hydration shell water around hydrophobic solutes with increasing size. The results of bulk water indicate that as the temperature increases, the capacity of forming four hydrogen bonds decreases gradually, and the probability of forming lesser (one, two, and three) hydrogen bonds increases. A similar trend in the distribution can be seen for the hydration shell water around hydrophobic solutes of varying size as expected. The difference of the distributions between the largest sphere (3 nm) and the 2D planar surface is quite visible even in this case, indicating that the specific topography of the surface plays an important role in dictating the local structuring of the hydration layer. The comparison between the temperature and size variation of the average number of H-bonds per water molecule has been shown in Figure S7, which has a behavior similar to that of ⟨Q⟩.

Figure 2

Distribution of the number of hydrogen bonds per water (NHB) in (a) bulk water with temperature variation from 273 to 373 K and (b) first hydration shell water around the hydrophobic molecules for single and multi LJ models with the size ranging from 0.4 to 3 nm, for the 2D hydrophobic surface and bulk water at 300 K. The size-dependent data presented in Figure b have been taken from our earlier work.[1]

(iii) Vibrational Power Spectrum

We have computed the vibrational power spectrum from the Fourier transform of the velocity autocorrelation (VAC) function for the relative velocities of a hydrogen atom with respect to the oxygen atom[36]where I, ω, v, and t are the intensity, frequency, velocity, and time, respectively. ⟨(0)·(t)⟩  is the VAC function. The full power spectrum for the bulk water with varying temperature is shown in Figure S1. Although the power spectrum is not expected to mimic the experimental IR/Raman spectra, the qualitative nature of the trend is consistent with prior works.[27,36] The vibrational part for the frequency range 3100–3800 cm–1 of the power spectrum is shown in Figure a. The intensity of the O–H bond stretching band at around ∼3400 cm–1 decreases with an increase in temperature, whereas the intensity of dangling OH bonds or the percentage of non-hydrogen-bonded water increases in the range 3550–3600 cm–1 as expected.

Figure 3

Vibrational spectra of (a) bulk water with temperature variation from 273 to 373 K and (b) first hydration shell water around the hydrophobic molecules for single and multi LJ models with size ranging from 0.4 to 4 nm for the 2D surface of the hydrophobic molecule and bulk water at 300 K.

The full spectrum for the hydration shell water around spherical hydrophobic solutes is shown in Figure S2. A cutoff distance of 0.585 nm has been used to define the hydration shell around the hydrophobic solutes. The frequency range 3100–3800 cm–1 of the power spectrum is shown in Figure b. In this case, we observed that the intensities of the OH stretching band of the hydration shell water fall with an increase in the size of solutes and simultaneously, the intensity of the non-hydrogen-bonded water molecules increases around 3550–3600 cm–1. A noticeable increase in this peak intensity can be seen at and beyond the hydrophobic solute size of 1 nm where the water order–disorder transition takes places as reported previously.[8] A very prominent peak at around ∼3600 cm–1 for the extreme case of lowest curvature of the hydrophobic solute (planar 2D surface) shows that the population of non-hydrogen-bonded water is maximum. The increase in the fraction of dangling H-bonds as a function of temperature and size of the solute has been shown in Figure S5. Thus, all of these results are consistent with the original proposal that changes in the hydration shell structure with increasing solute size closely mimic the increase in temperature for bulk water.

Based on our observations so far, we can safely conclude that there exists a nanometer scale order–disorder transition for spherical hydrophobic solutes, and the water around large solutes exhibits a significant increase in disorder. It is tempting to suggest the analogy that water around small solutes behaves as “cold water”, whereas water around large solutes behaves as “hot water”, but the increase in tetrahedrality is only marginal for small solutes (<1 nm), which may not be enough to explain the observed entropy reduction on dissolving methane in water. The excluded volume effect encountered by the surface water molecules seems to be a more likely explanation for that.[12]

Does the Length-Dependent Order–Disorder Crossover Exist in Linear Alcohols?

In this section, we turn our attention to the system of water-soluble hydrophobic solutes, namely, linear alcohols, that are routinely used in experimental studies as model systems to investigate hydrophobic hydration. We have focused on the spectral and structural properties of the hydration shell water around the hydrophobic part of the alcohols. The list of alcohols studied here is provided in Table S1 along with the effective concentration of a single molecule of alcohol in water. We have also simulated 0.5 M methanol, 1-propanol similar to the experimental studies by Ben-Amotz and co-workers[13,14] as well as 1-heptanol to mimic the experimental conditions closely.

(i) Tetrahedrality and Hydrogen Bonding Distribution

For alcohols, the hydration shell of only carbon atoms has been considered since it is expected that the water molecules around the non-polar part of the solute are expected to demonstrate strongest signatures of hydrophobic hydration, if any. The water molecules near the OH group are expected to have significant polar interactions (including H-bonds). The distribution of Q is demonstrated in Figure a for a single alcohol molecule as well as a 0.5 M concentration. The data for bulk water at 300 K are shown as a reference. For the single molecule as well as for the higher concentrations of alcohols, the distribution of Q appears to be more or less similar to bulk water except for the monomer of methanol. For the monomer of methanol, the population of Q at ∼0.8 is decreased considerably, which is in contrast to the experimental observation that tetrahedral ordering should have a slight increase around methane or methanol.[38] An earlier simulation work by Duboué-Dijon and Laage has reported a small reduction in Q around propanol as well.[33] Though methanol is small in size, the presence of the hydrophilic OH group may significantly alter the hydration shell structure as compared to methane. It has been suggested that the crossover temperature for methane should be significantly lower than that for methanol.[15] Moreover, we do not rule out the possibility of a higher statistical error for methanol since the number of water in the hydration shell is very small, leading to large fluctuations in the data. However, based on several earlier simulation studies, we must note here that the changes in the tetrahedral order parameter due to the presence of a small hydrophobic solute are quite small.[33,40] Therefore, it is difficult to attribute these rather subtle changes to the macroscopic thermodynamic implications of hydrophobic hydration.[41,42]

Figure 4

Distribution of (a) tetrahedral order parameter (Q) and (b) number of H-bonds (NHB) in the hydration shell water around alcohols with single molecules and a 0.5 M concentration. In the calculation of NHB, the OH group of the corresponding alcohol is considered if it is making hydrogen bond with the tagged water molecule.

We have also calculated the number of hydrogen bonds formed by individual hydration shell water molecules using the same conditions mentioned above. It should be noted that in the case of alcohols, we have additionally counted the hydrogen bonds if made between the hydration shell water and the OH group of alcohols. The distribution of number of hydrogen bonds is shown in Figure b. The distribution of monomeric alcohols overlaps with bulk water except for methanol, while a higher concentration of alcohols and methanol slightly decreases in the population of four hydrogen bonds. Overall, the hydrogen bond distribution shows a trend similar to that of water tetrahedrality, and thus, the increasing hydrophobic chain length up to 11 carbon atoms does not have a significant impact on the structural ordering and hydrogen bond distribution in the first hydration shell. This further confirms that the water tetrahedrality and hydrogen bonding do not get significantly affected by the increase in the hydrophobic chain length for linear alcohols at least up to 1-undecanol (n = 11).

(ii) Vibrational Power Spectrum

We have computed the vibrational power spectrum in the hydration shell of linear alcohols using the same protocol described earlier for spherical solutes. The full spectrum is shown in Figure S3, and the O–H stretch band of the spectrum is shown in Figure . The intensities of the spectrum almost overlap with each other irrespective of the increase in the hydrophobic chain length up to 11 carbon atoms for the alcohols. Thus, the experimentally observed increase in the percentage of water dangling OH bonds (in the frequency range of 3550–3660 cm–1) with the chain length of the alcohol could not be reproduced within our simulation framework.

Figure 5

O–H stretch vibrational band of the power spectrum of the first hydration shell water around the alcohols ranging from methanol to undecanol and bulk water at 300 K.

It would be pertinent to clarify here that we are computing the vibrational power spectrum and not experimentally observed Raman (or IR) spectra. Thus, we do not expect a quantitative agreement with the experimentally observed spectra. Moreover, non-polarizable force fields are likely to have quantitative differences with experimental data. The objective of the work reported here is not to reproduce the experimental spectra per se. Rather, we are interested in the qualitative trend of the size dependence of various order parameters and power spectrum for a given force field. Therefore, within the framework of a chosen water model, we argue that although the order–disorder crossover is visible for spherical hydrophobic solutes, it is clearly insignificant for linear alcohols. Hence, the nanometer-scale crossover for spherical solutes cannot be trivially extended to linear molecules. However, it is quite possible that some higher order effects such as polarizability might be responsible for the experimental observations (including Raman spectra). These might be more sensitive probes to the rather subtle changes in the hydration shell structure. Our future work will take up ab initio MD (AIMD) simulations in an attempt to investigate the origin of the trends observed in the experimental spectra.

Conclusions

We have carried out a systematic investigation of the nature of solute size-dependent order–disorder crossover in the hydration layer of hydrophobic molecules. We have used tetrahedral order parameter (Q) and number of hydrogen bonds per water to quantify the structural ordering of the hydration shell. A vibrational power spectrum has been computed using a flexible water model to be able to correlate with the experimental observations of vibrational spectra. All structural and dynamical properties capture the expected solute size-dependent order–disorder crossover around the nanometer length scale. Interestingly, the increase in disorder with increasing solute size is highly analogous to the increase in temperature in bulk water. Hence, in line with the Raman measurements from the Ben-Amotz group, we can establish the analogy that around smaller solutes, the hydration layer behaves as “cold water”, whereas around larger solutes and surfaces, it behaves as “hot water.” Admittedly, the similarity is more pronounced for larger solutes (hot water) rather than for small solutes since the increase in tetrahedrality around methane is rather small as repeatedly observed in several earlier studies as well. Nevertheless, the vibrational power spectrum shows the increasing population of dangling OH bonds with increasing size but only for the spherical solutes.

The solute size dependence for linear alcohols (chain length of the hydrophobic tail) shows a completely different behavior. In an earlier work, we have argued and demonstrated that linear hydrocarbon chains cannot exhibit the length-dependent crossover due to the presence of a sub-nanometer length scale along the cross section of the chain.[11] Thus, water molecules can wrap around the chain to almost maintain the hydrogen-bonded network as if it does not “feel” the longer dimension of the chain. We have shown that the emergence of disorder is not possible until the hydrocarbon starts to fold into a collapsed state with an increase in the length scale. It has also been reported by several groups earlier that n-alkanes cannot exist in the collapsed state for the number of carbons below 20.[43,44] Therefore, it is not surprising that our classical MD simulations do not find any evidence of size-dependent order–disorder transition for alcohols up to 11 carbon atoms.

This presents us with a puzzling question regarding the possible origin of the experimental observation of the increased population of dangling OH bonds around larger alcohols. While our analyses do not even capture the presence of such species under ambient conditions, it is quite possible that non-polarizable flexible water models are not able to capture such details. It is quite possible that the water molecules near an interface encounter a different local electrostatic environment that would alter the polarization effects on the interfacial water as compared to their bulk counterpart. Therefore, it is likely that AIMD simulations or polarizable water models might be able to reproduce these changes in the vibrational spectra.

Nevertheless, we speculate that although the vibrational spectra might be highly sensitive to the local interactions, that is not necessarily an ideal reporter of the formation of a global clathrate-like structure or significant increase in tetrahedrality as observed even in earlier AIMD simulations with hydrophobic solutes.[45] Havenith and co-workers have demonstrated through THz calorimetry that “not the tetrahedrally ordered component but the interstitial hydration water is found to be mainly responsible for the temperature-dependent change in ΔCp and ΔG”.[46] In an interesting computational work, Biswas and co-workers have shown that whether the O–H stretch band in the hydration shell will be blue- or red-shifted as compared to the bulk water would depend on the relative orientation of the water with respect to the solute.[20] The temperature dependence of solvation entropy has been correlated with two different types of hydrogen-bonded water in the hydration shell as well.[48] Thus, different kinds of experimental and computational observables may report different types of local and global structural and dynamical parameters in a highly context-dependent manner.[49] Hence, it would be crucial to put together information from different types of order parameters and thermodynamic variables in order to build a wholesome picture of the nature of hydrophobic hydration. For example, although some changes in the local environment might alter the vibrational dynamics in a subtle way, that does not necessarily correlate with the global thermodynamics of the hydration free energies which are often dominated and controlled by cancellation between enthalpic and entropic factors.[41,42,50] In summary, the hydrophobic effect continues to elude us, and our future work would further focus on building a unifying picture that can consistently explain the seemingly contrasting experimental and computational observations.