J. Phys. Chem. B 1997, 101, 4343-4348
4343
Hydrophobic Effect, Water Structure, and Heat Capacity Changes Kim A. Sharp* and Bhupinder Madan Department of Biochemistry & Biophysics, UniVersity of PennsylVania, 3700 Hamilton Walk, Philadelphia, PennsylVania 19104-6059 ReceiVed: January 16, 1997; In Final Form: April 1, 1997X
The hydration heat capacity (∆Chyd p ) of nine solutes of varying hydrophobicity was studied using a combination of a random network model equation of water and Monte Carlo simulations. Nonpolar solutes cause a concerted decrease in the average length and angle of the water-water hydrogen bonds in the first hydration shell, while polar and ionic solutes have the opposite effect. This is due to changes in the amounts, relative to bulk water, of two populations of hydrogen bonds: one with shorter and more linear bonds and the other with longer and more bent bonds. Heat capacity changes were calculated from these changes in water structure using a random network model equation of state. The calculated changes account for observed differences in ∆Chyd p for the various solutes. The simulations provide a unified picture of hydrophobic and polar hydration, and both a structural and thermodynamic explanation of the opposite signs of ∆Chyd p observed for polar and nonpolar solutes.
Introduction In recent years analysis of heat capacity changes has become of central importance in understanding the thermodynamics of protein folding, protein-protein binding, protein-nucleic acid binding, and the hydrophobic effect.1-7 There are several reasons for this. First, the improved sensitivity and accuracy of calorimeters has provided a wealth of heat capacity (Cp) data for these processes. Second, Cp relates the other three major thermodynamic quantitiessenthalpy (H), entropy (S), and free energy (G)sto each other
∂S ∂H ∂2G 〈∆H 〉 )T ) -T2 2 ) ∂T ∂T ∂T kT2 (i) (ii) (iii) (iv) 2
Cp )
(1)
where T is the temperature, k is Boltzmann’s constant, and 〈∆H2〉 is the mean-squared fluctuation in enthalpy. These four equivalent definitions of Cp illustrate the diverse thermodynamic implications for processes that involve changes in Cp. For example, unfolding of many proteins is accompanied by a large increase in Cp,1,2 implying (eqs 1i and 1ii) that the unfolding entropy and enthalpy depend strongly on temperature and implying (eq 1iii) that the stability curve (∆Gunfold vs T) will be an inverted U-shape with the consequent possibility of both heat and cold denaturation. Changes in Cp associated with protein folding and binding come almost entirely from changes in hydration heat capacity (∆Chyd p ) due to burial (desolvation) of polar and nonpolar groups.3,8,9 Burial of nonpolar groups also provides the dominant stabilizing interaction (the hydrophobic interaction) in protein folding and binding.10,11 The hydrophobic interaction is driven by the low solubility of nonpolar compounds in water due, at room temperature, to the unfavorable decrease in the entropy (increase in structure) of the hydrating waters. This has been described in many ways: as “icelike structures”, icebergs, flickering clusters, clathrates, and pentagonal clusters.12-14 It has become clear, however, that a description in terms of entropic changes alone does not explain the distinctive features of the hydrophobic effect.15 First, the large positive ∆Chyd for solvation of p * To whom correspondence should be addressed. X Abstract published in AdVance ACS Abstracts, May 1, 1997.
S1089-5647(97)00245-9 CCC: $14.00
nonpolar groups means that at higher temperatures (T > 80 °C) their insolubility results from unfavorable changes in enthalpy, not entropy. Second, the hydration of a majority of polar and ionic solutes is also accompanied by a decrease in entropy, and thus by some kind of structuring of water, but in this case is negative.7,16-20 Thus, analysis of Cp changes has ∆Chyd p become essential for a full understanding of polar and nonpolar hydration and the hydrophobic effect. In spite of the importance of ∆Chyd p changes in a wide range of processes involving solutes, proteins, and nucleic acids, basic features remain unexplained at the molecular level. From the data of Privalov and Gill2 the increase in heat capacity for nonpolar groups amounts to about 3 cal/(mol K) for each water in the first hydration shell or 20% of the total heat capacity of a single water molecule. There is no clear picture of how an increase in water structure can give rise both to this remarkably large positive ∆Chyd for nonpolar solutes and yet result in a p negative ∆Chyd for polar and ionic groups. Explanations for p the positive ∆Chyd of hydrophobic hydration in terms of the p induction of more ordered water which is “melted” at increasing temperatures15,21,22 do not explain the negative ∆Chyd p for polar and ionic groups. The latter phenomenon is perhaps the more puzzling because the entropy of polar and ionic hydration is usually negative. Thus, one deduces from eq 1ii that the ordering of water induced by these solutes increases with increasing temperature. This runs counter to the general expectation that it becomes harder to induce order (decrease entropy) at higher temperatures. Monte Carlo or molecular dynamics simulations of solutes in water using standard water and solute potential functions have provided much structural and thermodynamic information about hydration. However, calculation of ∆Chyd directly from such p simulations via eq 1 using standard water and solute potential function appears to be impractical without prohibitive size and length computer simulations.23 Even determination of ∆Chyd p by the most robust approach, using eq 1i, requires taking the difference of a difference in mean enthalpy from at least four separate simulationssbetween pure water and water plus one solutesand evaluating this difference at two or more temperatures. Direct determination of ∆Chyd p from simulations is thus subject to considerably more statistical uncertainty than the determination of the hydration free energy or enthalpy. © 1997 American Chemical Society
4344 J. Phys. Chem. B, Vol. 101, No. 21, 1997
Sharp and Madan
Figure 1. (a) Definition of H-bond length, d, and angle, θ. (b) Hydrogen bond structure and heat capacity changes. The solid surface shows the dependence of the RN model incremental heat capacity ∆Crn p (vertical axis) on the mean H-bond length, d, and the rms H-bond angle, θ (horizontal axes). Values of the bond length and angle (and the associated heat capacity change) for the first hydration shell of solutes obtained from Monte Carlo simulations are shown for the nonpolar (2), polar (9), and mixed (×) H-bond classes. The solutes are water (W), argon (A), ethanol (E), methane (M), K+ (K), Cs+ (C), N-methylacetamide (N), tetramethylammonium (T), I- (I), and benzene (B). (c) Linear regression analysis of H-bond length vs angle.
We show here how an equation of state for ∆Chyd derived p from a random hydrogen-bond (H-bond) network model of water can be combined with simulations of water around solutes to provide statistically reliable estimates of ∆Chyd p . Application of this approach to nine solutes of varying hydrophobicity provides a unified description of hydrophobic, polar and ionic hydration that explains at a molecular level how ordering of water around solutes can result in either an increase or a decrease in Cp.
with contributions from H-bond bending, stretching, vibrations, a zero point distortion term, and van der Waals interactions, strecth zpd vdw Cbend , Cvib p , Cp p , Cp , and Cp , respectively. These five contributions, expressed as cal/K per mole of H-bonds, are given by 2
Cbend ) 0.139e-θ /2 + 1.51e-5.41θ p
6
bend stretch (d,s,θ) + Cvib Crn p (d,s,θ) ) Cp (d,s,θ) + Cp p (d,s,θ) + vdw Czpd p (d,s,θ) + Cp (d,s,θ) (2)
(3)
Cstretch ) 2.10(1 - 3.72x) - 45.88σ2 p
Theory and Methods The random network (RN) model is based on the observation that liquid water retains much of the tetrahedral H-bonding network of ice I, albeit in a randomly distorted form due to thermal fluctuations, and that the essential features of this network can be described in terms of just three structural parameters: the mean and standard deviation in H-bond length, d, and s, respectively, and the root-mean-square (rms) H-bond angle, θ.24,25 The H-bond length and angle used in the RN model are defined in Figure 1a. The RN model provides an accurate description of the structure of pure liquid water in terms of the oxygen-oxygen radial distribution function, and it reproduces subtle properties of water, such as the bulk dielectric constant over the range 283-573 K, with remarkable accuracy.26-28 The RN model also provides good equations of state for the free energy, enthalpy and entropy of pure liquid water over a wide range of temperatures,26-28 from which one can derive an equation of state for the heat capacity in terms of the H-bond structural parameters23
2
Cvib p ) 0.9936T
∑ i)1
[
k Fi T
() Fi
coth
Fi -
T
T sinh
]
( )
()
2
k Fi
Fi T T
Fi
-
T2
2
-1.240θ -2.79x Czpd p ) 0.235e
) Cvdw p
(
)
(4)
(5)
(6)
0.242e-θ 1 1 + 12 θ (1 + x) (1 + x)6 12 6 0.558(1 - e-θ) (7) (1 + x)7 (1 + x)13
(
)
where
Fi )
hνi -1.24θ2-2.79x -0.381hνi -1.24θ2-2.79x e e , k Fi ) (8) 2k 2000k
where the mean and standard deviation in H-bond length are expressed in dimensionless form as x ) (d - 2.75)/2.75 and σ2 ) [(s/2.75)2 + x2], respectively, h is Planck’s constant, and νi ) 555, 660, 840, 60, 230, and 315 cm-1 for i ) 1-6, respectively, are the reference frequencies of the three transla-
Hydration Heat Capacity Simulations
J. Phys. Chem. B, Vol. 101, No. 21, 1997 4345
tional and three librational intermolecular vibrational modes of water. These equations were derived as described previously,23 but their presentation has been simplified somewhat by combining the different numerical constants and expressing the contributions per mole of H-bonds rather than per mole of water molecules. They are based on the Henn and Kauzmann25 form of the RN model which incorporates the constant pressure condition under which experiments and simulations are performed. In the RN model for pure water d, s, and θ are obtained selfconsistently at any given temperature by minimizing the total free energy with respect to these H-bond structural parameters.25 The presence of a solute, however, will perturb the H-bonding structure of the surrounding water. If do, so, and θo are the average H-bond structural parameters for water in the absence of solute, i.e., for bulk water, then the change in heat capacity from an H-bond with solute-perturbed values of d, s, and θ is given by rn rn ∆Crn p (d,s,θ) ) Cp (d,s,θ) - Cp (d0,s0,θ0)
(9)
The function ∆Crn p (d,s,θ) defined by eq 9 describes the incremental heat capacity per H-bond as a function of solute induced H-bond distortions. The dependence of the incremental heat capacity function on the two most sensitive RN variablessthe mean length, d, and the rms angle, θsis shown in Figure 1b. An increase in H-bond length and rms angle decreases the heat capacity contribution, while decreases in length and angle have the opposite effect. The heat capacity also depends to a lesser extent on s, decreasing as this variable increases. Using eq 9, the net Cp change produced by a solute can be obtained by summing the contributions from all of the hydrogen bonds it perturbs:
∆Chyd p )
∑i Ni∆Crnp (di,si,θi)
(10)
where ∆Crn p (di,si,θi) is the contribution to heat capacity change arising from a group of Ni perturbed H-bonds with average parameters di, si, and θi. The summation in eq 10 over different sets of H-bonds accounts for the fact that different H-bonds may be perturbed to different degrees, depending on their distance from the solute, what hydration shell they are in, and the polarity of the different solute groups neighboring the water. In addition to the H-bond contribution to ∆Chyd p , potential contributions can arise from (i) the fluctuation in the pressurevolume contribution to the hydration enthalpy, (ii) fluctuations in the intrasolute energy, and (iii) fluctuations in the mean solute-solvent interaction energy, Esw.23 At 1 atm pressure, the absolute value of the PV term even for the largest solute studied here is ,0.01 kcal/mol, so fluctuations in this term are negligible. Even for large flexible solutes such as proteins contribution ii is small,9 and so for the smaller, more rigid solutes here this term can also be ignored. Contribution iii can in principle be evaluated directly from simulations by obtaining the variation in mean solute-solvent interaction energy with temperature, but for most of the solutes shown here, this derivative cannot be distinguished statistically from zero in the simulations. In the one case where there is a statistically significant contribution (for argon23), it is a factor of 2 less than the water contribution. In the present study we restricted ourselves to an analysis of the direct water H-bond contribution to hydration heat capacity only. To evaluate d0, s0 and θ0 for pure water and di, si, and θi for water surrounding a solute, Monte Carlo simulations of pure water and a solute in water were carried out as described
previously,23 and eq 10 used to evaluate Cp changes caused by the solute-induced water distortions. Nine solutes encompassing a diverse group of properties were examined, including four ions (Cs+, K+, I-, and tetramethylammonium (TMA+)), two purely nonpolar solutes (argon and methane), an aromatic solute (benzene), and two solutes containing polar and nonpolar groups (ethanol and N-methylacetamide). Simulations were carried out with 216 water molecules, using the TIP4P potential29 which reproduces the structure and key thermodynamic properties of bulk water well, and is widely used in explicit water simulations. The solute and its interaction with water were modeled using the OPLS potential function which, in conjunction with the TIP4P water potential, provides accurate hydration free energies for a wide range of solutes.30 Sampling of water configurations was carried out using a Metropolis Monte Carlo scheme implemented in the program BOSS31 under constant temperature and pressure conditions (T ) 298 K, P ) 1 atm). The systems were equilibrated for 50 × 106 Monte Carlo steps, and then configurations sampled every 1000 steps from 10 consecutive runs of 10 × 106 steps to compute the desired quantities. Solute flexibility was included for multiatom solutes. The solute-water radial distribution function was obtained from the simulations and used to determine which hydration shell each water belonged to. Thus, waters lying within the first peak of this distribution function were classified as first shell, and so on. Each water was also classified as polar or nonpolar hydrating based on whether the nearest solute atom or group was polar or nonpolar respectively. All pairs of hydrogen-bonded waters were identified, using a distance cutoff (taken as 3.4 Å, the boundary of the first peak in the O-O radial distribution function of pure water), the hydrogen bond type was classified as polar, nonpolar or mixed according to the type of both waters and the hydrogen bond length and angle were evaluated. Mean values of d, s and θ, and their probability distributions were accumulated over the course of the simulation. Standard deviations were obtained from the variation in batch means between different Monte Carlo runs. Convergence was checked by monitoring the batch averages of energy and hydrogen-bond structural parameters. The mean structural parameters, d, s, and θ were found to converge satisfactorily with these size and length simulations, unlike direct estimates obtained via differences in enthalpy and from of ∆Chyd p fluctuations in enthalpy.23 The contribution of water to ∆Chyd for the ionic solutes is p potentially long-ranged, and to include the effect of water beyond the first hydration shell, a Born correction was used. Application of eq 1iii to the Born expression for hydration free energy yields
Cpborn )
(
[ ])
q2T 1 ∂2 2 ∂ 2acut 2 ∂T2 3 ∂T
2
(11)
where θ is the charge, is the dielectric constant of water, and acut is the cutoff radius for explicit determination of solvent contributions, taken at the outer boundary of the first hydration shell, at 3.5 Å for K+, 4 Å for Cs+, 6.4 Å for TMA+, and 3.5 Å for I-. The experimental temperature dependence of the water bulk dielectric constant was taken from ref 32. Experimental values of ∆Chyd for CsI, KI, and TMAI were p taken from refs 33 and 34, for argon from ref 2, and for ethanol, methane, N-methylacetamide, and benzene from ref 16. Results and Discussion Data on the structure of hydrogen bonds between two first hydration shell waters, between two second shell waters,
4346 J. Phys. Chem. B, Vol. 101, No. 21, 1997 between a first and a second shell water, etc., were collected separately. Statistically significant changes in d and θ were seen only for first shell-first shell H-bonds. The standard deviation in hydrogen-bond length, s, also showed no significant change for any group of H-bonds and thus has a negligible role in heat capacity changes. Significant changes in structure were only seen for mean values of d and θ for water-water H-bonds in the first hydration shell. The results for each solute are shown in Figure 1b. Three classes of H-bonds are distinguished in this figure: (i) “Nonpolar” if they are between waters that each hydrate a nonpolar solute or group; (ii) “polar” if they are between waters that each hydrate a polar (or ionic) solute or group; (iii) “mixed” if one water is hydrating a nonpolar group and the other a polar group (present in ethanol and Nmethylacetamide only). Waters hydrating ions or polar groups show an increase in mean H-bond length and rms angle compared to the case of pure water, while waters around nonpolar solutes or groups show decreased mean H-bond lengths and rms angles. H-bonds of the “mixed” type show smaller distortions because of the opposing influence of the polar and nonpolar groups, but in the same direction as the nonpolar solutes. Interestingly, distortions in H-bond length and angle take place in a highly concerted fashion, falling on the same straight line regardless of solute type (Figure 1c, d ) 2.74 Å + 0.067θ, correlation coefficient (r2) ) 0.94). This suggests that one can in fact collapse the length and angle distortions into a single structural distortion coordinate, with the hydrophobic groups at one end of the spectrum, and small ions at the other end. The H-bond length of 2.74 ( 0.02 Å corresponding to an angle of θ ) 0° (representing a hypothetical nonpolar solute of very large structuring power) is very close to the value of 2.75 Å seen in low-temperature (T < -175 K) ice I where also θ ≈ 0°,35 reinforcing the idea that temperature and solutes act in a similar way as perturbants of water structure. In fact, the concept of a “structural temperature” to describe hydration waters is often used in experimental studies.36 Figure 1b also shows that each set of d, θ values obtained from the simulations locates that group of H-bonds at particular point on the incremental heat capacity surface, ∆Crn p (d,s,θ), with the pure water point (d ) 2.95 Å, θ ) 29.4°) lying on the ∆Cp ) 0 contour. From these values and eq 10 the net ∆Chyd for each solute was calculated. For the ionic groups, p the values for the neutral iodide salt were obtained by summing the cation and anion contributions. This allows for comparison with experimental measurements which can directly access of neutral salts. Comparison of experimental and ∆Chyd p calculated hydration heat capacity changes (Figure 2) shows a very good correlation (linear regression ∆Chyd p (calc) ) 0.47∆ Chyd + 0.33 cal/(mol K), r2 ) 0.98). Figures 1b and 2 show p that the combination of random network model and water simulations captures the most striking feature of hydration heat capacities: the opposite behavior of polar and nonpolar groups. as one goes from The switch from negative to positive ∆Chyd p salts containing small inorganic ions such as Cs+ and K+ to larger hydrophobic ions such as TMA+ 17 is also seen in the simulations. Eisenberg and Kauzmann have identified two major contributions to the heat capacity of water: structural and vibrational.35 The relative importance of these two contributions is difficult to obtain from experiments, but estimates can be obtained from these simulations using eqs 2-8. For nonpolar solutes, large positive contributions to the hydration heat capacity come from the structural (bending and stretching) terms, with a smaller, negative contribution from the vibrational term (Table 1). For
Sharp and Madan
Figure 2. Comparison of experimental and calculated ∆Chyd p obtained using the data of Figure 1b combined with eq 10.
TABLE 1: Structural and Vibrational Contributions to the Hydration Heat Capacity % of total solute
bending
stretching
vibrational
K Cs+ TMA+ Ar CH4 benzene
-64 -53 41 69 60 54
-59 -78 82 47 54 63
23 31 -23 -16 -14 -17
+
polar solutes, the relative magnitudes of structural and vibrational contributions are similar but are of opposite sign. This should be contrasted with the relative contributions to the total heat capacity of pure water in the RN model, where only 33% comes from the structural part and 66% from the vibrational part. The zero-point and van der Waals terms make a negligible contribution to ∆Chyd p for all the solutes examined here. The results in Figure 1 are consistent with the standard view that water in some way becomes more “icelike” around nonpolar groups and less “ice-like” around polar groups. A more detailed analysis of the structural changes reveals a more complex picture. The H-bond angle probability distribution for the ionic solutes and polar groups shows that a large population have angles greater than 30o, with a prominent peak around 50°60°, while a smaller population has a more linear geometry with a peak centered around 12° (Figure 3A). This distribution is a consequence of the strong electrostatic solute-solvent interaction which tends to align the water dipoles along the solutewater axis (inset, Figure 3A), so that many, but not all, H-bonds within the first shell are substantially bent compared to bulk water. For nonpolar solutes, the first hydration shell geometry is best understood in comparison to pure water. The H-bond angle probability distribution for pure water shows a large peak again centered around 12° (Figure 3B). There is, however, a small but significant second population centered at about 50°. Pure water is known to have a coordination number of 4.4-4.8.35 Examination of water configurations from the computer simulations shows that the low angle peak arises from the majority of water molecules that are typically coordinated in a quasi-
Hydration Heat Capacity Simulations
J. Phys. Chem. B, Vol. 101, No. 21, 1997 4347
Figure 4. Heat capacity as a function of the difference in energy levels, ∆E, in the two-state model, eq 12. The energy difference and relative heat capacity corresponding to water in bulk solvent and around hydrophobic and polar solutes are indicated schematically in the plot.
population, that decreases θ and thus increases Cp. In this analysis, one might say that nonpolar groups do not so much increase the ordering of water as decrease the disordering of water. The division between the two H-bond classes occurs around θ ) 35° (Figure 3A,B), which we note is similar to the angle used to distinguish broken and made H-bonds in recent neutron scattering studies of water38 and in a computer simulation of water H-bond dynamics.39 The overall picture that emerges from our structural analysis is one of two populations of H-bond angles, with nonpolar groups decreasing the population with more distorted geometry and polar groups increasing their population. The same trend is also seen in the H-bond length probability distribution (not shown), but it is less pronounced. In fact, a two-state model, which is the simplest thermodynamic model with a nonzero heat capacity, may be usefully employed to interpret the simulation results. For a system with only two energy levels, separated by a gap ∆E, application of eq 1iv gives Figure 3. H-bond angle probability distribution functions for pure water (solid line) and the first hydration shell of (A) Cs+ (b) and the hydroxyl group of ethanol (0). (B) Methane (2) and the methyl group of ethanol (+).
tetrahedral fashion around another water molecule (inset, Figure 3B), while the second peak arises from an occasional extra or mismatch water37 that, due to thermal fluctuations, approaches within H-bonding distance. This water is constrained by the other waters, however, to approach from the direction of the “face of the tetrahedron” where it must make a bent H-bond with the central water at an average angle, it turns out close to half the tetrahedral angle (≈109°/2). A nonpolar solute, which can interact only weakly and nondirectionally with water through van der Waals interactions, must nevertheless coexist with the considerable amount of structure that persists in liquid water. It thus tends to occupy the same position with respect to the quasi-tetrahedral lattice that the more weakly interacting water would, in effect competing for this mismatch site (inset, Figure 3B). It is the reduction of the high angle population, not a decrease in the average angle of the more linear H-bond
Cp )
[
]
∆E2 ∆E2 1 e-∆E/kT [p p ] ) 2 2 1 2 -∆E/kT kT kT 1 + e 2 + e-∆E/kT
(12)
where p1 and p2 are the probabilities of the system being in the lower and higher energy states, respectively. For there to be a significant heat capacity there must be a sizable fluctuation in energy, so the two states must be of different energy and yet significantly populated. Thus if the energy gap is very small, the heat capacity will be small, as the ∆E2/kT term is small. Conversely, if the energy gap is very large, then the probability of populating the upper energy level, p2, is very small, and again the fluctuation, and thus the heat capacity, is small. Increasing ∆E at first will increase Cp as the first factor in eq 12 increases, but further increases in ∆E will decrease Cp because p2 decreases. The two-state heat capacity function plotted in Figure 4 illustrates this behavior. The results of our random network/ explicit water simulation may thus be interpreted in terms of a two-state model as follows. Nonpolar solutes weakly perturb the surrounding water, creating a greater disparity in energy of different water microstates, but not so great that the waters
4348 J. Phys. Chem. B, Vol. 101, No. 21, 1997 cannot fluctuate between these states. The predominant effect is to increase the energy difference term, resulting in a positive heat capacity of hydration. In constrast, polar and small ionic solutes strongly perturb the surrounding water. Although this creates a still larger disparity in energy different water states, the hydrating waters are so strongly associated with the solute that the waters cannot fluctuate so easily between these states. This can result in a net decrease in heat capacity due to a large decrease in the probability term. This description of the effect of nonpolar and polar solutes is illustrated schematically in the plot of the two-state heat capacity function, Figure 4. Of course, the random network/explicit water model describes a more complex and realistic potential surface with many more than two states, but the qualitative effects of solute pertubation are made clear from the two-state analogy. The idea of rapid exchange between somewhat energetically different water states, perhaps represented by the two populations of linear and bent hydrogen bonds seen in Figure 3, relates back to the idea of flickering clusters40,41 and is consistent with a recent computer simulation of hydrogen-bond kinetics in pure water39 where rate constants for making and breaking hydrogen bonds between neighboring waters by rotational motions were determined to be kf ) 1 ps-1 and kb ) 0.7 ps-1, respectively, which corresponds to a small effective energy of formation, given by ∆E ) -kT ln (kf/kb), of -0.21 kcal/mol. In summary, the use of the random network model now makes it feasible to obtain statistically reliable estimates of hydration heat capacities for both polar and nonpolar solutes from explicit water simulations and to interpret the results directly in terms of changes in water structure. It should be emphasized that although the explicit water simulations described here require considerably amounts of computer time, they represent the only method of which we are aware where statistically reliable values of ∆Chyd p can be obtained without prohibitive length simulation and where, as importantly, the results can be interpreted directly in structural terms. Nevertheless, our calculated values of show a systematic underestimation of the experimental ∆Chyd p numbers by about 50%, and investigation of this discrepancy is one future direction for further study. One reason for the underestimate may be that there are more subtle perturbations in water structure beyond the first hydration shell which, because the number of hydrogen bonds involved is large, could provide significant additional heat capacity contributions. This is particularly true for charged solutes, and the longer range contributions were estimated here using a continuum solvent model. A continuum model for heat capacities is limited in accuracy, however. Significantly larger and longer simulations would be needed to reliably detect small changes in water structure beyond the first shell. An alternative avenue for the improvement of heat capacity calculations may be by reparametrization. Because the random network and TIP4P water models were developed using different sets of experimental data, they may need to be reparametrized for optimal accuracy. This is an inherent possibility in any hybrid model, although no parametrization was attempted here based on the limited number of solutes studied. Independent of future developments in these type, of heat capacity calculations, the current model already provides an
Sharp and Madan explanation for several important features of small solute hydration heat capacity behavior, it provides a unified picture of polar and nonpolar solvation, and it provides new insights into the phenomenon of hydrophobicity. This approach should open up opportunities for further study of hydration heat capacities of small solutes and larger molecules such as proteins. Acknowledgment. We thank Bill Jorgensen, Erin Duffy, and Dan Severance for advice on using BOSS, Paul Axelson and William De Grado for a critical reading of the manuscript, and Hugh Gallagher for help with Figure 1. Financial support is acknowledged from NIH (GM54105) and the E. R. Johnson Research Foundation. References and Notes (1) Baldwin, R. Proc. Natl. Acad. Sci. U.S.A. 1986, 83, 8069-8072. (2) Privalov, P. L.; Gill, S. J. AdV. Protein Chem. 1988, 39, 191-234. (3) Makhatadze, G.; Privalov, P. J. Mol. Biol. 1990, 213, 385-391. (4) Spolar, R.; Livingstone, J.; Record, M. T. Biochemistry 1992, 31, 3947-55. (5) Spolar, R. S.; Record, M. T. Science 1994, 263, 777-84. (6) Murphy, K. P.; Privalov, P. L.; Gill, S. J. Science 1990, 247, 559561. (7) Murphy, K.; Freire, E. AdV. Protein Chem. 1992, 43, 313-361. (8) Sturtevant, J. Proc. Natl. Acad. Sci. U.S.A. 1977, 74, 2236-40. (9) Velicelebi, G.; Sturtevant, J. M. Biochemistry 1979, 18, 1180-6. (10) Kauzmann, W. AdV. Protein Chem. 1959, 14, 1-63. (11) Dill, K. A. Biochemistry 1990, 29, 7133. (12) Frank, H. S.; Evans, M. W. J. Chem. Phys. 1945, 13, 507. (13) Nemethy, G.; Scheraga, H. J. Chem. Phys. 1962, 36, 3382-3400. (14) Tanford, C. H. The Hydrophobic Effect; John Wiley and Sons: New York, 1980. (15) Kauzmann, W. Nature 1987, 325, 763-4. (16) Cabani, S.; Gianni, P.; Mollica, V.; Lepori, L. J. Solution Chem. 1981, 10, 563-595. (17) Marcus, Y. Biophys. Chem. 1994, 51, 111-128. (18) Murphy, K.; Gill, S. Thermochim. Acta 1990, 172, 11-20. (19) Privalov, P.; Makhatadze, G. J. Mol. Biol. 1993, 232, 660-679. (20) Makhatadze, G.; Privalov, P. J. Mol. Biol. 1990, 213, 375-384. (21) Gill, S.; Dec, S.; Olofsson, G.; Wadso, I. J. Phys. Chem. 1985, 89, 3758-61. (22) Nemethy, G.; Scheraga, H. J. Chem. Phys. 1962, 36. (23) Madan, B.; Sharp, K. A. J. Phys. Chem. 1996, 100, 7713-21. (24) Rice, S. A.; Sceats, M. G. J. Phys. Chem. 1981, 85, 1108-1119. (25) Henn, A. R.; Kauzmann, W. J. Phys. Chem. 1989, 93, 3770-3783. (26) Sceats, M. G.; Rice, S. A. J. Chem. Phys. 1979, 72, 3248-62. (27) Sceats, M. G.; Stavola, M.; Rice, S. A. J. Chem. Phys. 1979, 70, 3927-38. (28) Sceats, M. G.; Rice, S. A. J. Chem. Phys. 1980, 72, 6183-91. (29) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J Chem. Phys. 1983, 79, 926. (30) Jorgensen, W. L.; Tirado-Rives, J. J. Am. Chem. Soc. 1988, 110, 1657-1666. (31) Jorgensen, W. L. BOSS, Version 3.3, Yale UniVersity, New HaVen, CT, 1992. (32) Lide, D. R. CRC Handbook of Chemistry and Physics; CRC Press: Boca Raton, FL, 1990. (33) Lewis, G. N.; Randall, M.; Pitzer, K. S.; Brewer, L. In Thermodynamics; McGraw-Hill: New York, 1961; pp 224-226, 244-249. (34) Friedman, H. L.; Krishnan, C. V. In Thermodynamics of Ion Hydration; Franks, F., Ed.; Plenum Press: New York, 1973. (35) Eisenberg, D.; Kauzmann, W. The Structure and Properties of Water; Oxford University Press: Oxford, 1969. (36) Verrall, R. E. In Infrared Spectroscopy of Aqueous Electrolyte Solutions; Franks, F., Ed.; Plenum Press: New York, 1973. (37) Matubayasi, N. J. Am. Chem. Soc. 1994, 116, 1450-6. (38) Teixeira, J.; Bellisent-Funel, M. C.; Chen, S. H. J. Phys.: Condens. Matter 1990, 2 (Suppl. A), 105. (39) Luzar, A.; Chandler, D. Nature 1996, 379, 55-57. (40) Frank, H. S.; Wen, W. Y. Discuss. Faraday Soc. 1957, 24, 133. (41) Nemethy, G.; Scheraga, H. J. Chem. Phys. 1964, 41, 680.
