Molecular Solvation in Water−Methanol and Water−Sorbitol Mixtures

J. Phys. Chem. B 2007, 111, 4467-4476

4467

Molecular Solvation in Water-Methanol and Water-Sorbitol Mixtures: The Roles of Preferential Hydration, Hydrophobicity, and the Equation of State Prateek P. Shah and Christopher J. Roberts* Department of Chemical Engineering, UniVersity of Delaware, Newark, Delaware 19716 ReceiVed: December 22, 2006; In Final Form: February 27, 2007

Molecular dynamics simulations of aqueous mixtures of methanol and sorbitol were performed over a wide range of binary composition, density (pressure), and temperature to study the equation of state and solvation of small apolar solutes. Experimentally, methanol is a canonical solubilizing agent for apolar solutes and a protein denaturant in mixed-aqueous solvents; sorbitol represents a canonical “salting-out” or protein-stabilizing cosolvent. The results reported here show increasing sorbitol concentration under isothermal, isobaric conditions results in monotonic increases in apolar solute excess chemical potential (µex 2 ) over the range of experimentally relevant temperatures. For methanol at elevated temperatures, increasing cosolvent composition ex results in monotonically decreasing µex 2 . However, at lower temperatures µ2 exhibits a maximum versus cosolvent concentration, as seen experimentally for Ar in ethanol-water solutions. Both density anomalies and hydrophobic effectsscharacterized by temperatures of density maxima and apolar solute solubility minima, respectivelysare suppressed upon addition of either sorbitol or methanol at all temperatures and compositions simulated here. Thus, the contrasting effects of sorbitol and methanol on solute chemical potential cannot be explained by qualitative differences in their ability to enhance or suppress hydrophobic effects. Rather, we find µex 2 values across a broad range of temperatures and cosolvent composition can be quantitatively explained in terms of isobaric changes in solvent densitysi.e., the equation of statesalong with the corresponding packing fraction of the solvent. Analysis in terms of truncated preferential interaction parameters highlights that care must be taken in interpreting cosolvent effects on solvation in terms of local preferential hydration.

I. Introduction Solvation of apolar solutes in water and mixed-aqueous solvents is of widespread importance in a variety of industries, including pharmaceuticals,1-5 food science,6-10 and chemical separations.11-15 Models of apolar solvation in mixed-aqueous solvents have been proposed at varying levels of detail, depending to some extent on the size and complexity of the solutes of interest relative to the solvent components.16-28 Despite significant effort and progress, quantitative and in some cases even qualitative prediction of solvation thermodynamics in mixed-aqueous solvents based on molecular interactions and/ or structural attributes of the cosolvent species remains an outstanding challenge.3,4,16-19,21,22,25-31 Mixed-aqueous solvents can be generally considered as liquid mixtures composed of water and one or more highly miscible cosolvents (or solutes) such as alcohols, polyols, carbohydrates, and organic or inorganic salts.3,32 We adopt here the commonly used Scatchard notation33 for mixed-aqueous solvents in which water is component 1, the cosolvent or cosolute is component 3, and the added solute of interest is component 2. Often, component 2 is either a low molecular weight compound or a relatively high molecular weight biological molecule such as a protein. In both cases, component 2 (the solute) often has relatively limited solubility on a molar or mole fraction basis, and therefore its infinite-dilution excess chemical potential (µex 2 ) determines its distribution and reactivity within a given cosolvent mixture.18,34,35 The work reported here focuses on * Corresponding author. E-mail: [email protected].

spherical, apolar solutes of atomic and molecular dimensions (0.3-0.5 nm diameters) as the simplest realistic models of low molecular weight species such as methane and apolar amino acid side chains such as alanine.36 The two cosolvents selected here are methanol and sorbitol. Experimentally, methanol-water mixtures are typically found to increase the solubility of apolar compounds (relative to pure water),20 as well as to denature proteins.37,38 The latter effect is sometimes argued to be due to increased affinity of methanol for hydrophobic side chains from the interior of folded proteins.38 The former effect has the added complexity that it may be qualitatively different at low versus high alcohol content; e.g., Ar solubility at 280 K and atmospheric pressure, relative to that in pure water, decreases upon addition of ethanol until ca. 15 mol %, reaches a minimum, and then increases upon further addition of ethanol.39 Conversely, sorbitol-water mixtures, like many sugar-water solutions, are typically found to decrease the solubility of macromolecular solutes40 and to promote folding of proteins.41,42 From a chemical perspective, sorbitol and methanol are quite similar. Sorbitol is essentially what one would obtain through the thought experiment of linearly bonding the carbons of six methanol molecules while maintaining their sp3 hybridization, with a particular choice of stereochemistry in the resulting polyol. Bearing in mind the experimental behavior described above, it is rather striking that despite relatively similar chemical natures these water-cosolvent mixtures display such qualitatively different behavior in terms of solvation of apolar solutes at similar w/w% cosolvent. In addition, the relatively simple,

10.1021/jp0688714 CCC: $37.00 © 2007 American Chemical Society Published on Web 04/06/2007

4468 J. Phys. Chem. B, Vol. 111, No. 17, 2007

Shah et al.

nonionic structures of sorbitol and methanol make them amenable to detailed molecular simulation studies over a broad range of temperatures, densities, and compositions. Together, these factors help to motivate our selection of these cosolvents for study here. A characteristic thermodynamic property describing mixedaqueous solvents is the isothermal, isobaric transfer free energy (∆µtr2 ) for a small apolar solute (component 2) between a solvent composed of pure water and a mixed aqueous solvent of a given molality (m3) of component 3. Formally, ∆µtr2 is given by eq 1 for a particular choice of temperature (T), pressure (p), and composition (m3) ex ∆µtr2 (T,p,m3) ) µex 2 (T,p,m3) - µ2 (T,p,m3 ) 0)

(1)

If addition of a given cosolvent results in a decrease in ∆µtr2 , it is considered a solubilizing or a salting-in agent. Conversely, if ∆µtr2 is increased, the cosolvent is typically referred to as a salting-out agent.43,44 For proteins, it is often argued that salting-in or chaotropic cosolvents decrease ∆µtr2 for both folded and unfolded states but more so for the unfolded state.3 As a result they (net) act as denaturants, shifting the folding-unfolding equilibrium toward unfolded conformations. A similar argument leads to the qualitative heuristic that salting-out or kosmotropic cosolvents are expected to stabilize folded conformations by increasing ∆µtr2 for unfolded conformations to a larger extent than ∆µtr2 increases for folded conformations.3,45 In considering the thermodynamics of solvation and transfer between water and mixed-aqueous solvents, a common approach is to categorize cosolvents in terms of the preferential interaction parameter (Γ32) between components 2 and 3. Purely thermodynamically, Γ32 can be defined on a molality scale as46,47

Γ32 )

( ) ∂m3 ∂m2

)-

µ3,T,p

( ) ∂µ2 ∂µ3

)-

m2,T,p

and is related to ∆µtr2 via

∆µtr2 (T,p,m3) ) - kBT

∫0m

3

( ) ( )

( )

Γ32

∂lna3 ∂m3

∂µ2 ∂m3

m2,T,p

∂µ3 ∂m3

m2,T,p

T,p,m2)0

dm3

(2)

(3)

In eq 3, a3 is the thermodynamic activity of component 3 in the binary 1-3 solution (m2 ) 0 for infinite dilution of 2) at a given (T, p, m3), kB is Boltzmann’s constant, and Γ32 is retained inside the integral because it can in principle be a function of m3 at the specified T and p. On the basis of thermodynamic stability criteria for a single-phase mixture,48 (∂lna3/∂m3)T,p,m2)0 must be positive. Therefore, the sign of Γ32 necessarily determines the sign of the derivative of ∆µtr2 with respect to m3. Provided Γ32 does not change sign throughout the range of m3 in the integral in eq 3, the sign of ∆µtr2 will also be dictated by the sign of Γ32. A rigorous, formally exact statistical mechanical interpretation of Γ32 is offered by Kirkwood-Buff (KB) theory.49 Specifically, Γ32 with component 2 at infinite dilution is given by eq 4 at a given T, volume (V), and for the equilibrium density, composition, and pressure corresponding to specific µ1 and µ3 values in

an open system.25,30

Γ32 ) F3(G23 - G21) )

(

〈N2N3〉

〈N2〉〈N3〉

-

〈N2N1〉 〈N2〉〈N1〉

)

(4)

T,V,µ1,µ3

In the above equation, G23 and G21 are defined as

G2j )

∫0∞ [gj2j(r) - 1]4πr2 dr

j ) 1,3

(5)

with gj2j(r) defined as the center-center radial distribution function (rdf) between components 2 and j in the open system (T, V, µ1, µ3, m2 f 0). Ni denotes the number of molecules of component i, and the brackets indicate an equilibrium average in the grand canonical ensemble. Equation 4 and analogous expressions49,50 show how composition fluctuations are correlated in mixtures. Γ32 quantifies the relative magnitude of 2-3 correlations in composition compared to 2-1 correlations. If Γ32 is positive, addition of component 2 (solute) to the control volume will tend to promote addition of component 3 (cosolvent) over the addition of component 1 (water) and vice versa. Physically, this is often interpreted as indicating that a positive (negative) Γ32 corresponds to preferential accumulation (exclusion) of cosolvent in (from) the immediate vicinity of the solute. Rigorously, however, preferential interaction parameters represent a sum of interactions over many solvation shells, not simply those in the immediate proximity of component 2.21,30,51 This result will prove particularly relevant in the context of methanol-water solutions in Section III. An alternative view of cosolvents that is historically related to the practice of classifying cosolvents and cosolutes as structure-makers and structure-breakers is one in which addition of cosolvents is thought to enhance or suppress the inherent hydrophobic properties of water.52-54 In this view, structuremakers are cosolvents that increase the rigidity and strength of the hydrogen bond network in water and therefore increase the (entropic) hydrophobic penalty for inserting small, apolar solutes.44,55 Historically, the classification kosmotrope has typically been considered synonymous with structure-making. Conversely, structure-breakers have historically been argued to disrupt the hydrogen bond network in water and therefore suppress the hydrophobic effect for small, apolar solutes.44,55,56 Finally, a recently developed viewpoint regarding hydrophobic solvation is one in which the bulk solvent equation of state and the intrinsic, local density fluctuations predominantly determine solvation thermodynamics of small, apolar solutes.57-60 In this view, although one may be able to identify certain hydrogen bond network patterns surrounding a solute (component 2), µex 2 is dictated to a large degree by the reversible work of forming a cavity (wcav(σHS)) with a diameter σHS equal to the effective hard core diameter for the solute at a given temperature and pressure. wcav(σHS) is determined solely by the intrinsic probability (p0(σHS)) of spontaneously observing a (spherical) region of diameter σHS that is devoid of any solvent molecules,30,60,61 and p0(σHS) is in turn determined by the equation of state of the solvent, i.e., its density at a given temperature and pressure, as well as its isothermal compressibility.60,61 The equation-of-state (EoS) approach also has important thermodynamic implications beyond what is outlined above. Continuing an EoS-based analysis, one finds that a number of the characteristic experimental signatures of the hydrophobic effectse.g., solubility minima vs temperature along isobarssare intimately related to the presence of thermodynamic anomalies in pure water such as density maxima vs temperature along isobars.62

Solvation in Water-Methanol and -Sorbitol Mixtures TABLE 1: System Specifications

TABLE 2: Lennard-Jones Parameters and Partial Charges

water molecules

cosolvent molecules

densities (g/cm3)

water

250

0

0.94-1.08

4 w/w% 8 w/w% 16 w/w% 30 w/w% 60 w/w%

240 228 204 168 84

1 2 4 7 14

0.96-1.06 0.98-1.08 0.98-1.08 1.06-1.16 1.30-1.35

6 12 24 42 84

0.94-1.04 0.92-1.02 0.88-0.98 0.82-0.92 0.72-0.82

system

sorbitol

methanol 4 w/w% 8 w/w% 16 w/w% 30 w/w% 60 w/w%

240 228 204 168 84

J. Phys. Chem. B, Vol. 111, No. 17, 2007 4469

Although the EoS perspective of apolar solvation has thus far primarily been employed with essentially pure water as the solvent, we include it in the above discussion because its basic tenets will prove useful in the discussion of our results in Section III. Specifically, we find that analysis in the spirit of the EoS view of water provides a means to quantitatively collapse ∆µtr2 for a given apolar solute across a range of conditions on a single, common curve that is independent of the temperature, density, cosolvent composition, or cosolvent type (methanol or sorbitol). II. Methods Molecular Dynamics Simulations. Liquid water along with binary mixtures of water and either sorbitol or methanol were simulated at selected compositions (0, 4, 8, 16, 30, and 60 w/w% cosolvent) using canonical ensemble molecular dynamics (MD). Periodic boundary conditions63 were employed with system sizes of approximately 250 water molecules or united-atom units (i.e., H2O and OH or CHn)1,2,3 groups). The number of molecules of water and each cosolvent type, and the range of densities for each composition, are shown in Table 1. MTK dynamics,64 an extension of the original Nose-Hoover dynamics,65 were employed to maintain temperature while conserving a modified Hamiltonian. A multiple time step algorithm was used to handle the higher frequency bond stretching, angle bending, and torsional modes.66 The equations of motion were integrated using the velocity verlet algorithm employing a 0.2 fs short time step (for intramolecular motions and interactions only) and a 2 fs long time step. Long-range electrostatic interactions were treated with the reaction field method67 as done previously for monosaccharide-water simulations.68,69 All interactions were pairwise with a cutoff distance of 7.9 Å. Model Potentials. Water was modeled with a flexible SPC/E potential,70-72 which treats the molecule as a single LennardJones (LJ) sphere centered on the oxygen atom, with three embedded partial atomic chargesstwo equal, partial positive charges on the center of the hydrogen atoms and a single, partial negative charge on the oxygen center. The parameter values72 are given in Table 2. Flexibility was added through harmonic oscillators of the form E ) k(r - req)2 for the H-O-H bond angle as well as O-H bond stretch.71 The intramolecular parameters71 for the spring constants and for equilibrium bond distances and angles are shown in Table 3. A united-atom approach was used to model methanol and sorbitol. Specifically, carbon and oxygen atoms are each given a LJ center. All hydrogens lie within the repulsive core of the respective carbon or oxygen centers. Hydroxyl H atoms possess an explicit partial atomic charge, as do the carbon and oxygen

atom type

 (kJ/mol)

σ (Å)

q (units of e)

Owater Hwater CH CH2 CH3 O H

0.65017

3.1656

0.3344 0.49324 0.86526 0.7106

3.85 3.905 3.775 3.07

-0.8476 +0.4238 +0.265 +0.265 +0.265 -0.700 +0.435

TABLE 3: Intramolecular Potential Parameters bond stretch

k (kcal/mol Å2)

req (Å)

O-H (water) CH3-OH CH2-OH CH-OH CH-CH2 CH-CH O-H

1109 386 386 386 260 260 553

1.00 1.425 1.425 1.425 1.526 1.526 0.96

bond angle

k (kcal/mol rad2)

θeq (o)

H-O-H CH3-O-H CH2-O-H CH-O-H CH-CH-OH CH-CH2-OH CH2-CH-OH CH2-CH-CH CH-CH-CH

92 55 55 55 80 80 80 63 63

109.47 108.5 108.5 108.5 109.5 109.5 109.5 111.5 111.5

torsional angle

n

Vn/2 (kcal/mol)

OH-CH2-CH-OH

2 3 2 3

0.5 1.0 0.5 0.5

OH-CH-CH-OH

atoms, but no LJ parameters. The parameter values were taken from the OPLS model73 and are given in Table 2. Intramolecular interactions were treated with harmonic oscillators for bond stretching and angle bending modes74 analogous to that described above for flexible SPC/E water. In addition to flexible bond lengths and angles, torsional potentials of the form E ) ∑n)2,3 (Vn/2)(1 + cos(nφ)) were applied in order to maintain the stereochemistry of the sorbitol chains. All intramolecular parameters are given in Table 3.74 Although the OPLS force field was originally optimized with the TIP4P model of water,73,75 SPC/E was employed here for convenience and for comparison with previous results from one of us.68 From a qualitative and semiquantitative perspective, SPC/E and all other commonly employed, classical potentials for liquid water incorporate the same physics of water-water and water-cosolvent interactions, and each produces water’s characteristic solvation thermodynamics and structure to a similar degree.76,77 As the analyses and conclusions in this report do not rely upon direct quantitative comparison to experimental data, quantitative disparities that may arise from using SPC/E rather than TIP4P are not expected to be significant. Simulation Procedure. In order to avoid numerical instabilities due to particle overlaps or excessive bond strains in initial configurations with high concentrations of (sorbitol) cosolvents, each system was initially prepared as a rarified arrangement at 400 K and a density of approximately 0.1 g/cm3. The density of the system was successively incremented in units of 0.01 g/cm3 by adjusting the molecular center-of-mass positions. Simulations were run for periods of 20 ps to allow temperature, potential energy, and pressure to stabilize between density increments. After reaching the desired density, an additional

4470 J. Phys. Chem. B, Vol. 111, No. 17, 2007 equilibration period of at least 200 ps was performed at 400 K to relax the system and remove artifacts of the densification procedure. The resulting configurations for each density at each species composition were used as the initial configurations for simulations performed at 400 K. Although this procedure was practically necessary for only the concentrated sorbitol-water solutions, it was applied to all systems for the sake of consistency. For each state point, an MD simulation was first equilibrated until the mean squared displacement of the less mobile species was greater than 100 Å2, or until 500 ps had been simulated, provided the mean squared displacement was at least 10 Å2. In general, these equilibration periods were sufficiently long for the energy and pressure to reach steady values (data not shown). Ensemble averaging was then performed for a period of 1 ns with configurations being saved every 0.5 ps. The final configuration from one simulation was used as the initial configuration for the next one at the same density and composition, but 20 K lower in temperature. The procedure of reequilibration and ensemble averaging continued for temperatures stepped sequentially to as low as 240 K in 20 K increments. For 30 and 60 w/w% sorbitol, initial relaxation took place at 800 K and simulations were performed to temperatures as low as 240 K, in 40 K increments, in order to span a larger range of average energies and a range of water diffusivities comparable to more dilute solutions. The first equilibration criterion, based on the displacement of the least mobile species, dominated at higher temperatures. The second criterion of a fixed equilibration period dominated at lower temperatures. Due to relatively low mobilities at lower temperatures, equilibration criteria were not met below 280 K for 30 w/w% sorbitol and below 320 K for 60 w/w% sorbitol. No state points were included in the subsequent analysis if they did not meet at least one of the equilibration criteria. Data Analysis. For each simulated state point, average virial and ideal components of the pressure were obtained using standard methods.63 Widom insertion78 was used to obtain excess chemical potential results for methane (LJ) and hard sphere (HS) solutes through 103 randomly placed insertions per saved configuration at each simulated state point. Details regarding statistical uncertainties and convergence of the average excess chemical potential values are provided as part of the Supporting Information. The values of µex 2 and average pressure (〈p〉) were each interpolated versus mass density (Fm) at fixed temperature (T) and cosolvent weight fraction (w3). A quadratic fit was used for 〈p〉 vs Fm, and the chemical potential was essentially linear in Fm (data not shown). This procedure yielded 〈p〉 (Fm;T,w3) and µex 2 (Fm;T,w3) for each T and w3; 2 ) HS or LJ solute, and the semi-colons in 〈p〉 (Fm;T,w3) ex and µex 2 (Fm;T,w3) indicate T and w3 were held fixed. 〈p〉 or µ2 was interpolated versus Fm. Inverting 〈p〉(Fm;T,w3) provided Fm(〈p〉;T,w3). µex 2 (〈p〉;T,w3) was then obtained parametrically from Fm(〈p〉;T,w3) and µex 2 (Fm;T,w3). Subsequent analyses, e.g., of the behavior of solute excess chemical potential and ∆µtr2 , were performed along isobars at fixed temperature as a function of composition, mass density (irrespective of cosolvent type), as well as packing fraction. For a given state point, an average packing fraction (η) was calculated using the LJ diameters from Table 1 as the size of the effective hard core of each united atom or water molecule. Monte Carlo integration provided the fractional volume in a given configuration that lay within the effective hard cores of the constituent molecules. This fractional volume was averaged over a large number of configurations (∼103) for each state point

Shah et al. to provide η(Fm;T,w3). These values were then interpolated and combined parametrically to provide η(〈p〉;T,w3) in the same manner described above for the excess chemical potential. Hard sphere rdfs for use in Kirkwood-Buff integrals were calculated from the positions of the hypothetical HS solutes that were accepted during Widom insertions across all saved configurations at a given state point. The rdfs were calculated within a subvolume of the simulation box that was smaller than the box length but as large as possible to allow the rdfs to approach unity at large separations. As noted previously,21,30 this procedure is adopted in order to allow particle numbers to fluctuate within the subvolume and thereby attempt to approximate behavior in the grand canonical ensemble while actually simulating in the computationally more expedient canonical ensemble. As it cannot be assumed a priori that the rdfs will converge sufficiently close to unity at long pairwise distances within the finite cutoffs needed, we calculated truncated Kirkwood-Buff integrals (G ˜ ij) as a function of cutoff distance (Rc),

G ˜ ij(Rc) ) 4π

R [gjij(r) - 1]r2 dr ∫r)0 c

(6)

The corresponding truncated preferential interaction parameter (γ32) is defined as

γ32 ) F3(G ˜ 23 - G ˜ 21)

(7)

In the limit of large Rc values, gjij(r) will reach unity, and G ˜ ij will equal Gij. Similarly, γ32 approaches Γ32 as Rc becomes sufficiently large. Hard Sphere Insertions. Widom test insertion of HS solutes of diameter σHS were performed with the same sets of configurations as those used in test insertions of methane. However, rather than neglecting the repulsive core of each H2O, CHn, and OH group, each such group was treated as having a HS with a diameter corresponding to the first maximum in the respective self-self rdf. This procedure was used to more realistically mimic the distribution of available space when compared to inserting soft-body solutes (e.g., methane) into configurations using the full (soft-body) potentials for H2O, CHn, and OH groups. As a result, the size of HS solutes that could be practically accommodated was significantly lower than if the more common practice was employed in which an interaction between an HS solute and a soft-body particle is considered zero until the center of the soft body is within 0.5σHS of the HS center.30,61 Therefore, the effective HS diameter when expressing µex HS in terms of the work of cavity formation depends on which solvent component one is considering (e.g., for σHS ) 2.8 Å, we have σeff,OH ) σeff,H2O ) 5.6 Å and σeff,CHn ) 6.4 Å). III. Results and Discussion Interpolated densities for each temperature and solvent composition along an atmospheric pressure isobar are shown in Figure 1. Flexible SPC/E water exhibits a density maximum near 260 K. This temperature is significantly below the experimental value of 277 K but is in good agreement with values for the temperature of maximum density (TMD) obtained by others using the (rigid) SPC/E model of water.68,79-81 Upon addition of methanol, the TMD shifts systematically to lower temperatures, below the lowest temperature reported here (Figure 1a). Addition of sorbitol also suppresses the TMD (Figure 1b), albeit to a somewhat lesser extent for a given w/w% cosolvent. The hydrophobic effect, as indicated by the temperature of maximum βµex (where β ) 1/kBT) or minimum

Solvation in Water-Methanol and -Sorbitol Mixtures

Figure 1. Density as a function of temperature for (A) methanolwater and (B) sorbitol-water mixtures. Symbols represent different compositions: 0 (circles), 4 (squares), 8 (diamonds), 16 (×), 30 (+), and 60 w/w % (triangles). Note that data for water (circles) are the same in both panels and can be seen more clearly in panel B.

solubility for both methane (Figure 2) and hard sphere solutes (Figure 3), is also suppressed upon addition of either methanol or sorbitol. Uncertainties in βµex are no greater than 5% of the average value reported, based on the standard deviation of the distribution of e-β∆E (∆E ) insertion energy) across insertions at a given state point (Supporting Information). The excess chemical potential results for methane and hard spheres in pure water agree quantitatively or semiquantitatively with experiment82 and previous simulation61,83-85 once differences in solvent density are accounted for. For example, previous studies84,85 with SPC/E water report an excess chemical potential for methane between 8.5 and 10 kJ/mol at atmospheric pressure (1.00 ( 0.02 g/cm3) and 300 K. At this temperature and bulk density, our results give an interpolated value for µex of 8.1 ( 0.4 kJ/mol. The value of 8.1 kJ/mol just noted is different from the value shown at 300 K in Figure 2A because of the difference

J. Phys. Chem. B, Vol. 111, No. 17, 2007 4471

Figure 2. Excess chemical potential of methane as a function of temperature in (A) methanol-water and (B) sorbitol-water mixtures. Symbols represent different compositions: 0 (circles), 4 (squares), 8 (diamonds), 16 (×), 30 (+), and 60 w/w % (triangles). The lines serve as a guide to the eye.

in liquid densities between the previous and current work at atmospheric pressure. Within statistical uncertainty, our results are also in agreement with the average experimental value of 8.1 kJ/mol at atmospheric pressure and 298 K (experimental density 0.997 g/cm3).82 Hard sphere µex 2 values in pure water agree with previous simulation data61,83 up to the largest attempted size of σeff,H2O ) 7.6 Å (Supporting Information). To the best of our knowledge, no comparable results for methane or hard sphere solvation in sorbitol-water and/or methanolwater mixtures have been reported under conditions sufficiently similar to those examined here. Figures 1-3 indicate that addition of either sorbitol or methanol suppresses canonical signatures of the hydrophobic effect and density anomalies to lower temperatures. This behavior argues against the hydrophobic-enhancement or -suppression view for the case of these cosolvents and by inference argues against that view for cosolvents with similar chemistries such as other polyols and primary alcohols.

4472 J. Phys. Chem. B, Vol. 111, No. 17, 2007

Figure 3. Excess chemical potential of hard spheres (σHS ) 1.6 Å, σeff,OH ) σeff,H2O ) 4.4 Å, and σeff,CHn ) 5.2 Å) as a function of temperature along 1 atm isobar in (A) methanol-water and (B) sorbitol-water mixtures. Symbols represent different compositions: 0 (circles), 4 (squares), 8 (diamonds), 16 (×), 30 (+), and 60 w/w % (triangles). The lines serve as a guide to the eye.

With the above behaviors in mind, it is worth noting that one of the original experimental criteria for classifying a cosolvent as “structure-making” is that it increases the TMD,52,54 presumably via strengthening the three-dimensional hydrogen bond network that underlies the existence of density maxima and the location of the TMD locus for pure water.86,87 Within that context, it is obvious that both methanol and sorbitol would be classified as “structure-breakers” despite the fact that polyols are themselves hydrogen bond network formers and are highly hydrogen bonded in solution.68,88,89 As shown later in this section, despite sorbitol suppressing the hydrophobic signatures and the TMD of water, its addition clearly results in positive ∆µtr2 and negative Γ32. Methanol is even more complex in that its addition can result in either positive or negative ∆µtr2 and Γ32. Together, these results lead to the conclusion that even the qualitative nature of apolar solvation in this class of mixed aqueous solvents is not

Shah et al. accurately accounted for by arguments based on perspectives such as structure-making/structure-breaking or enhancement/ suppression of the hydrophobic effect. Figure 4 shows the transfer free energies (eq 1) for methane and HS solutes as a function of cosolvent composition (sorbitol, Figure 4A; methanol, Figure 4B) at fixed temperature and atmospheric pressure. The results in Figure 4C are plotted for ∆µtr2 per unit water accessible surface area (ASA) in order to place them on a more convenient scale for comparison to experimental data with a larger solute (argon) in ethanol-water solutions at atmospheric pressure and 280 K. The use of a per unit ASA basis is purely for convenience and is not intended to indicate any mechanistic interpretations that can be inferred from the scaling dependence of solute chemical potentials with solute dimensions.90-92 The results in Figure 4 illustrate a number of contrasting features of apolar solvation in sorbitol-water versus methanolwater solutions. Sorbitol-water solutions exhibit near-linear, monotonically increasing ∆µtr2 with increasing sorbitol concentration for both methane and HS solutes. The ∆µtr2 vs m3 profiles are only weakly dependent on temperature, implying a relatively small difference in excess partial molar entropy upon transfer from water to sorbitol-water solvents relative to methanol-water solvents. Note, however, this should not be simply interpreted as indicating solute-water and solutecosolvent interactions dominate ∆µtr2 ; the results later in this section and in the Supporting Information clearly indicate this is not the case. Rather, this behavior is due to changes in enthalpy caused by disruption of favorable water-water, watercosolvent, and cosolvent-cosolvent interactions upon cavity formation. In contrast, the ∆µtr2 vs m3 profiles for methane and HS solutes in methanol-water solvents are significantly temperature dependent and also display qualitatively different composition dependence compared to sorbitol-water. Specifically, ∆µtr2 exhibits a maximum as a function of m3 at low temperatures; with increasing temperature, this maximum is suppressed until only a monotonic, near-linear decrease in ∆µtr2 with increasing m3 is observed at high temperatures. For methanol-water solutions, both enthalpic and entropic contributions to ∆µtr2 are appreciable depending on which temperature is considered (not shown). In terms of experimental observables, the results in Figure 4 indicate that addition of sorbitol disfavors dissolution or solvation of small apolar solutes at essentially all experimentally accessible sorbitol concentrations and across a broad range of temperatures. Addition of methanol, however, disfavors apolar solvation at low methanol concentrations and temperatures but favors solvation at higher concentrations or temperatures. It is interesting to note that the behavior for methanol-water systems is qualitatively and semiquantitatively in agreement with that observed experimentally for Ar solubility in ethanol-water solutions (Figure 4C).39 As mentioned in Section I, Kirkwood-Buff integrals can in principle be used with the pair distribution functions to calculate Γ32, which in turn determines the sign of (∂∆µtr2 /∂m3)T,p (eq 3). Due to poor statistics for rdfs between the water/cosolvent molecules and larger diameter hard spheres, the results for a smaller hard sphere (σHS ) 0.6 Å, σeff,OH ) σeff,H2O ) 3.4 Å, and σeff,CHn ) 4.2 Å) are used to illustrate the behavior. The rdfs for larger σHS values exhibited the same qualitative behavior for all cases shown. Figure 5A shows the solute-water and solute-cosolvent rdfs for a σHS ) 0.6 Å hard sphere solute in methanol-water

Solvation in Water-Methanol and -Sorbitol Mixtures

J. Phys. Chem. B, Vol. 111, No. 17, 2007 4473

Figure 5. (A) Radial distribution functions and (B) corresponding truncated interaction parameters (γ32) as a function of cutoff distance (Rc) for hard spheres (σHS ) 0.6 Å, σeff,OH ) σeff,H2O ) 3.4 Å, and σeff,CHn ) 4.2 Å) in methanol-water mixtures as a function of composition at 320 K and approximately atmospheric pressure. Arrows indicate increasing composition of cosolvent (0, 4, 8, 16, 30, 60 w/w%). Inset shows γ32 evaluated to longer distances to show finite size effects. Note that no HS-cosolvent rdf or γ32 can be defined for pure water.

Figure 4. Transfer free energy of methane (connected by solid lines) and hard spheres (connected by dotted lines) as a function of (A) sorbitol and (B) methanol composition. σHS is the same as in Figure 3. Overlay of experimental data39 for argon solvation in ethanol-water mixtures at 280 K (solid line) with simulation results for methane solvation in methanol-water mixtures. x-Axis in (C) provided as mole fraction of CHn + OH groups for easier comparison between alcohols of different molecular weights. Symbols represent different temperatures: 360 (circles), 340 (squares), 320 (diamonds), 300 (×), 280 (+), and 260 K (triangles) at 1 atm. The lines serve as a guide to the eye.

mixtures from simulations at 320 K and densities corresponding to an approximately atmospheric isobar. Relative to water, methanol is displaced from the immediate surface of the solute. The solute-sorbitol rdf is significantly below one for the first few solvation shells (Supporting Information). Quantitatively, methanol molecules are displaced to much shorter distances within the first solvation shell compared to sorbitol molecules. This behavior is qualitatively anticipated simply on the basis of the larger excluded volume of sorbitol relative to methanol. Overall, this agrees with qualitative arguments relating configurational entropy to the larger excluded volumes of both cosolvent species relative to water.93-95 However, in order to accurately interpret cosolvent behavior and thermodynamics of

4474 J. Phys. Chem. B, Vol. 111, No. 17, 2007 solvation in the framework of preferential accumulation or exclusion (i.e., KB theory), one must combine and integrate the rdfs for both cosolvent-solute and water-solute to produce γ32(Rc) or Γ32 (see also eqs 4 and 7 and associated discussion in Section I). An example of the truncated preferential interaction parameter, γ32 (eq 7), as a function of cutoff distance (Rc) is given in Figure 5B for a HS solute in methanol-water mixtures. Despite the displacement of methanol from the HS surface relative to displacement of water, and the lack of direct solute-cosolvent or solute-water interactions, the sum over the first few solvation shells leads to positive γ32 or preferential solvation by the cosolvent at high methanol concentrations. At intermediate methanol concentrations, an analysis based purely on γ32 is inconclusive because it is not straightforward to determine even the sign of the value to which γ32 will converge at large Rc. This is compounded by the well-known issues21,30 with convergence of γ32 at large Rc in finite-sized simulations, due to the r2 factor in the integrand of eqs 5 and 6, and the fact that rdfs do not practically achieve a value of exactly 1 in finite simulations.22,25,26 Similar results regarding ambiguity of the sign of Γ32 were also found for methane in methanol-water solutions in larger simulation boxes. For sorbitol-water solutions (Supporting Information), the sign of γ32 and Γ32 is less ambiguous, but a quantitative value of Γ32 still cannot be reliably extracted due to the above-noted convergence issues at large Rc. The preceding results indicate: (1) qualitative, but not quantitative, determination of Γ32 via molecular simulation is potentially possible for small apolar solutes in concentrated mixed-aqueous solvents, provided |Γ32| is significantly nonzero; (2) even though methanol-water solutions clearly display (d∆µtr2 /dm3)T,p < 0 for small apolar solutes at high m3, and thus Γ32 must be positive, it is incorrect to interpret this thermodynamic behavior as indicative of preferential accumulation in the immediate vicinity of the solutesrather, accumulation occurs only beyond the first one or two solvation shells. The last result illustrates that KB theory and all rigorous definitions of preferential exclusion and accumulation of cosolvents16,25,30,47,49,96 require one to consider correlations in solute and solvent composition over a large (ideally, a macroscopic) and open control volume, not just the first few solvation shells of a solute. Turning now to an alternative analysis of our results, the ∆µtr2 results from Figure 4 for methane in both sorbitol-water and methanol-water solvents over a range of temperatures are plotted in Figure 6A as a function of the reduced packing fraction, η/ηw, defined as follows. Each point in Figure 6A corresponds to a particular cosolvent composition, temperature, and density from Figure 4. The density for each point is that required to maintain atmospheric pressure at that temperature and cosolvent composition. The average packing fraction for that cosolvent state point (T, p ) 1 atm, w3 or m3) is denoted by η. The average packing fraction for pure water at that T and p ) 1 atm is denoted by ηw. What is perhaps most striking from Figure 6 is that the ∆µtr2 results from a broad range of temperatures and compositions collapse onto a common curve that is also independent of whether the cosolvent is methanol or sorbitol. This is particularly surprising when one considers that ∆µtr2 is monotonically increasing, with a weak temperature dependence for sorbitolwater, but is strongly temperature-dependent and non-monotonic for methanol-water solutions (cf. Figure 4). Note that different densities correspond to different cosolvent types and/or concentrations. At fixed composition, η is essentially linear in density over the relevant ranges for either sorbitol-water or

Shah et al.

Figure 6. (A) Transfer free energy of methane as a function of reduced packing fraction (defined in text). Different symbols denote different cosolvents and temperatures as indicated in the legend. Pure water is represented by an asterisk at (1,0). (B) Reduced packing fraction as a function of reduced mass density for methanol-water (open symbols) and sorbitol-water (closed symbols) along two selected isotherms (360 and 260 K). Symbols represent same conditions in panel A. Pure water is denoted by an asterisk at (1,1), and lines are guides to the eye.

methanol-water solutions (results not shown). However, changes in solvent composition along with the corresponding change in density along an isobar can lead to the packing fraction being nonlinear in density or composition. This behavior is evident in Figure 6B, where η/ηw is shown as a function of the reduced solvent mass density (Fm/Fm,water) for the highest and lowest temperature data from the main panel. Considering that η/ηw is not always a monotonic function of Fm (composition not fixed), we conclude that Fm/Fm,water is not a complete surrogate for η/ηw in Figure 6A. Comparison between Figure 4, Figure 6A, and Figure 6B indicates: (1) packing fraction along an isobar is a non-simple function of solvent density due to changes in solvent composition; (2) the temperature dependence of ∆µtr2 can be quantitatively accounted for by the temperature dependence of the

Solvation in Water-Methanol and -Sorbitol Mixtures solution density at the given pressure, so long as it is combined with the corresponding packing fraction η; (3) the qualitative and quantitative differences in ∆µtr2 as a function of composition across both sorbitol-water and methanol-water solvents can be accounted for by converting from a composition or density scale to one of packing fraction. In addition, in Figure 6A, one observes that positive ∆µtr2 (“salting-out” or unfavorable apolar solvation relative to water) requires η/ηw > 1, and negative ∆µtr2 (“salting-in” or favorable apolar solvation) requires η/ηw < 1. Figure 6B shows that η/ηw > 1 (“salting-out”) is always obtained for sorbitol-water solutions if the overall mass density is greater than that of pure water. For methanol-water solutions at high temperatures, a lower overall mass density relative to pure water leads to η/ηw < 1 (“salting-in”), but the relationship fails at lower temperatures. Thus, the overall mass density of the solvent as a function of composition at fixed T and p, i.e., its equation of state, is a primary but not necessarily the sole factor that determines whether apolar solvation will be favorable or unfavorable compared to pure water. Recalling the earlier analysis of the truncated preferential interaction parameters (Figure 5B), the behavior shown in Figure 6 is perhaps not unexpected to some degree because methanol is not strongly associated within the first neighbor solvation shell with either methane or HS solutes, even up to high methanol concentrations (60 w/w%). Thus, the mechanism by which methanol solubilizes or stabilizes apolar solutes in aqueous solutions is anticipated to be fundamentally different from classic, strong chemical denaturants and solubilizing agents such as urea.3,97 This analysis further highlights that caution should be exercised when interpreting positive Γ32, or equivalently (∂∆µtr2 /∂m3)T,p < 0, to indicate “binding” or strong association of the cosolvent molecules with the molecular surface of a solute. Rather, one must also consider the effects of the cosolvent on the overall solvent equation of state, as this may dominate over or effectively compete with excluded volume effects and direct solute-cosolvent interactions that are typically the primary focus of a number of models for Γ32.16-18,24 Finally, although not discussed in detail as part of our analysis here, it should be noted that the equation of state contributions to ∆µtr2 and Γ32 are inherently determined to a large extent by the binary water-cosolvent (or 1-2) interactions. Thus, another way to consider the results provided in this report is to conclude that accurate models of apolar solvation in mixed-aqueous solvents cannot rely solely upon accounting for direct solutecosolvent (or 2-3) interactions and excluded volume effects. They also must implicitly or explicitly treat the water-cosolvent interactions in an accurate manner to avoid neglecting the effects of those interactions on both overall density and the packing arrangements inherent in the binary solvent. For much larger solutes than those considered here, the latter effects may arguably be of lesser importance, but the former effects will remain essential. IV. Summary and Conclusions Using MD simulations, we studied the equation of state and solvation of small apolar solutes (methane and hard spheres) in methanol-water and sorbitol-water systems as a function of composition, density (pressure), and temperature. The presence of either cosolvent results in a number of important changes in the solvent properties, including a shift to lower temperatures for the TMD and the temperature of minimum solubility for apolar gases. Both of these results indicate a disruption of characteristic properties of water as a solvent, i.e., the hydro-

J. Phys. Chem. B, Vol. 111, No. 17, 2007 4475 phobic effect. Combined with the facts that addition of sorbitol disfavors apolar solute dissolution, and methanol can either favor or disfavor dissolution, our results indicate classification schemes such as structure-making/-breaking and hydrophobic enhancement/suppression are of questionable utility for these and chemically similar types of water-cosolvent mixtures. Kirkwood-Buff analysis based on simulated radial distribution functions clearly indicates that, despite Γ32 > 0 being required at high methanol concentrations for thermodynamic consistency, it is incorrect to interpret positive Γ32 values as necessitating preferential accumulation of methanol in the first few solvation shells of small apolar solutessrather, methanol is preferentially accumulated only at longer distances from the solute surface. At lower methanol concentrations, KB analysis is not reliable due to |Γ32| not being sufficiently greater than zero. Similar qualitative and quantitative limitations are found with sorbitol-water solvents. In contrast to the limited success of the above approaches to interpret apolar solvation thermodynamics in methanol-water and sorbitol-water mixtures, accounting for the EoS and local packing in the binary mixtures provides a means to quantitatively describe transfer free energies for a variety of conditions (T, p, composition, and cosolvent type) in terms of a common profile that requires only knowledge of the overall packing fraction at a given state point and in pure water. This result is particularly striking in light of the strong, nonmonotonic dependence of ∆µtr2 on methanol composition at low temperatures. For sorbitol-water solutions, packing fraction is a monotonic function of density, across different compositions, over the range of conditions studied here. Therefore, for sorbitol-water solutions the EoS contribution dominates and essentially dictates the sign and magnitude of ∆µtr2 for apolar solutes. Methanolwater solutions at high temperatures display behavior analogous to that for sorbitol-water solutions. However, at low temperatures packing fraction is a non-monotonic function of solvent density for methanol-water systems, and so both the EoS and the local packing or solvent structure play a role in dictating the sign and magnitude of ∆µtr2 . On the basis of the observation that the majority of polyols and carbohydrates can be reasonably approximated as primarily composed of methanol subunits, the results reported here suggest similar findings may hold across a range of naturally and industrially relevant cosolvents that are derived from these chemical building blocks. Acknowledgment. The authors gratefully acknowledge the donors of the American Chemical Society Petroleum Research Fund and the University of Delaware Research Foundation for support of this research. Supporting Information Available: Additional details are available regarding Widom insertion and estimates of uncertainty, comparison of HS µex with literature data, and analysis of incomplete preferential interaction parameters for sorbitolwater and methanol-water mixtures. This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Gould, P.; Goodman, M.; Hanson, P. Int. J. Pharm. 1984, 19, 149. (2) Ruelle, P.; Kesselring, U. W. J. Pharm. Sci. 1998, 87, 998. (3) Timasheff, S. N. Annu. ReV. Biophys. Biomol. Struct. 1993, 22, 67. (4) Shulgin, I.; Ruckenstein, E. J. Chem. Phys. 2005, 123. (5) Ruckenstein, E.; Shulgin, I. Int. J. Pharm. 2003, 267, 121.

4476 J. Phys. Chem. B, Vol. 111, No. 17, 2007 (6) Blandamer, M. J.; Engberts, J.; Gleeson, P. T.; Reis, J. C. R. Chem. Soc. ReV. 2005, 34, 440. (7) Slade, L.; Levine, H. Crit. ReV. Food Sci. Nutr. 1991, 30, 115. (8) Chirife, J.; Buera, M. D. Crit. ReV. Food Sci. Nutr. 1996, 36, 465. (9) Beuchat, L. R. J. Food Prot. 1983, 46, 135. (10) Sperber, W. H. J. Food Prot. 1983, 46, 142. (11) Shulgin, I.; Ruckenstein, E. Ind. Eng. Chem. Res. 2002, 41, 6279. (12) Campanella, E. A.; Mathias, P. M.; Oconnell, J. P. AIChE J. 1987, 33, 2057. (13) McKetta, J. J.; Katz, D. L. Ind. Eng. Chem. 1948, 40, 853. (14) Dhima, A.; de Hemptinne, J. C.; Moracchini, G. Fluid Phase Equilib. 1998, 145, 129. (15) Dhima, A.; de Hemptinne, J. C.; Jose, J. Ind. Eng. Chem. Res. 1999, 38, 3144. (16) Schellman, J. A. Biopolymers 1978, 17, 1305. (17) Schellman, J. A. Biophys. J. 2003, 85, 108. (18) Zimmerman, S.; Minton, A. Annu. ReV. Biophys. Biomol. Struct. 1993, 22, 27. (19) Minton, A. P. J. Pharm. Sci. 2005, 94, 1668. (20) Yaacobi, M.; Ben-Naim, A. J. Phys. Chem. 1974, 78, 175. (21) Ben-Naim, A. Water and Aqueous Solutions; Plenum: New York, 1974. (22) Ben-Naim, A. J. Chem. Phys. 1977, 67, 4884. (23) Ben-Naim, A. J. Phys. Chem. 1989, 93, 3809. (24) Parsegian, V. A.; Rand, R. P.; Rau, D. C. Proc. Natl. Acad. Sci. U.S.A. 2000, 97, 3987. (25) Chitra, R.; Smith, P. E. J. Phys. Chem. B 2001, 105, 11513. (26) Chitra, R.; Smith, P. E. J. Phys. Chem. B 2002, 106, 1491. (27) Baynes, B. M.; Trout, B. L. J. Phys. Chem. B 2003, 107, 14058. (28) Auton, M.; Bolen, D. W. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 15065. (29) Fredenslund, A.; Jones, R. L.; Prausnitz, J. M. AIChE J. 1975, 21, 1086. (30) Ben-Naim, A. Statistical Thermodynamics for Chemists and Biochemists; Plenum Press: New York, 1992. (31) Shulgin, I. L.; Ruckenstein, E. J. Phys. Chem. B 2006, 110, 12707. (32) Chi, E.; Krishnan, S.; Randolph, T.; Carpenter, J. Pharm. Res. 2003, 20, 1325. (33) Scatchard, G. J. Am. Chem. Soc. 1946, 68, 2315. (34) Sandler, S. I. Chemical & Engineering Thermodynamics, 3rd ed.; John Wiley & Sons, Inc.: New York, 1999. (35) O’Connell, J.; Prausnitz, J. Ind. Eng. Chem. Fundam. 1964, 3, 347. (36) Smith, A. V.; Hall, C. K. Proteins: Struct., Funct., Genet. 2001, 44, 344. (37) Fernandez, A.; Sinanoglu, O. Biophys. Chem. 1985, 21, 163. (38) Herskovits, T.; Gadebeku, B.; Jaillet, H. J. Biol. Chem. 1970, 245, 2588. (39) Brandts, J. F.; Hunt, L. J. Am. Chem. Soc. 1967, 89, 4826. (40) Ben-Naim, A.; Yaacobi, M. J. Phys. Chem. 1974, 78, 170. (41) Gekko, K. J. Biochem. (Tokyo) 1982, 91, 1197. (42) Gekko, K.; Morikawa, T. J. Biochem. (Tokyo) 1981, 90, 39. (43) Hribar, B.; Southall, N. T.; Vlachy, V.; Dill, K. A. J. Am. Chem. Soc. 2002, 124, 12302. (44) Collins, K.; Washabaugh, M. Quart. ReV. Biophys. 1985, 18, 323. (45) Xie, G. F.; Timasheff, S. N. Protein Sci. 1997, 6, 211. (46) Casassa, E. F.; Eisenberg, H. AdV. Protein Chem. 1964, 19, 287. (47) Anderson, C.; Courtenay, E.; Record, M. J. Phys. Chem. B 2002, 106, 418. (48) Debenedetti, P. G. Metastable Liquids: Concepts and Principles; Princeton University Press: Princeton, New Jersey, 1996. (49) Kirkwood, J. G.; Buff, F. P. J. Chem. Phys. 1951, 19, 774. (50) Timasheff, S. N. Biochemistry 1992, 31, 9857. (51) Timasheff, S. N. Biochemistry 2002, 41, 13473. (52) Frank, H. S.; Evans, M. W. J. Chem. Phys. 1945, 13, 507. (53) Oro, J. R. D. J. Biol. Phys. 2001, 27, 73. (54) Sarma, T.; Ahluwalia, J. C. J. Phys. Chem. 1970, 74, 3547.

Shah et al. (55) Dill, K.; Truskett, T.; Vlachy, V.; Hribar-Lee, B. Annu. ReV. Biophys. Biomol. Struct. 2005, 34, 173. (56) Batchelor, J. D.; Olteanu, A.; Tripathy, A.; Pielak, G. J. J. Am. Chem. Soc. 2004, 126, 1958. (57) Pratt, L. R. Annu. ReV. Phys. Chem. 2002, 53, 409. (58) Bishop, W. H. Biophys. Chem. 1987, 27, 197. (59) Berne, B. J. Proc. Natl. Acad. Sci. U.S.A. 1996, 93, 8800. (60) Hummer, G.; Garde, S.; Garcia, A.; Pohorille, A.; Pratt, L. Proc. Natl. Acad. Sci. U.S.A. 1996, 93, 8951. (61) Hummer, G.; Garde, S.; Garcia, A.; Paulaitis, M.; Pratt, L. J. Phys. Chem. B 1998, 102, 10469. (62) Ashbaugh, H.; Truskett, T.; Debenedetti, P. J. Chem. Phys. 2002, 116, 2907. (63) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, U.K., 1987. (64) Martyna, G.; Klein, M.; Tuckerman, M. J. Chem. Phys. 1992, 97, 2635. (65) Hoover, W. Phys. ReV. A 1985, 31, 1695. (66) Martyna, G.; Tuckerman, M.; Tobias, D.; Klein, M. Mol. Phys. 1996, 87, 1117. (67) Steinhauser, O. Mol. Phys. 1982, 45, 335. (68) Roberts, C.; Debenedetti, P. J. Phys. Chem. B 1999, 103, 7308. (69) Caffarena, E.; Grigera, J. J. Chem. Soc., Faraday Trans. 1996, 92, 2285. (70) Toukan, K.; Rahman, A. Phys. ReV. B 1985, 31, 2643. (71) Wallqvist, A.; Teleman, O. Mol. Phys. 1991, 74, 515. (72) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269. (73) Jorgensen, W. J. Phys. Chem. 1986, 90, 1276. (74) Weiner, S.; Kollman, P.; Case, D.; Singh, U.; Ghio, C.; Alagona, G.; Profeta, S.; Weiner, P. J. Am. Chem. Soc. 1984, 106, 765. (75) Jorgensen, W.; Chandrasekhar, J.; Madura, J.; Impey, R.; Klein, M. J. Chem. Phys. 1983, 79, 926. (76) Francesco Sciortino, S. S.; Poole, P. H. Molecular Dynamics Simulations of Water; World Scientific: New Jersey, 1995. (77) Guillot, B. J. Mol. Liq. 2002, 101, 219. (78) Widom, B. J. Chem. Phys. 1963, 39, 2808. (79) Harrington, S.; Poole, P. H.; Sciortino, F.; Stanley, H. E. J. Chem. Phys. 1997, 107, 7443. (80) Baez, L. A.; Clancy, P. J. Chem. Phys. 1994, 101, 9837. (81) Bryk, T.; Haymet, A. D. J. Mol. Simul. 2004, 30, 131. (82) Ben-Naim, A.; Marcus, Y. J. Chem. Phys. 1984, 81, 2016. (83) Postma, J. P. M.; Berendsen, H. J. C.; Haak, J. R. Faraday Symp. Chem. Soc. 1982, 55. (84) Paschek, D. J. Chem. Phys. 2004, 120, 6674. (85) Konrad, O.; Lankau, T. J. Phys. Chem. B 2005, 109, 23596. (86) Stanley, H.; Barbosa, M.; Mossa, S.; Netz, P.; Sciortino, F.; Starr, F.; Yamada, M. Phys. A 2002, 315, 281. (87) Errington, J. R.; Debenedetti, P. G. Nature (London) 2001, 409, 318. (88) Barone, G.; Castronuovo, G.; Doucas, D.; Elia, V.; Mattia, C. A. J. Phys. Chem. 1983, 87, 1931. (89) Franks, F.; Grigera, J. R. Water Sci. ReV. 1990, 5, 187. (90) Lum, K.; Chandler, D.; Weeks, J. D. J. Phys. Chem. B 1999, 103, 4570. (91) Huang, D. M.; Geissler, P. L.; Chandler, D. J. Phys. Chem. B 2001, 105, 6704. (92) Ashbaugh, H. S.; Kaler, E. W.; Paulaitis, M. E. J. Am. Chem. Soc. 1999, 121, 9243. (93) Bhat, R.; Timasheff, S. N. Protein Sci. 1992, 1, 1133. (94) Arakawa, T.; Timasheff, S. N. Biochemistry 1985, 24, 6756. (95) Schachman, H. K.; Lauffer, M. A. J. Am. Chem. Soc. 1949, 71, 536. (96) Schellman, J. A. Biopolymers 1975, 14, 999. (97) Bennion, B.; Daggett, V. Proc. Natl. Acad. Sci. U.S.A. 2003, 100, 5142.