Potential of Mean Force of Hydrophobic Association: Dependence on

J. Phys. Chem. B 2007, 111, 10765-10774

10765

Potential of Mean Force of Hydrophobic Association: Dependence on Solute Size Emil Sobolewski,† Mariusz Makowski,†,‡ Cezary Czaplewski,‡ Adam Liwo,‡ Stanisław Ołdziej,‡ and Harold A. Scheraga*,‡ Faculty of Chemistry, UniVersity of Gdan´ sk, ul. Sobieskiego 18, 80-952 Gdan´ sk, Poland, and Baker Laboratory of Chemistry and Chemical Biology, Cornell UniVersity, Ithaca, New York 14853-1301 ReceiVed: January 23, 2007; In Final Form: July 5, 2007

The potentials of mean force (PMFs) were determined for systems involving formation of nonpolar dimers composed of methane, ethane, propane, isobutane, and neopentane, respectively, in water, using the TIP3P water model, and in vacuo. A series of umbrella-sampling molecular dynamics simulations with the AMBER force field was carried out for each pair in either water or in vacuo. The PMFs were calculated by using the weighted histogram analysis method (WHAM). The shape of the PMFs for dimers of all five nonpolar molecules is characteristic of hydrophobic interactions with contact and solvent-separated minima and desolvation maxima. The positions of all these minima and maxima change with the size of the nonpolar molecule, that is, for larger molecules they shift toward larger distances. The PMF of the neopentane dimer is similar to those of other small nonpolar molecules studied in this work, and hence the neopentane dimer is too small to be treated as a nanoscale hydrophobic object. The solvent contribution to the PMF was also computed by subtracting the PMF determined in vacuo from the PMF in explicit solvent. The molecular surface area model correctly describes the solvent contribution to the PMF together with the changes of the height and positions of the desolvation barrier for all dimers investigated. The water molecules in the first solvation sphere of the dimer are more ordered compared to bulk water, with their dipole moments pointing away from the surface of the dimer. The average number of hydrogen bonds per water molecule in this first hydration shell is smaller compared to that in bulk water, which can be explained by coordination of water molecules to the hydrocarbon surface. In the second hydration shell, the average number of hydrogen bonds is greater compared to bulk water, which can be explained by increased ordering of water from the first hydration shell; the net effect is more efficient hydrogen bonding between the water molecules in the first and second hydration shells.

1. Introduction Hydrophobic association plays an important role in many chemical and biological phenomena in aqueous solution. These include formation and stability of various self-assembly aggregates and biological structures in an aqueous environment: proteins, micelles, biological membranes, and macromolecular complexes.1-8 Hydrophobic interactions are also important in surfactant aggregation, coagulation, complexation, detergency, and gas clathrates.9 Hydrophobic association of nonpolar solutes in water arises from solvent-induced interactions, that is, indirect interactions among solutes, caused by specific changes of the structure of water in the vicinity of nonpolar molecules.2 Solvent-induced hydrophobic interactions are described by the free energy of association, or changes of the free energy as a function of the coordinates of the solute molecules, that is, the potential of mean force (PMF) that acts between two or more nonpolar molecules in an aqueous environment. The tendency for hydrophobic association of nonpolar solutes in water can be regarded as a partial reversal of the thermodynamically unfavorable process of hydrophobic hydration.2,8,10 Near room temperature, hydrophobic association is characterized by a positive enthalpy and a positive entropy of association. The favorable entropy of association dominates and provides a net attraction between nonpolar solutes in water. * To whom correspondence should be addressed. Phone: (607) 255 4034; fax: (607) 254 4700; e-mail: [email protected]. † University of Gdan ´ sk. ‡ Cornell University.

Experimental methods provide fundamental thermodynamic data for many hydrophobic phenomena11,12 and help explain the microscopic origins of hydrophobic hydration. X-ray scattering provides a direct measure of the oxygen-oxygen distribution function of water around nonpolar solutes.10,13 Neutron scattering with isotope substitution provides a direct determination of the hydrogen-hydrogen distribution function.14 Radial distribution functions of water around simple nonpolar solutes such as argon15 and methane16 in heavy water have been determined by neutron scattering. Experimentally measured second virial coefficients of mixtures of water with nonpolar substances17,18 also shed light on hydrophobic hydration. All these experiments help to describe the microscopic origins of hydrophobic hydration, but direct experimental information on the microscopic origins of hydrophobic association is not available. Because of the low solubility of nonpolar compounds in water, hydrophobic association of model compounds cannot be measured experimentally, except for a homologous series of molecules with polar head groups and an increasing size of the nonpolar tails;8 therefore, theoretical studies remain the main source of our knowledge about this phenomenon for simple model systems. Most computer simulations of hydrophobic association of simple nonpolar solutes treat dilute aqueous solutions, that is, they focus on the PMF for two nonpolar particles in aqueous solution,19-23 and its dependence on the temperature,2,10,24-27 pressure,28-30 or ionic strength of the solution.31 In our recent investigations, we studied the distance dependence of the

10.1021/jp070594t CCC: $37.00 © 2007 American Chemical Society Published on Web 08/22/2007

10766 J. Phys. Chem. B, Vol. 111, No. 36, 2007

Sobolewski et al.

multibody contributions to the PMF for three and four nonpolar molecules,32-37 and we have shown that the distance dependence of the cooperative multibody contribution is in good agreement with that of the molecular surface area, which confirmed the observation that molecular surface area provides a good description of hydrophobic hydration.2,10,38,39 The effects of solute size and the strength of the interactions between the solute molecules upon hydrophobic association of the nonpolar side chains of amino acids are explained in ref 2. These properties are expressed in eq 11 of ref 2 by the terms ∆YS and ZR. Section 5 of ref 2 also provides a semiquantitative argument for a weaker solvent-separated minimum configuration compared to a stronger contact-minimum configuration. Rank and Baker studied the distance-dependent PMF for a united-atom methane dimer (modeled by a Leonard-Jones potential with parameters σ ) 3.71 Å and ) 0.294 kcal/mol) and a larger united-atom solute dimer (modeled by a LeonardJones potential with parameters σ ) 5.18 Å and ) 0.420 kcal/ mol),39 where the Lennard-Jones (LJ) potential is defined by eq 1

[(σr ) - (σr ) ]

ELJ(r) ) 4

12

6

(1)

with σ being the value of the distance (r) at which the LJ energy is zero and being the energy at the minimum of the LJ potential. Rank and Baker found that the free-energy minimum of the PMF, obtained from simulations in water, is less deep (0.55 kcal/mol) for large solutes than for a methane pair (0.75 kcal/mol) despite the increase in the strength of the solutesolute LJ interactions used in the simulations ( ) 0.420 kcal/ mol compared to ) 0.294 kcal/mol for methane). The PMF is less deep for large solutes because of stronger solute-water interactions. The depth of the PMF depends on direct solute-solute interactions and solvent-induced hydrophobic interactions. For larger particles with consequently larger values of , both the direct solute-solute and solute-water interactions are stronger. This is caused by a larger number of interaction centers for larger particles. Stronger solute-water interactions lead to better solubility. When nonpolar molecules pass from the gas phase to aqueous solutions, the solubility is determined by the balance of two contrasting factors: the excluded volume entropy change due to cavity creation in the solvent and the direct solutesolvent van der Waals interactions.40 The excluded volume entropy change due to cavity creation in the solvent depends only on the solute size and leads to poor solubility of nonpolar compounds in water. Solute size is not the only factor that determines solubility; the strength of the direct solute-solvent van der Waals and induced-dipole interactions is also important. Graziano40 confirmed this assertion by comparison of the solubility of n-alkanes, alk-1-enes, and alk-1-ynes. The size of the three molecules is very similar, but the strength of their van der Waals and induced-dipole interactions with water is very different, leading to different solubility as shown by comparison of experimental values of Gibbs hydration freeenergy changes at 298 K: 1.98, 1.26, and -0.50 kcal/mol for propane, propene, and propyne, respectively.40 Solvent-induced hydrophobic interactions depend on the same factors determining solubility, discussed above, because hydrophobic association can be regarded as a partial reversal of the thermodynamically unfavorable process of hydrophobic hydration.2,8,10 Two other studies, comparing hydrophobic association of methane and isobutylene, were carried out by Tsai and co-

workers,41,42 but they were concerned with the total free energy of association and not with the distance dependence of the PMF. In our recent study,34 we calculated the PMF for the hydrophobic association of pairs and triplets of nonpolar solute molecules with different sizes and strengths of LJ interactions: united-atom methane and two models of large united-atom solutes. As in the work of Rank and Baker,39 we also found that the depth of the PMF minimum for large solutes (σ ) 5.18 Å) depends on the strengths of the LJ interactions: it is less deep for stronger solute-solute LJ interactions ( ) 0.420 kcal/ mol), and it is deeper for weaker solute-solute LJ interactions such as those between united-atom methane ( ) 0.294 kcal/ mol). These results show that, for the united-atom model, the depth of the PMF in water depends on solute-solute and solute-solvent interactions, as discussed in a previous paragraph. Lorentz-Berthelot rules43 were used for water-solute interactions, as in eq 2, where ii and ij refer to homogeneous and heterogeneous pairs, respectively:

ij ) xii‚jj σij )

σii + σjj 2

(2)

In models using mixing rules, an increase in the strength of solute-solute interactions (for fixed strength of the water-water interactions) causes an increase of the strength of solute-water interactions. Such a stronger tendency toward hydration of the hydrophobic particles in water impairs hydrophobic interactions. In effect, the contact minimum in the PMF for a pair of hydrophobic particles is less deep because of stronger solutewater interactions. The effect of solute size on the PMF of hydrophobic association can be investigated with various methods. Recently, we compared two approaches to calculate potentials of mean force: the particle insertion method and the weighted histogram analysis method (WHAM), using neon, united-atom methane, and large united-atom solutes as model systems.37 Huang et al. computed radial and orientational distribution functions of water around three different molecules: argon, methane, and neopentane.44 They also calculated the PMF between two neopentane molecules processing the results of umbrella-sampling simulations by imposing the requirement that the values of the probability distributions are the same at some point in the overlap regions of neighboring windows of the simulation.45 This technique differs from that of WHAM used in our and other work.32,37,46 Southall and Dill47 studied aqueous solvation of a nonpolar solute as a function of its radius using a simplified statistical mechanical [Mercedes Benz (MB)] model of two-dimensional water. They confirmed the classical picture of hydrophobicity: the transfer of a small solute into water involves a large positive heat capacity. Hydrophobic hydration is characterized by a large negative entropy at low temperature and a large positive enthalpy at high temperature. By contrast to the behavior for small nonpolar solutes, the transfer of large nonpolar solutes into water involves no such large changes in heat capacity or entropy.47 In a follow-up paper, Southall and Dill calculated the solute size dependence of the PMF48 and found that the solute size determines the relative preference of the contact and solventseparated minima. Larger solutes prefer the contact state more strongly, while smaller solutes have a greater tendency to become solvent separated,48 as proposed in ref 2. The MB water model has two-dimensional geometrical limitations; for instance,

Potential of Mean Force of Hydrophobic Association

J. Phys. Chem. B, Vol. 111, No. 36, 2007 10767

it did not find a significant maximum in the change of the heat capacity for the desolvation barrier configuration, as postulated by Graziano.49 The nature of hydrophobicity at room temperature is dominated by entropic factors only for a small particle size and becomes dominated by energetic factors for a large particle size.47 Near room temperature, the hydration process for small particles is dominated by the entropy term, which results in a positive free energy of solution and accounts for the low solubility. The solubility of small nonpolar solute species in water decreases when the solution is heated and increases when the system is cooled near room temperature.2,50 This is directly related to the temperature dependence of the solubility of hydrocarbons in water (see Figure 9 of ref 50). For larger particles, that is, mesoscopic or macroscopic size nonpolar particles, for example, planar nonpolar surfaces, long-range hydrophobic interactions can be induced by a dewetting transition or capillary cavitation51 or by the formation of vicinal microbubbles.52,53 Capillary cavitation is closely related to hydrophobic dewetting phenomena.54 More specifically, dewetting is linked to the enhanced probability of observing large cavities in water driven by the strong cohesive interactions between water molecules.54 The crossover of hydrophobic interactions from a small molecular scale to a large macroscopic scale was predicted to occur at about 10 Å length scale in the size of hydrophobic particles.55 In the present paper, we address only small-range hydrophobicity by comparing the PMFs of hydrophobic association of pairs of methane, ethane, propane, isobutane, and neopentane molecules, respectively, in water. The largest molecule studied, neopentane, has a diameter, d, of 5.7 Å which is still below the 10 Å length scale and is qualitatively similar to other smaller solutes.56 2. Methods Sampling of the configurational space, necessary in the umbrella-sampling method, was carried out by using molecular dynamics (MD) simulations for the following all-atom models of hydrophobic dimers, that is, two methane molecules, two ethane molecules, two propane molecules, two isobutane molecules, and two neopentane molecules. The dimers were immersed in TIP3P57 water boxes with periodic boundary conditions with a box side of about 31, 33, 35, 36, and 38 Å for methane, ethane, propane, isobutane, and neopentane, respectively. Molecular dynamics simulations were carried out in two steps. In the first step, each system was equilibrated in the NPT ensemble (constant number of particles, pressure, and temperature) at 298 K for 100 ps. Production MD simulations were then run in the NVT ensemble (constant number of particles, volume, and temperature). In all MD simulations, the integration step was 2 fs. A 9 Å cutoff was used for all nonbonded interactions, and the electrostatic energy was evaluated using the particle-mesh Ewald summation.58 For each system, a series of 20 windows of 10 ns simulations per window were run with different harmonic-restraint potentials imposed on the distance (ξ) between the two atoms (one from each particle) closest to the center of the mass of each of the molecules (indicated by the stars in Figure 1) as given by eq 3.

V(ξ) ) k(ξ - do)2

(3)

with the force constant k ) 2 kcal/mol/Å2 and the “equilibrium” distance (do) for a given window equal to 3.5, 4.0, ..., and 13.0 Å for every nonpolar dimer studied in this work, for windows from 1 to 20. Snapshots from the MD simulations were saved

Figure 1. Partial atomic charges (in electron charge units) of the methane (a), ethane (b), propane (c), isobutane (d), and neopentane (e) molecules calculated by using the RESP method61 on the basis of HF/ 6-31G* calculations carried out with GAMESS60 used in the calculations with the AMBER force field.59 The atoms are labeled with standard AMBER atom-type symbols. The restraints during the MD simulations were imposed at the distances between the atoms labeled with stars.

every 0.2 ps. A total number of 50 000 configurations were collected for each window. The charges on the atoms of the solute molecules needed for the AMBER 7.059 force field were determined for each system using a standard procedure for fitting the point-charge electrostatic potential to the molecular electrostatic potential computed by using the electronic wave function calculated at the restricted Hartree-Fock (RHF) level with the 6-31 G* basis set. The program GAMESS60 was used to carry out the quantummechanical calculations, while the program RESP61 of the AMBER 7.0 package was used to compute the fitted charges. The charges and the AMBER atom types are shown in Figure 1. To determine the potentials of mean force (PMFs) for the systems studied, we processed the results from all windows of the restrained MD simulations for each system using the weighted histogram analysis method (WHAM).46,62 For a given system, one-dimensional histograms, dependent only on the distance between the geometric centers of interacting dimers, were constructed. This means that our PMF curves were averaged over all possible orientations. The bin dimension applied in the WHAM calculations of the PMF was equal to 0.1 Å for all systems. The calculated PMFs should tend to zero with increasing distance (after subtracting the constant factor accounting for the hydrophobic-hydration free energy of an isolated solute molecule). The region 11.5-13.0 Å was used as a baseline to superimpose the PMFs calculated by WHAM. To determine the PMFs for pairs of the nonpolar dimers in vacuo, a similar computational procedure was used as described


Figure 2. PMF curves for the neopentane dimer for two different numbers of configurations obtained from the umbrella-sampling/ WHAM method using the TIP3P water model. The simulations were carried out at a temperature of 298 K. The dotted and solid lines refer to 20% (10 000 configurations) and 100% (50 000 configurations), respectively, of the total number of generated configurations.

above. The only difference was the use of shorter MD simulations, that is, 2 ns per window. Molecular surface area was calculated on the basis of the scheme and analytical equations given by Rank and Baker.39 We used five different radii corresponding to the particles studied in this paper, assuming their shape to be spherical. The values of the radii were 1.950, 2.175, 2.400, 2.625, and 2.850 Å for the methane, ethane, propane, isobutane, and neopentane monomers, respectively. The radius of the water probe was 1.4 Å. Figure 2 shows the dependence of the PMF for the neopentane dimer, obtained using the umbrella-sampling/WHAM method, on the number of configurations collected from each window using the TIP3P water model calculations as an example. As seen from Figure 2, convergence is achieved as the number of configurations increases. The same convergence was observed for the remaining dimers. We also analyzed the packings, interactions, and orientations of water molecules in the vicinity of the solutes. In the umbrellasampling simulations, only the Cartesian coordinates of the solute molecules were stored for the PMF calculations. To determine the distribution, orientation, and number of hydrogen bonds of water molecules around the interacting solutes, we ran an additional 20 ns of MD simulations for the methane and neopentane dimers, respectively, which are formed from the smallest and the largest molecules, respectively, considered in this study; harmonic restraints were imposed to keep the monomers at selected distances (namely, the contact minimum, the desolvation barrier, and the solvent-separated minimum distances for both systems), and the Cartesian coordinates of all atoms were stored every 0.2 ps for further analysis. Because the dimer systems have cylindrical symmetry, we expressed the distributions in the cylindrical coordinates h (the axis passing through the line linking the centers of the solute molecules) and r (an axis perpendicular to h) (see Figure 3); the distributions are averaged over the azimuthal angle θ. We computed the normalized distribution function of the water molecule density, the average number of hydrogen bonds between water molecules, and the average values of the cosine of the angle ω between the dipole moment of a water molecule and the r axis (see Figure 3) at all points of the two-dimensional grid in h

Sobolewski et al.

Figure 3. Schematic illustration of the definitions for the cylindrical distributions: the cylindrical axis h links the solute molecules, r is an axis perpendicular to it, and distributions are averaged over the azimuthal angle θ. The monomer solutes, represented by circles, are shown at the contact distance. A water molecule is presented (at an arbitrary orientation) to define the angle ω between the r axis and the dipole moment vector of the water molecule. The H-O‚‚‚O angle R used to define a hydrogen bond is presented using two water molecules in the hydrogen-bonded dimer geometry.

Figure 4. PMF curves for the dimer systems investigated in this study determined by using the TIP3P model of water: methane (thin solid line), ethane (dashed line), propane (thick solid line), isobutane (dotteddashed line), and neopentane (dotted line).

and r around the dimers. Two-dimensional maps for all distributions were prepared with a 0.2 Å grid. Two water molecules were considered to be hydrogen bonded63,64 if the oxygen‚‚‚oxygen distance was not greater than 3.5 Å and, simultaneously, the H-O‚‚‚O angle R was not greater than 30° (see Figure 3). In addition, we calculated the distributions of the average number of nearest water neighbors for water molecules around the methane and neopentane dimers, using the same 3.5 Å oxygen‚ ‚ ‚oxygen distance cutoff, and the distributions of average binding energies between these nearestneighbor water molecules, normalized per molecule or per pair of interacting molecules. 3. Results and Discussion Potential of Mean Force. Figure 4 shows the results of the calculations of the PMF for methane (thin solid line), ethane (dashed line), propane (thick solid line), isobutane (dotteddashed line), and neopentane (dotted line) dimers, using the umbrella-sampling/WHAM method with the TIP3P water model. The PMF curves have characteristic shapes with one deep minimum each at about 3.9, 4.5, 4.9, 5.3, and 5.8 Å, corresponding to methane, ethane, propane, isobutane, and neopen-

Potential of Mean Force of Hydrophobic Association tane, respectively. This minimum is referred to as the contact minimum. The deepest contact minima are observed for the isobutane (dotted-dashed line) and the neopentane (dotted line) dimers with depths of about -0.9 and -1.05 kcal/mol, respectively. Less deep contact minima are observed for the remaining nonpolar dimers and are about -0.75 kcal/mol. The depths of the contact minima for the ethane (dashed lines) and propane dimers (thick solid line) are slightly smaller than for the methane dimers. This feature of all-atom simulations differs from that discussed in the Introduction for united-atom simulations and might be caused by the fact that ethane and propane are prolate and not spherical molecules, their PMFs having been averaged over all possible orientations. The PMF curves shown in Figure 4 also contain the second minima at distances of about 7, 7.6, 8.2, 8.8, and 9.2 Å for methane, ethane, propane, isobutane, and neopentane, respectively. These so-called solvent-separated minima correspond to distances at which precisely one water molecule can enter the space between the two solutes.2,65 The first maximum in the PMF curve between the contact and solvent-separated minima is referred to as the desolvation maximum. The lowest desolvation maximum (approximately 0.2 kcal/mol) at the shortest distance is observed for the methane dimer (thin solid line), and the highest is observed at the longest distance for the neopentane dimer (dotted line). As seen from Figure 4, the positions and heights of these maxima change according to the size of the solute. The same variations of positions are observed for the solvent-separated minimum of dimers investigated in this study. All PMFs of the nonpolar dimers presented in this study have similar depths of the solventseparated minimum, being about -0.15 kcal/mol. For the neopentane dimer, Huang et al.44 observed a much deeper contact minimum (about -2.75 kcal/mol). Huang et al. postulated that the hydration of the neopentane dimer was similar to the hydration of a macroscopic nonpolar object and was different from that of argon and methane dimers. However, we do not observe any large differences between the PMF curves of the neopentane dimers and the dimers of the other smaller molecules. With the simulation method that they used (described in the Introduction), Huang et al. showed results for only 400 ps per window. Their PMF curve for neopentane is quite rugged, and the convergence of their PMF calculations was not discussed. Our results are in agreement with experimental results reported by Graziano,56 who presented some thermodynamic data for the hydration process of argon, methane, and neopentane. He did not observe any qualitative differences in thermodynamic features between methane and neopentane; he observed only quantitative ones. Figure 5 shows the PMFs obtained in vacuo for the nonpolar particles investigated in this study. All PMFs obtained in vacuo have the characteristic shape of the Lennard-Jones-like potential. As with the curves obtained in explicit solvent, the positions and depths of the minima change with the size of the solutes. The minimum with the smallest depth is observed for the methane dimers (thin solid line), being about -0.3 kcal/mol at 4.1 Å. For the ethane (dashed line) and propane (thick solid line) dimers, the depths of the minima are comparable, being about -0.45 kcal/mol and -0.43 kcal/mol at positions 5.0 and 5.4 Å, respectively. The deepest minima in the PMFs obtained in vacuo are for the isobutane (dotted-dashed line) and neopentane (dotted line) dimers, being about -0.88 kcal/mol and -0.95 kcal/mol at positions 5.6 and 6.0 Å, respectively. The positions of these minima in the PMF curves for all the dimers are shifted in vacuo toward greater distances with


Figure 5. The PMF curves for the dimer systems investigated in this study determined in vacuo: methane (thin solid line), ethane (dashed line), propane (thick solid line), isobutane (dotted-dashed line), and neopentane (dotted line).

Figure 6. The solvent contribution to the PMF determined as a difference between the PMFs in solvent and in vacuo: methane (thin solid line), ethane (dashed line), propane (thick solid line), isobutane (dotted-dashed line), and neopentane (dotted line).

increasing size and differ in magnitude from the positions of the contact minima obtained in water. For all dimers, the contact minimum in water occurs at shorter distances than in vacuo. The reason for this behavior is the stronger hydrophobic interactions in water than LJ interactions in vacuo between these particles. The largest shift of the position of the contact minimum (along the abscissa) toward smaller distances in water (approximately 0.5 Å) is observed for the ethane and propane dimers, being about 0.35 Å for the isobutane and neopentane dimers, and the smallest shift is observed for the methane dimer, being about 0.2 Å. For all dimers except isobutane, the contact minimum in water is deeper than in vacuo, demonstrating that the influence of water is more important than that of van der Waals forces in hydrophobic association.2 In Figure 6, we present the solvent contribution to the PMF determined as a difference between the PMFs in solvent and those in vacuo. As seen from Figure 6, the solvent contribution to the PMF has one maximum corresponding to the position of the desolvation maximum in the PMF curves for the nonpolar dimers in water (see Figure 4). The lowest maxima are observed for the methane (thin solid line), ethane (dashed lines), and propane (thick solid line) dimers, being about 0.3, 0.4, and 0.38


Figure 7. Molecular surface area calculated with the analytical equations derived by Rank and Baker39 with different radii of the monomeric solute particles equal to 1.950 (thin solid line), 2.175 (dashed line), 2.400 (thick solid line), 2.625 (dotted-dashed line), and 2.850 Å (dotted line).

kcal/mol, respectively. We observe the highest maxima for the solvent contribution to the PMFs for the isobutane (dotteddashed line) and neopentane dimers (dotted line), being about 0.68 kcal/mol. On the basis of the scheme and analytical equation given by Rank and Baker,39 we computed the molecular surface area for pairs of nonpolar solutes of radii 1.950 Å (thin solid line), 2.175 Å (dashed line), 2.400 Å (thick solid line), 2.625 Å (dotteddashed line), and 2.850 Å (dotted line). The plots are shown in Figure 7. It can be seen that the positions and heights of the maxima increase with the radius of the solute and are the lowest for the smaller values of the radius and are the highest for the solute with the largest value of the radius. A similar situation is observed for the solvent contribution to the PMFs shown in Figure 6, which is attractive at the contact minimum in Figures 6 and 7 and hence goes to large negative values at the contact minimum and at smaller distances. However, the molecular surface area was calculated for ideal spherical systems, while only two (methane and neopentane) out of five solutes used in our study are close to spherical in shape. The solvent contributions were averaged over all possible orientations for all five dimers (see Methods). The dependence of the molecular surface area on the distance within the dimer provides a good approximation to the solvent contribution to the PMF. However, it cannot capture such details as the solvent-separated minima or the variations of the desolvation barrier for nonspherical ethane and propane. The distance at which the desolvation barrier occurs is described correctly by the molecular surface area. The molecular surface area shows the correct tendency of a higher barrier for larger molecules. Distribution of Water Molecules. An analysis of the structure of water in the vicinity of the solute particles can provide new physical insight about the molecular origin of hydrophobic interactions. We calculated two-dimensional distribution functions describing the normalized water density, the average number of hydrogen bonds between water molecules, and the average value of the cosine of the angle ω of the orientation of water molecules with respect to the r axis (see Figure 3). Two-dimensional maps of the normalized distribution function of the water molecule density around the neopentane and methane dimers at three selected distances between the centers

Sobolewski et al. of the solute molecules (3.90, 5.95, and 7.20 Å for the methane dimer and 5.80, 7.85, and 9.20 Å for the neopentane dimer) which correspond to the contact minimum, the desolvation barrier, and the solvent-separated minimum configurations, respectively, are shown in Figure 8. The normalized distribution function of the water density increases in the first and second solvation spheres around the solute molecules compared to the bulk water density. It can be seen that the increase is not uniformly distributed (yellow and red) within the first solvation sphere. Water accumulates (red) close to the bisector plane (h ) 0) between the neighboring neopentane molecules, where the first solvation spheres of two molecules overlap as in Paschek’s study of the Xe-Xe pair.30 The normalized distribution function of the water density around the methane dimer, shown in Figure 8 at corresponding distances, is very similar; the only difference is that the methane molecule and the solvation sphere around it are approximately 25% smaller than those for neopentane. Figure 9 shows the cylindrical distributions of the average number of hydrogen bonds between water molecules around the neopentane and methane dimers. In both systems, the water molecules in the first solvation shell (green) have more broken hydrogen bonds than in bulk water. The average number of hydrogen bonds between water molecules close to the solutes is in good agreement with previous studies of water distributions around single methane and neopentane molecules and is smaller for neopentane (blue) than for methane as observed by Huang et al.44 (however, contrary to the results of that study,44 we observed no qualitative difference between the PMF of neopentane and that of smaller molecules). A smaller number of hydrogen bonds in the first solvation sphere can easily be explained because these water molecules are in contact with the solute and, hence, are no longer completely surrounded by other water molecules, that is, water in the first solvation shell cannot form hydrogen bonds with the solutes. A slightly larger average number of hydrogen bonds in the second solvation sphere (orange) arises because of a more-ordered structure in the first solvation sphere compared to that in bulk water. The average number of hydrogen bonds is not distributed uniformly within the first solvation shell of the dimer and is always smallest in the region between the solute molecules where two solvation spheres overlap. The average number of hydrogen bonds depends both on the average number of nearest water neighbors and on the ordering of the water orientations necessary to form hydrogen bonds. Comparison of the distributions of the average number of hydrogen bonds with the average number of nearest water neighbors around methane and neopentane dimers, calculated using the same oxygen‚‚‚oxygen distance cutoff (3.5 Å) to define a hydrogen bond, shows that ordering of the orientations has only a secondary effect. In all layers, the distribution of the average number of hydrogen bonds always has lower values than the average number of nearest water neighbors, but the shapes of the distributions are almost identical (see Figures 9d and 10a). The average number of nearest-neighbor water neighbors is lower in the first hydration shell than in bulk water. However, the orientation ordering is important for the pairwise water-water binding energies, that is, the strength of the hydrogen bonds. Figure 10b shows the distribution of the average water-water interaction energy between nearest-neighbor water molecules around the neopentane dimer, and Figure 10c shows the average water-water interaction energy normalized per water-water contact. The pairwise water-water binding energies are more negative in the first solvation shell of hydrophobic solute, that is, the hydrogen bonds are stronger in the first hydration shell.



Figure 8. Normalized distribution functions of the water molecule density in the vicinity of the methane dimer at monomer separation distances ∆h of (a) 3.90, (b) 5.95, and (c) 7.20 Å and the neopentane dimer at monomer separation distances ∆h of (d) 5.80, (e) 7.85, and (f) 9.20 Å which correspond to the contact minimum, the desolvation barrier, and the solvent-separated minimum configurations, respectively. The color scale is shown above the panels, and the bulk water density is displayed in white. The solute is in gray. The figure reveals that there is both a clear first hydration shell (maximum in distribution function next to solute in green and yellow with even higher values in red between two solutes) and a weaker second hydration shell (maximum in cyan with values only slightly higher than for bulk water) separated by regions of lower density than bulk water (in blue and violet).

Paschek30 also observed both a reduction of the average number of water neighbors and enhanced population of strongly bound water molecules in the first hydration sphere of a dimer of xenon atoms. A smaller average number of hydrogen bonds between water molecules close to the solutes is not in contradiction with the traditional explanation of the hydrophobic effect (i.e., the reorganization of the structure of water around nonpolar solutes which creates ordered, low-entropy, flickering clathrate-like structure), as in the model of Nemethy and Scheraga,66 because these water molecules close to the solutes have more ordered orientations; even if they participate in a smaller number of hydrogen bonds, these hydrogen bonds are stronger. This observation is clearly illustrated in Figure 11 displaying the average value of the cosine of the angle ω between the dipolemoment vector of a water molecule and the normal to the axis of the methane or neopentane dimers. It can be seen that the water molecules close to the solutes are not oriented at random as in bulk water (which would correspond to cos ω ) 0), but the maximum value of cos ω is ≈0.05 for the methane dimer (yellow) and is ≈0.1 for the neopentane dimer (red); these values of cos ω imply a small net component of the dipole moment pointing outward from the molecular surface. This ordering of the water molecules in the first hydration sphere is reflected in the increase of the strength of the hydrogen bonds (Figure 10)

and the increase of the number of hydrogen bonds of the water molecules from the second hydration sphere with respect to bulk water (Figure 9). The model of Nemethy and Scheraga66 predicted an increase in the average number of water-water hydrogen bonds in the first hydration shell of a hydrocarbon molecule. However, this model was based on a discrete definition of the first hydration shell and of a hydrogen bond without considering the variations of the strengths of the hydrogen bonds. Also, this model did not consider situations in which the hydrogen bonds cannot be defined because of distorted geometry in which the water molecules nevertheless still participate in a low-energy structure. The definition of energy levels in terms of hydrogen bonding in the model of Nemethy and Scheraga could result in overemphasizing the ordering of water molecules, which is suggested by the fact that the changes of heat capacity upon addition of hydrocarbon to water predicted by this model were greater than the experimental values (Table 7 in ref 66). On the other hand, the TIP3P model used in our study underestimates the ordering of water molecules, which could result in underestimating the number of hydrogen bonds. 4. Conclusions We determined the potentials of mean force for five pairs of nonpolar solutes, that is, methane, ethane, propane, isobutane,


Sobolewski et al.

Figure 9. Distribution of the average number of hydrogen bonds between water molecules in the vicinity of the methane dimer at monomer separation distances ∆h of (a) 3.90, (b) 5.95, and (c) 7.20 Å and the neopentane dimer at monomer separation distances ∆h of (d) 5.80, (e) 7.85, and (f) 9.20 Å which correspond to the contact minimum, the desolvation barrier, and the solvent-separated minimum configurations, respectively. The color scale is shown above the panels, and the average number of hydrogen bonds for bulk water density is displayed in white.

Figure 10. Distribution of (a) the average number of nearest water neighbors, calculated using the same oxygen‚‚‚oxygen distance cutoff (3.5 Å) to define a hydrogen bond as for the average number of calculated hydrogen bonds presented in Figure 9; (b) the average water-water interaction energy between nearest-neighbor water molecules; (c) the average water-water pair energy calculated by dividing the average water-water interaction energy by the average number of nearest water neighbors, i.e., normalized per water-water contact. All distributions are calculated around the neopentane dimer at 5.80 Å distance. The color scale is shown above the panels, and the average values for bulk water are displayed in white.

and neopentane, in water as functions of the distance between their geometric centers. The shape of the PMFs determined by using the TIP3P model of water is characteristic of hydrophobic interactions with contact and solvent-separated minima and desolvation maxima. We demonstrated that the neopentane dimer is too small to be treated as a macroscopic hydrophobic

object. The results of the PMF calculations are sensitive to the parameters of umbrella sampling (number of windows and number of snapshots collected in each window) and to the method for processing the umbrella-sampling results to remove a biased potential.32-37,46 The shape of the PMFs determined in vacuo is characteristic of the Lennard-Jones-like interaction



Figure 11. Distribution of the cosine of the angle ω between the vector normal to the axis of the dimer and the water dipole moment vector in the vicinity of the methane dimer (see Figure 3) at monomer separation distances ∆h of (a) 3.90, (b) 5.95, and (c) 7.20 Å and the neopentane dimer at monomer separation distances ∆h of (d) 5.80, (e) 7.85, and (f) 9.20 Å. The color scale is shown above the panels, and the value of the average cosine of the angle ω of bulk water is displayed in white.

for nonpolar molecules. The changes in the height and maximum positions of the solvent contribution to the PMF are qualitatively well described by the molecular surface area. An analysis of the structure of water in the vicinity of the solute particles shows a weakly ordered structure in the first hydration shell. The average number of hydrogen bonds for the water molecules close to the solutes is smaller, but the bonds are stronger than for bulk water; also, the orientations of these water molecules are not as random as those in bulk water. The first solvation shell is not static as in crystalline clathrates, and experimental results support some similarity between water in a clathrate and water around nonpolar solutes: 13C chemical shifts of aqueous methane upon magic-angle spinning NMR spectra are virtually identical to the chemical shifts of methane in the dodecahedral cage of a crystalline methane hydrate, and the hydration number of aqueous methane is 20.67 The hydration shell of aqueous methane can be described as one in which rather long-lived hydration shells with 20 water molecules dominate the methane hydration process but in which the methane and water molecule configurations rapidly change so that the dodecahedral cages of a rigid methane hydrate are not present.67 On the other hand, a direct comparison of the first solvation shell of krypton atoms with the actual crystalline cage observed by extended X-ray absorption fine structure (EXAFS) shows considerable differences between this crystalline cage structure and the first hydration shell in the liquid.68 Water ordering effects were found, but they were attributed to the second solvation shell of the krypton atom. The coordination number

of Kr in liquid water solution is significantly lower than in crystalline clathrates,69 which is not in agreement with our results. We are not in position to reconcile the differences between the NMR and EXAFS results except to point out the large size of methane (used in the NMR experiment) compared to the smaller size of Kr (used in the EXAFS experiment). In our opinion, the traditional explanation of the hydrophobic effect pointing to the reorganization of the structure of water around nonpolar solutes which creates an ordered, low-entropy first solvation shell66 is not in contradiction to the view of the importance of cavity creation and the small size of water molecules.70 These two views emphasize different aspects of the same phenomenon and, consequently, are complementary to each other. Nevertheless, to which extent the different reorganization of the structure of water around nonpolar solutes and the excluded volume terms contribute to hydrophobic interaction is actually the subject of an ongoing debate.30,71-73 Ben-Naim74 argues that the reorganization of the structure of water around nonpolar solutes does not affect the Gibbs free energy of hydrophobic interaction but affects only the entropy and enthalpy changes, which results in heat capacity changes. Similarly, Lee75 argues that only the excluded volume term plays a fundamental role in the Gibbs free energy changes induced by hydrophobic hydration. Acknowledgment. This work was supported by grants from the U.S. National Science Foundation (MCB05-41633), the U.S.

10774 J. Phys. Chem. B, Vol. 111, No. 36, 2007 National Institutes of Health (GM-14312), the National Institutes of Health Fogarty International Center Grant TW007193, and the Polish Ministry of Science and Education (1 T09A 099 30 and BW/8000-5-0396-7). Mariusz Makowski was supported by a grant from the “Homing” program of the Foundation for Polish Science (FNP). This research was conducted by using the resources of (a) our 818-processor Beowulf cluster at the Baker Laboratory of Chemistry and Chemical Biology, Cornell University, (b) the National Science Foundation Terascale Computing System at the Pittsburgh Supercomputer Center, (c) our 45-processor Beowulf cluster at the Faculty of Chemistry, University of Gdansk, (d) the Informatics Center of the Metropolitan Academic Network (IC MAN) in Gdansk, and (e) the Interdisciplinary Center of Mathematical and Computer Modeling (ICM) at the University of Warsaw. References and Notes (1) Kauzmann, W. AdV. Protein Chem. 1959, 14, 1. (2) Ne´methy, G.; Scheraga, H. A. J. Phys. Chem. 1962, 66, 1773. Erratum: J. Phys. Chem. 1963, 67, 2888. (3) Scheraga, H. A.; Ne´methy, G.; Steinberg, I. Z. J. Biol. Chem. 1962, 237, 2506. (4) Poland, D. C.; Scheraga, H. A. J. Phys. Chem. 1965, 69, 2431. Erratum: J. Phys. Chem. 1965, 65, 4425. (5) Tanford, C. The hydrophobic effect: formation of micelles and biological membranes; Wiley: New York, 1973. (6) Ben-Naim, A. Hydrophobic interactions; Plenum Press: New York, 1980. (7) Ravishanker, G.; Beveridge, D. L. Theoretical studies of the hydrophobic effect. In Theoretical chemistry of biological systems; Na´raySzabo´, G., Ed.; Elsevier: Amsterdam, 1986; Chapter 7, pp 423-494. (8) Scheraga, H. A. J. Biomol. Struct. Dyn. 1998, 16, 447. (9) Blokzijl, W.; Engberts, J. B. F. N. Angew. Chem., Int. Ed. Engl. 1993, 32, 1545. (10) Griffith, J. H.; Scheraga, H. A. J. Mol. Struct. (THEOCHEM) 2004, 682, 97. (11) Blandamer, M. J.; Fox, M. F. Spectroscopic properties. In Water: A ComprehensiVe Treatise; Franks, F., Ed.; Plenum: New York, 1973; Vol. 2, Chapter 8, p 459. (12) Enderby, J. E.; Neilson, G. W. X-Ray and Neutron Scattering by Aqueous Solutions of Electrolytes. In Water: A ComprehensiVe Treatise; Franks, F., Ed.; Plenum: New York, 1979; Vol. 6, Chapter 1, p 1. (13) Hura, G.; Russo, D.; Glaeser, R. M.; Head-Gordon, T.; Krack, M.; Parrinello, M. Phys. Chem. Chem. Phys. 2003, 5, 1981. (14) Turner, J.; Soper, A. K. J. Chem. Phys. 1994, 101, 6116. (15) Broadbent, R. D.; Neilson, G. W. J. Chem. Phys. 1994, 100, 7543. (16) De Jong, P. H. K.; Wilson, J. E.; Neilson, G. W.; Buckingham, A. D. Mol. Phys. 1997, 91, 99. (17) Rigby, M.; Prausnitz, J. M. J. Phys. Chem. 1968, 72, 330. (18) Abdulgatov, I. M.; Bazaev, A. R.; Ramazanova, A. E. J. Chem. Thermodyn. 1993, 25, 249. (19) Owicki, J. C.; Scheraga, H. A. J. Am. Chem. Soc. 1977, 99, 7413. (20) Rapaport, D. C.; Scheraga, H. A. J. Phys. Chem. 1982, 86, 873. (21) Smith, D. E.; Haymet, A. D. J. J. Chem. Phys. 1993, 98, 6445. (22) van Belle, D.; Wodak, S. J. J. Am. Chem. Soc. 1993, 115, 647. (23) Young, W. S.; Brooks, C. L., III. J. Chem. Phys. 1997, 106, 9265. (24) Lu¨demann, S.; Schreiber, H.; Abseher, R.; Steinhauser, O. J. Chem. Phys. 1996, 104, 286. (25) Chau, P. L.; Mancera, R. L. Mol. Phys. 1999, 96, 109. (26) Shimizu, S.; Chan, H. S. J. Chem. Phys. 2000, 113, 4683. (27) Ghosh, T.; Garcia, A. E.; Garde, S. J. Am. Chem. Soc. 2001, 123, 10997. (28) Ghosh, T.; Garcia, A. E.; Garde, S. J. Chem. Phys. 2002, 116, 2480. (29) Paschek, D. J. Chem. Phys. 2004, 120, 6674. (30) Paschek, D. J. Chem. Phys. 2004, 120, 10605. (31) Mancera, R. L. Chem. Phys. Lett. 1998, 296, 459. (32) Czaplewski, C.; Rodziewicz-Motowidło, S.; Liwo, A.; Ripoll, D. R.; Wawak, R. J.; Scheraga, H. A. Protein Sci. 2000, 9, 1235.

Sobolewski et al. (33) Czaplewski, C.; Rodziewicz-Motowidlo, S.; Liwo, A.; Ripoll, D. R.; Wawak, R. J.; Scheraga, H. A. J. Chem. Phys. 2002, 116, 2665. (34) Czaplewski, C.; Ripoll, D. R.; Rodziewicz-Motowidło, S.; Liwo, A.; Wawak, R. J.; Scheraga, H. A. Int. J. Quantum Chem. 2002, 88, 41. (35) Czaplewski, C.; Rodziewicz-Motowidlo, S.; Dabal, M.; Liwo, A.; Ripoll, D. R.; Scheraga, H. A. Biophys. Chem. 2003, 105, 339. Erratum: Biophys. Chem. 2004, 111, 267. (36) Czaplewski, C.; Liwo, A.; Ripoll, D. R.; Scheraga, H. A. J. Phys. Chem. B 2005, 109, 8108. (37) Czaplewski, C.; Kalinowski, S.; Liwo, A.; Scheraga, H. A. Mol. Phys. 2005, 103, 3153. (38) Tun˜o´n, I.; Silla, E.; Pascual-Ahuir, J. L. Protein Eng. 1992, 5, 715. (39) Rank, J. A.; Baker, D. Protein Sci. 1997, 6, 347. (40) Graziano, G. J. Chem. Soc., Faraday Trans. 1998, 94, 3345. (41) Tsai, J.; Gerstein, M.; Levitt, M. Protein Sci. 1997, 6, 2606. (42) Raschke, T. M.; Tsai, J.; Levitt, M. Proc. Natl. Acad. Sci. U.S.A. 2001, 98, 5965. (43) Allen, M. P.; Tildesley, D. G. Computer simulation of liquids; Oxford University Press: New York, 1987; Chapter 1.4.2, p 21. (44) Huang, X.; Margulis, C. J.; Berne, B. J. J. Phys. Chem. B 2003, 107, 11742. (45) Pangali, C.; Rao, M.; Berne, B. J. J. Chem. Phys. 1979, 71, 2975. (46) Kumar, S.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A.; Rosenberg, J. M. J. Comput. Chem. 1992, 13, 1011. (47) Southall, N. T.; Dill, K. A. J. Phys. Chem. B 2000, 104, 1326. (48) Southall, N. T.; Dill, K. A. Biophys. Chem. 2002, 101, 295. (49) Graziano, G. J. Chem. Phys. 2005, 123, 034509. (50) Ne´methy, G.; Scheraga, H. A. J. Chem. Phys. 1962, 66, 3401. (51) Carambassis, A.; Jonker, L. C.; Attard, P.; Rutland, M. W. Phys. ReV. Lett. 1998, 80, 5357. (52) Hummer, G.; Garde, S. Phys. ReV. Lett. 1998, 80, 4193. (53) Christenson, H. K.; Claesson, P. M. Science 1988, 239, 390. (54) Ben-Amotz, D. J. Chem. Phys. 2005, 123, 184504. (55) Lum, K.; Chandler, D.; Weeks, J. D. J. Phys. Chem. B 1999, 103, 4570. (56) Graziano, G. J. Phys. Chem. B 2004, 108, 9371. (57) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (58) Darden, T.; York, D.; Pedersen, L. J. Chem. Phys. 1993, 98, 10089. (59) Case, D. A.; Pearlman, D. A.; Caldwell, J. W.; Chaetham, T. E., III; Wang, J.; Ross, W. S.; Simmerling, C. L.; Darden, T. A.; Merz, K. M.; Stanton, R. V.; Cheng, A. L.; Vincent, J. J.; Crowley, M.; Tsui, V.; Gohlke, H.; Radmer, R. J.; Duan, Y.; Pitera, J.; Massova, I.; Seibel, G. L.; Singh, U. C.; Weiner, P. K.; Kollman, P. A. AMBER 7; University of California: San Francisco, CA, 2002. (60) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. A.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.; Windus, T. L.; Dupuis, M.; Montgomery, J. A. J. Comput. Chem. 1993, 13, 1347. (61) Bayly, C. I.; Cieplak, P.; Cornell, W. D.; Kollman, P. A. J. Phys. Chem. 1993, 97, 10269. (62) Kumar, S.; Rosenberg, J. M.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A. J. Comput. Chem. 1995, 16, 1339. (63) Luzar, A.; Chandler, D. Nature 1996, 379, 55. (64) Starr, F.; Nielsen, J.; Stanley, H. Phys. ReV. E 2001, 62, 579. (65) Geiger, A.; Rahman, A.; Stillinger, F. H. J. Chem. Phys. 1979, 70, 263. (66) Nemethy, G.; Scheraga, H. A. J. Chem. Phys. 1962, 38, 3401. (67) Dec, S. F.; Bowler, K. E.; Stadterman, L. L.; Koh, C. A.; Sloan, E. D., Jr. J. Am. Chem. Soc. 2006, 128, 414. (68) Finney, J. L.; Bowron, D. T.; Daniel, R. M.; Timmins, P. A.; Roberts, M. A. Biophys. Chem. 2003, 105, 391. (69) Ashbaugh, H. S.; Asthagiri, D.; Pratt, L. R.; Rempe, S. B. Biophys. Chem. 2003, 105, 323. (70) Southall, N. T.; Dill, K. A.; Haymet, D. J. J. Phys. Chem. B 2002, 106, 521. (71) Ashbaugh, H. S.; Truskett, T. M.; Debenedetti, P. G. J. Chem. Phys. 2002, 116, 2907. (72) Graziano, G. J. Chem. Phys. 2003, 119, 10448. (73) Ashbaugh, H. S.; Truskett, T. M.; Debenedetti, P. G. J. Chem. Phys. 2003, 119, 10450. (74) Ben-Naim, A. Biopolymers 1975, 14, 1337. (75) Lee, B. Biopolymers 1991, 31, 993.

Potential of Mean Force of Hydrophobic Association: Dependence on

Recommend Documents