INTRODUCTION

Top Abstract Introduction Results & Discussion Conclusions References

The thermodynamics of self-assembly in aqueous solution is governed in large part by the microscopic structure and organization of water around solutes having distinct chemical architectures. Perhaps the simplest self-assembly process driven by solvent-mediated forces is the pairwise association of methane in water. Explicit water simulations have been applied extensively to characterize the methane-methane potential of mean force in water (New and Berne, 1995; Pangali et al., 1979; Smith and Haymet, 1993; van Belle and Wodak, 1993; Young and Brooks, 1997), which in turn has been useful for inferring the influence of hydrophobic interactions on the self-assembly of more complicated solutes. The effects of temperature (Lüdemann et al., 1996, 1997; Skipper et al., 1996), pressure (Hummer et al., 1998; Payne et al., 1997), and the shape of simple, idealized hydrophobic solutes (Garde et al., 1996; Wallqvist and Berne, 1995a,b) on hydrophobic interactions have been studied with this approach to provide insights into the molecular mechanisms of, for example, protein folding (Kauzmann, 1959) or surfactant self-assembly (Hunter, 1987; Israelachvili, 1992) in aqueous solution. Explicit water simulations of the hydration of simple n-alkanes (Garde et al., 1996; Jorgensen, 1982; Jorgensen and Buckner, 1987; Kaminski et al., 1994; Rosenberg et al., 1982; Tobias and Brooks, 1990; Wallqvist and Covell, 1995, 1996) have likewise proven useful for examining the influence of hydrophobic interactions in stabilizing certain molecular conformations in aqueous solution.

For applications involving macromolecular solutes, explicit water simulations have only limited utility because of the large number of waters of hydration and the large number of solute internal degrees of freedom that must be taken into account. The statistical precision required for these simulations places inordinate demands on computational resources and limits the system sizes that can be realized. An alternative approach is to use implicit water models, such as those embodied in continuum models of hydration. Solvent-mediated driving forces are accounted for within the context of the continuum model, but the inherent differences between water organization around hydrophobic and hydrophilic solutes lead to disconnected theoretical treatments of hydrophobic and hydrophilic hydration. For example, the hydrophobic contribution to the free energy of hydration is usually taken to be proportional to solute surface area (Chothia, 1974; Hermann, 1972; Reynolds et al., 1974), although the definition of surface area is not unique (Lee and Richards, 1971; Richards, 1977). Electrostatic contributions, on the other hand, are calculated assuming the solvent responds as a macroscopic dielectric medium on all length scales (Honig and Nicholls, 1995; Jackson, 1975; Nakamura, 1996; Rogers, 1986). Consequently, phenomenological parameters contained in the continuum model must be fit to selected experimental data. Alanine dipeptide conformational equilibria (Marrone et al., 1996; Ösapay et al., 1996; Schmidt and Fine, 1994) and self-assembly processes, such as -helix propagation (Yang and Honig, 1995a) and -sheet formation (Yang and Honig, 1995b), have been described successfully by using the continuum model with parameters fit to the solubilities of small nonpolar and polar solutes in water (Schmidt and Fine, 1994; Sitkoff et al., 1994, 1996). However, the physical significance of these parameters is not evident, and the consequences of applying these parameters to the description of complex self-assembly phenomena are uncertain.

Continuum model parameters can also be fit to explicit water simulations as a means of evaluating the intrinsic macroscopic approximations to "microscopic reality" in the continuum model. For example, a consequence of treating water as a macroscopic dielectric medium is that the solvent-induced electrostatic potential varies linearly with charge on the solute (Åqvist and Hansson, 1996; Figueirido et al., 1994; Pratt et al., 1994). Thus the electrostatic contribution to the hydration free energy is proportional to the square of solute charge. Extensive simulations of both ionic and polar solutes in explicit water have been carried out to test the validity of the continuum approximation, which for ion hydration is found to hold over a wide range of ionic charge (Garde et al., 1998; Hummer et al., 1996b; Jayaram et al., 1989; Kalko et al., 1996; Straatsma and Berendsen, 1988). In contrast, electrostatic interactions between polar solutes and water frequently do not exhibit a linear response (Åqvist and Hansson, 1996). In particular, when the polar solute is water itself, the solvent response is distinctly nonlinear (Hummer et al., 1995; Rick and Berne, 1994).

The central question we address in this paper is: Can the continuum model describe simultaneously the vacuum-to-water transfer free energies of simple hydrophobic and amphiphilic solutes, as well as their conformational equilibria in aqueous solution? To this end, we consider idealized hydrophobic and hydrophilic solutes at infinite dilution in water and calculate their free energies of transfer from vacuum to water, as well as the change in free energy of hydration for well-defined conformational transitions. The normal alkanes serve as benchmarks for hydrophobic hydration. We consider as well n-butane with a charge of ±1e (e is the fundamental unit of charge) added to one terminal methyl group. These ionic tetramers were selected for two reasons. First, the effect of electrostatic interactions with water can be separated from the underlying hydrophobic contributions, because the hydration thermodynamics of n-butane has been studied extensively by simulation. Second, these solutes represent idealized models of amphiphilic molecules for which we can adjust unambiguously the balance between hydrophobic and hydrophilic interactions.

	FREE ENERGY PERTURBATION USING EXPLICIT WATER SIMULATIONS

The free energy of hydration, A_hyd, for a monatomic solute is defined as the free energy of transferring the solute from vacuum to water (Ben-Naim, 1978; Ben-Naim and Marcus, 1984). For molecular solutes, A_hyd also depends on solute conformation, which for linear alkane chains we define by the backbone dihedral angles,  = (₁, ₂, ... , _n). The probability of observing a specific solute conformation in water is dictated by A_hyd() and the intramolecular energy of this solute conformation in the ideal gas, E_int():

(1)

where ¹ = k_BT. The normalized distribution of conformations depends only on relative energy differences and not the absolute value of either A_hyd() or E_int(). In this work, we neglect the intramolecular contributions and focus only on the hydration contribution; i.e., A_hyd() as a function of .

The n-alkane hydration free energies were calculated as the sum of hydration free energies for incrementally transforming C_n1 into C_n,

(2)

where n is the hydrocarbon chain length and C₀ denotes pure water. The A_hyd(C_i1  C_i) values are determined from the free energy perturbation (FEP) expression (Zwanzig, 1954),

(3)

where  refers to the reference state and  +  to the perturbed state of the system, E() is the potential energy in state , and the brackets ... denote averaging over solvent configurations in state . For the transformation C_i1  C_i, the Lennard-Jones (LJ) interaction parameters,  and , of the individual methyl groups are scaled linearly with . This involves "growing in" the ith methyl group at a chain end and simultaneously transforming the LJ parameters of the existing groups as required by the OPLS parameterization described below, i.e.,  = _final + (1  )_initial and  = _final + (1  )_initial. Perturbations of  = ±0.05 were performed at fixed values of  from 0.05 to 0.95 in increments of 0.10. A_hyd(C_i1  C_i) is then calculated as the sum of A_hyd(   + ) from  = 0 to 1. These calculations were carried out for n-butane, n-pentane, and n-hexane in the all-trans conformation ( = 180°).

The free energy of hydration as a function of conformation for n-butane and the charged tetramers was calculated by replacing the perturbation variable  with . The difference in hydration free energy between the cis ( = 0°) and trans ( = 180°) conformations, A_hyd(0°  180°), was determined at fixed values of  from 7.5° to 172.5°, in increments of 15.0° with  = ±7.5°.

The free energy of charging the ionic tetramers in water was calculated by scaling the solute charge, q_s, by . This free energy expanded to second order in q_s is (Hummer et al., 1996b; Pratt et al., 1994)

(4)

For a system that is truly infinite in size, the ion-water electrostatic energy is directly proportional to q_s multiplied by the electrostatic potential at the charge site,   E/q, and ²E/q_s² is zero. However, for the finite, periodic systems employed in our simulations, self-interactions arise (Hummer et al., 1996b) and must be accounted for in calculating the electrostatic energy. The self-interaction energy is quadratic in q_s, and thus ²E/q_s² is nonzero.

The only difference between positive and negative ions in Eq. 4 arises from E/q_s0, the electrostatic potential at zero charge. Previous studies have shown that this quantity is positive, thereby favoring the hydration of anions (Ashbaugh and Wood, 1997; Hummer et al., 1996b, 1997; Pratt et al., 1994). The fluctuations in the electrostatic potential, however, depend strongly on the sign of q_s (Garde et al., 1998; Hummer et al., 1996b), and in Eq. 4 these fluctuations are evaluated only at zero charge. Therefore, a more accurate expression for A_hyd(0  q_s), which incorporates this information at charge state q_s, is (Hummer and Szabo, 1996)

(5)

which is exact to fourth order in q_s. It should be noted that this expression requires simulation results only at the initial and final charge states to obtain an accurate estimate of A_hyd(0  q_s).

Canonical ensemble Monte Carlo (MC) simulations (Allen and Tildesley, 1987) of one solute molecule at infinite dilution in water were performed at a temperature of 25°C and a water density of 0.997 g/cm³. The details of these simulations are given in Table 1. Water was modeled using the simple point charge (SPC) potential (Berendsen et al., 1981). The n-alkanes were modeled using the OPLS united-atom LJ potential parameters (Jorgensen et al., 1984). Solute-water cross-interaction parameters were obtained from the geometric mean combining rules, i.e., _sw = (_ssww)^1/2 and _sw = (_ssww). Potential and structural parameters for water and the n-alkanes are given in Table 2. LJ parameters for the charged tetramers were assumed to be the same as those for uncharged n-butane. Minimum image LJ interactions were truncated on a site-by-site basis at half the simulation box length, L/2. Solute perturbations and averages were evaluated on the fly once every MC pass (see Table 1). Statistical uncertainties in the free energy were estimated from block averages obtained by dividing the simulation runs into five equal blocks that were assumed to be independent.

Electrostatic interactions were evaluated by the generalized reaction-field (GRF) method (Hummer et al., 1994). Although Ewald summation (De Leeuw et al., 1980) is generally regarded as the best available technique for calculating the electrostatic energy of a periodic system of charges, the method is slowed by the calculation of long-range Fourier space contributions. The GRF method, on the other hand, neglects long-range contributions to the energy and therefore can be evaluated more rapidly. Previous simulation studies have shown that the two methods give essentially the same pair correlation structure of aqueous electrolyte solutions (Hummer et al., 1994, 1996b). Most importantly, the two methods give identical charging free energies of both ionic (Hummer et al., 1996b) and polar solutes (Hummer et al., 1995).

The GRF electrostatic energy of a system of N_w water molecules and one solute molecule is

(6)

where q_w and q_s are the water and solute partial charges, and  is the number of solute charge sites. A factor of 1/4₀ is neglected in Eq. 6 for notational simplicity, where ₀ is the permittivity of free space. The effective GRF electrostatic pair potential depends only on the minimum image distance r between charges and has a cutoff distance r_c:

(7)

where (x) is the Heaviside unit-step function. The self-interaction potential is defined as (Hummer et al., 1995, 1996b) _GRF(r) = _GRF(r)  1/r. An electrostatic cutoff of r_c = L/2 is used throughout this work.

For a solute with one charge site ( = 1), the derivatives of the electrostatic energy with respect to q_s are

(8a)

and

(8b)

where

(8c)

_GRF is the direct GRF electrostatic potential at the solute charge site due to the solvent water molecules. _GRF(0) corrects for finite system size and potential cutoff effects on the electrostatic energy. For r_c = L/2, the GRF self-interaction term is _GRF(0) = 24/5L  /20L. In the limit of an infinite simulation size, it is easily confirmed that _GRF(0) = 0 and _GRF corresponds to the true 1/r electrostatic potential. Inclusion of the self-interaction energy makes the free energy of charging an ion robust such that the limiting infinite dilution value is approached for simulations of moderate size (N_w > 64) without further correction (Hummer et al., 1996b).

	CONTINUUM MODEL OF HYDRATION

The thermodynamic cycle employed by the continuum model to calculate A_hyd is depicted in Fig. 1. In the first step the solute charges are turned off in vacuum, A_vac(q_s  0). In the second step the uncharged solute is transferred from the vacuum to the condensed aqueous phase, A_H. Finally, the solute charges are turned back on in solution, A_aq(0  q_s). The work required to perform all three steps is the hydration free energy,

(9)

The individual contributions to A_hyd are evaluated as follows.

View larger version (33K): [in this window] [in a new window]   FIGURE 1   The thermodynamic path used to evaluate the hydration free energy of a charged solute in a fixed conformation.

The free energy of transferring n-alkanes (excluding methane) from vacuum to water is found empirically to be linear with the molecular surface area of the solute (Hermann, 1972; Chothia, 1974; Reynolds et al., 1974). This observation suggests the following expression for hydrophobic contributions to the free energy of hydration (step 2 in Fig. 1):

(10)

where  is solute surface area,  is the free energy/surface area coefficient, and b is the free energy intercept. The free energy/surface area coefficient is frequently referred to as a surface tension, although its connection to a macroscopic oil/water surface tension is merely anecdotal. Several definitions of solute surface area have been proposed: the van der Waals (vdW), solvent-accessible (SAS), and molecular surface areas. Differences between these definitions are discussed in detail elsewhere (Lee and Richards, 1971; Richards, 1977). All three surface areas were calculated using Connolly's Molecular Surface Program (Connolly, 1981, 1983) and a methyl/methylene group vdW radius of 1.9 Å and a water probe radius of 1.4 Å.

To evaluate the electrostatic contributions to the free energy (steps 1 and 3 of Fig. 1), the solute is treated as a low dielectric cavity imbedded within a high dielectric solvent (water) or in a vacuum. The resulting Poisson equation (Jackson, 1975) for the electrostatic potential is

(11a)

and

(11b)

where r is the spatial coordinate, _i and _e are the electrostatic potentials in the interior and exterior regions of the solute denoted by B_i and B_e respectively, _i is the dielectric constant of the solute interior,  is the number of interior solute charges, {r} is the set of charge positions, and (r) is the Dirac delta function. _i and _e are subject to the following boundary conditions at the solute/solvent interface:

(12a)

and

(12b)

where B denotes the interfacial surface bordering regions B_i and B_e, _e is the dielectric constant of the solvent, and /n is the derivative with respect to the outward unit surface normal of B.

Poisson's equation must be solved numerically for polyatomic solutes. This is accomplished by using the boundary element method (BEM) to solve the surface integral form of Poisson's equation for _i and _i/n on B (Horvath et al., 1996; Pratt et al., 1997; Rashin, 1990; Yoon and Lenhoff, 1990; Zauhar and Morgan, 1985). The electrostatic potential at a position r in the solute interior is then obtained from

(13a)

where

(13b)

and d(r') is a differential area element of the solute surface at r'  B. The first term on the right-hand side of Eq. 13a is due to the direct electrostatic interactions with the interior solute charges, and _surf(r) is the solvent contribution to the potential arising at the solute surface. The electrostatic free energy of transferring the solute from vacuum to aqueous solution is

(14)

where A_elec = A_vac(q_s  0) + A_aq(0  q_s), _aq is the dielectric constant of water, and _e = 1 is the dielectric constant of a vacuum. _i is assumed to be equal to 1 because the simulations neglect molecular polarizability. In this case there is no dielectric discontinuity at the solute boundary in a vacuum and A_vac(q_s  0) is zero. _aq is set equal to 65, the dielectric constant of SPC water at ambient conditions (Hummer et al., 1995, 1996b). For large _e the continuum free energy is insensitive to the value of _e.

For a spherical ion in solution, Poisson's equation can be solved analytically. In this case, the charging free energy is given by the Born equation (Born, 1920):

(15)

where R is the Born radius of the ion. The Born radius is typically treated as an adjustable parameter that reflects the effects of water structure in the vicinity of the ion. We model the molecular solutes as a collection of spheres with Born radii that depend on the final charge of each carbon site. The radius of an uncharged methyl/methylene group, R₀, is taken to be the vdW radius, 1.9 Å. The Born radii of positively (R₊) and negatively (R) charged methyl groups are fit to free energies obtained from explicit water simulations. The solute interior, B_i, is defined by the molecular surface (Richards, 1977), which is calculated by rolling a 1.4-Å-radius probe over the Born radii of the solute in a given conformation.

The surface discretization of B was performed using Zauhar's SMART (Smooth MoleculAR Triangulator) program (Zauhar, 1995). Rather than approximating each triangluar surface element as flat, this algorithm generates a more accurate curvilinear representation of the solute surface. An angle parameter of 25° was used to triangulate the solute molecular surface (Zauhar, 1995). Using a finer element size did not affect the results significantly. The SMART discretized surface was input into the BEM Poisson solver developed by Neal and Lenhoff (Neal, 1997). _i and _i/n were assumed to vary linearly over each element, and the surface integrals were evaluated using Gaussian quadrature.

	RESULTS AND DISCUSSION

Top Abstract Introduction Results & Discussion Conclusions References

n-Alkanes in aqueous solution

The calculated n-alkane free energies of hydration in SPC water are reported in Table 3. Previously reported simulation results for TIP4P water (Kaminski et al., 1994) and experimental values from the literature (Ben-Naim and Marcus, 1984; McAuliffe, 1966) are also included in this table. Not surprisingly, A_hyd is positive for all of the n-alkanes, indicative of their low solubilities in water. The free energies obtained from the SPC water simulations are in very good agreement with those calculated from the TIP4P water simulations. The simulations also predict, in accord with experiment, that A_hyd decreases from methane to ethane, but then increases with increasing carbon number for longer n-alkanes. However, the simulations significantly overestimate the experimental values of A_hyd. The disparities between simulations and experiment are due in large part to the discrepancy in the incremental free energy of adding a methylene group to ethane: A_hyd(C₂  C₃) is 0.21 kcal/mol from experiment and 0.73 ± 0.04 kcal/mol from the SPC water simulations. Differences between simulation and experiment also occur to some extent as a consequence of systematic differences in the OPLS united-atom parameters for methane, ethane, and propane (Table 2). For propane and the longer n-alkanes, the united-atom parameters are independent of chain length. Consequently, the incremental free energy of hydration of adding a methylene group to these alkane chains obtained from simulation is approximately constant (~0.2 kcal/mol) and is in good agreement with experiment. Plotting the cumulative values of A_hyd(C_i1  C_i) as a function of  clearly illustrates this behavior (Fig. 2). The incremental free energy profiles for methane, ethane, and propane (i = 1, 2, 3) show little resemblence to one another. However, for n-butane, n-pentane, and n-hexane (i = 4, 5, 6), these profiles are practically indistinguishable.

View larger version (23K): [in this window] [in a new window]   FIGURE 2   Cumulative values of Ahyd(Ci1  Ci) as a function of . Symbols correspond to i = 1 (), i = 2 (), i = 3 (), i = 4 (), i = 5 (), and i = 6 (). Results for i = 1 are plotted in the inset.

Free energies of hydration are plotted as a function of n-alkane SAS area in Fig. 3. Excluding methane, the experimental values show a linear dependence on SAS area. A linear regression of the data gives _SAS^expt = 6.8 cal/(mol Ų) and b_SAS^expt = 0.6 kcal/mol. The values obtained from simulation do not show the same linear dependence below the SAS area for propane, which can be attributed to the systematic differences in the OPLS united-atom parameters discussed above. For the simulation results, _SAS^sim = 7.4 cal/(mol Ų) and b_SAS^sim = 1.8 kcal/mol, which is obtained by fitting SAS areas for propane and the longer chain alkanes. Free energies of hydrophobic hydration are typically correlated with solute SAS area (Chothia, 1974; Hermann, 1972; Reynolds et al., 1974; Schmidt and Fine, 1994; Sitkoff et al., 1994, 1996). The vdW surface area or the molecular surface area could be used as well. Indeed, all three surface areas correlate linearly with A_hyd equally well, which is not surprising, because all three surfaces are linearly correlated with each other for the n-alkanes. The calculated values of and b for the three surface areas are given in Table 4. In each case, the value of  obtained from simulation is somewhat greater than the experimental value. Apart from the approximate nature of the intermolecular interactions used in the simulations, a major difference between simulation and experiment is the constraint of fixed all-trans conformation for the n-alkanes imposed on the simulations. The experimental results, on the other hand, are naturally averaged over all available conformations. The value of  will not change appreciably, however, if the simulation constraint is relaxed to sample more compact conformations with lower surface areas, because the maximum change in solute surface area for the allowed conformational changes is small compared to the total area. For example, the SAS area of n-butane in the cis and trans conformations differs by only 4%.

View larger version (15K): [in this window] [in a new window]   FIGURE 3   Ahyd as a function of n-alkane SAS area. , The SPC water simulation results from this work. , TIP4P water simulation results (Kaminski et al., 1994). , Experimental values (Ben-Naim and Marcus, 1984; McAuliffe, 1966). The solid lines are the best fits of Eq. 10 to the results for the longer chain alkanes. The line through the experimental data is Ahyd = 6.8 cal/(mol Å2) + 0.6 kcal/mol, and the line through the simulation data is Ahyd = 7.4 cal/(mol Å2) + 1.8 kcal/mol. , The experimental value for Ahyd of cyclohexane.

View this table: [in this window] [in a new window]   TABLE 4   Values of  and b (Eq. 10) for the n-alkanes in water calculated from SPC water simulations and from experimental values for n-alkane solubilities in water

Conformational free energy of n-butane

The simulation results for A_hyd() as a function of n-butane conformation are shown in Fig. 4. The increase in A_hyd() with increasing  indicates that hydration destabilizes the elongated trans conformation of n-butane relative to the more compact cis conformation. This behavior is consistent with the notion that hydrophobic interactions tend to minimize the extent of oil-water contact. The A_hyd() profiles are also in qualitative agreement with previous simulation (Jorgensen, 1982; Jorgensen and Buckner, 1987; Kaminski et al., 1994; Rosenberg et al., 1982; Tobias and Brooks, 1990) and theoretical (Hummer et al., 1996a; Pratt and Chandler, 1977; Zichi and Rossky, 1986) studies of the free energy of hydration of different n-butane conformations in water. In Fig. 5, A_hyd() is plotted as a function of n-butane surface area rather than the backbone dihedral angle. All three surface area correlations show qualitatively a linear dependence of A_hyd on surface area, as suggested by Eq. 10. We calculate  = 109, 49.3, and 111 cal/(mol Ų) for the vdW, solvent-accessible, and molecular surfaces, respectively. The vdW surface area correlation gives the best fit of the simulation results based on root mean square difference (see Fig. 5). The vdW surface area correlation also captures most accurately the plateau in A_hyd as a function of  for  > 90° (Fig. 4). The other two surface area correlations give plateaus closer to = 180° (trans conformation). See, for example, the fit of Eq. 10 to the simulation results, using the SAS area in Fig. 4. The success of the vdW surface area correlation is surprising in light of recent studies indicating that the molecular surface area best describes the hydrophobic interactions between methane pairs in water (Jackson and Sternberg, 1993, 1994; Pitarch et al., 1996; Rank and Baker, 1997). Interestingly, the suggested value of for methane-methane pair interactions, calculated using the molecular surface area, lies between 104 and 110 cal/(mol Ų), in very good agreement with the  value for the molecular surface area correlation in Fig. 5.

View larger version (18K): [in this window] [in a new window]   FIGURE 4   Ahyd as a function of n-butane conformation referenced to Ahyd(0°) = 0. , SPC water simulation results. , The fit of Eq. 10 to the simulation results, using the vdW surface ( = 109 cal/(mol Å2)). , The fit of Eq. 10, using the SAS area ( = 49.3 cal/(mol Å2)). Error bars on the simulation results denote one standard deviation.

View larger version (19K): [in this window] [in a new window]   FIGURE 5   Ahyd as a function of n-butane conformation as correlated with solute surface area referenced to (0°)  0 (cis conformation). , vdW surface area ( = 109 cal/(mol Å2)); , molecular surface area ( = 111 cal/(mol Å2)); , SAS area ( = 49.3 cal/(mol Å2)). The plots for molecular surface and vdW surface areas are shifted down by 0.25 and 0.5 kcal/mol, respectively, from the plot for SAS area. The root mean square difference between the simulation results and the linear fit of Eq. 10 are 0.01, 0.03, and 0.03 kcal/mol for the vdW, solvent-accessible, and molecular surfaces, respectively.

The values of  extracted from the n-alkane free energies of hydration (Table 4) are an order of magnitude less than the values extracted from fits of A_hyd as a function of n-butane surface area (Fig. 5), regardless of how solute surface area is defined. The justification for surface area scaling of the free energy embodied in Eq. 10 is based on the expectation that the free energy of forming a macroscopic surface scales with its surface area. It is not evident, however, that this scaling should also hold for molecular length scales, such as those associated with the surface areas of simple linear alkanes. Previous simulation studies of n-butane in water have shown that the entropy of hydrophobic hydration per water molecule depends on solute conformation (Ashbaugh and Paulaitis, 1996). Simulation studies of the hydration of idealized spherical and ellipsoidal solutes have likewise shown that differences between surface area derivatives of the hydration free energy could be reconciled only if solute curvature was taken into account (Wallqvist and Berne, 1995b). These results suggest that surface area scaling breaks down for molecular length scales.

Experimental evidence for this breakdown can be found by comparing values of A_hyd derived from the measured solubilities of cycloalkanes in water to those derived from the solubilities of their n-alkane counterparts. The cycloalkanes can be thought of in this comparison as "collapsed" linear analogs of the n-alkanes constrained to ring-forming conformations. A_hyd = 0.80, 1.20, and 1.20 kcal/mol at 25°C for cyclopropane, cyclopentane, and cyclohexane (Horvath et al., 1996), respectively. These values are much less than A_hyd for ethane (Table 3), even though ethane is the most soluble n-alkane and has a lower surface area than the cycloalkanes. Moreover, as shown in Fig. 3, A_hyd for cyclohexane falls well below the linear correlation for experimental hydration free energies of the n-alkanes. Based on the experimental A_hyd and calculated SAS area differences betweeen n-hexane and cyclohexane,  is 34 cal/(mol Ų). This free energy/surface area coefficient is much closer in magnitude to the value we would calculate from the conformational equilibria. Similar conclusions can be reached for the other cycloalkanes.

We propose the following generalization of Eq. 10 to explain the large differences in free energy/surface area coefficients derived from hydration free energies and from conformational equilibria. The reversible work required to form a hydrophobic surface in water is equal to A, where  is the free energy/surface area coefficient calculated from n-alkane free energies of hydration as a function of solute surface area (Table 4). The change in the free energy of hydration associated with a conformational change for n-butane is given by

(16)

The first term on the right-hand side of this equation is related to the change in hydration free energy due to a change in solute surface area. The second term reflects the change in hydration as a function of solute conformation. The free energy/surface area coefficient evaluated from conformational equilibria is simply

(17)

Thus the two free energy/surface area coefficients are equivalent when / = 0. However, when /  0,  and  can be substantially different, even if / is small, because the second term in Eq. 17 is this derivative multiplied by the total solute area, a comparatively large quantity. Based on the vdW surface values ( = 12.0 cal/(mol Ų),  = 109 cal/(mol Ų), and  = 98.1 Å²), we estimate that / = 1.0 cal/(mol Å⁴). Thus  varies from 8.3 cal/(mol Ų) for cis n-butane to 12.0 cal/(mol Ų) for trans n-butane.

Free energy of charging trans-butane

The simulation averages for the electrostatic potential and fluctuations in the electrostatic potential required to calculate the free energy of charging the terminal methyl group of trans-butane are reported in Table 5. From Eq. 5, A_hyd(0e  +1e) = 56.8 ± 0.8 kcal/mol and A_hyd(0e  1e) = 94.8 ± 1.1 kcal/mol. Based on averages for the electrostatic potential and fluctuations in the electrostatic potential recently reported for simulations of charging methane in SPC water (Hummer et al., 1996b), we calculate A_hyd(0e  +1e) = 60.3 kcal/mol and A_hyd(0e  1e) = 106.6 kcal/mol, using Eq. 5. The free energies of charging methane are slightly more negative, as expected, because the major contribution to the free energy is the electrostatic interaction with water molecules in the first hydration shell of the ion. Methane has more water molecules in its first hydration shell than does the terminal methyl group of n-butane, because of the excluded volume of the attached propyl group for the latter. In both cases, the free energy of charging the anion is considerably more negative than that for the cation.

The difference in anion and cation hydration can be assigned in part to the nonzero electrostatic potential at zero charge: E/q_s0 = 8.4 ± 0.4 kcal/(mol e) (Table 5). Subtracting q_sE/q_s0 from A_hyd effectively removes this contribution from the total free energy of charging, which reduces the asymmetry in ion hydration: q_sE/q_s0 accounts for 10-15% of the total free energy of charging these ions.¹ Nonetheless, the difference is still significant: 65.2 ± 0.8 kcal/mol and 86.4 ± 1.1 kcal/mol for charging the cation and the anion, respectively. This remaining difference is due in large part to differences in the microscopic structure of water around oppositely charged ions, which is a manifestation of the asymmetrical spatial distribution of charges in the water molecule. That is, anion hydration is favored over cation hydration because the positively charged hydrogens of water molecules in the first hydration shell can reside closer to the anion than can the negatively charged oxygens of water molecules in the first hydration shell of the cation (Garde et al., 1998; Hummer et al., 1996b).

The Born radii for the charged methyl groups obtained by fitting the unadjusted free energies (i.e., those including a nonzero value of E/q_s0) are R₊ = 2.82 Å and R = 1.56 Å, which are in reasonable agreement with the Born radii obtained for charged methane: R₊ = 2.71 Å and R = 1.54 Å. The larger R₊ relative to R is consistent with the asymmetry in ion hydration attributed above to water structure. A consequence of assuming that the solvent responds as a macroscopic dielectric is that the free energy of charging is proportional to q_s² and therefore is independent of the sign of q_s (e.g., Eq. 15). Hence E/q_s0 from the continuum model mus