J. Phys. Chem. 1995, 99, 13118- 13125
13118
Free-Energy Cost of Bending n-Dodecane in Aqueous Solution. Influence of the Hydrophobic Effect and Solvent Exposed Area A. Wallqvist" and D. G . Cove11 Frederick Cancer Research and Development Center, National Cancer Institute, Science Applications International Corporation, Frederick, Maryland 21 702 Received: March 30, 1995@
Computer simulations of the free energy associated with different conformations of a single n-dodecane molecule in liquid water are reported. The alkane chain was monitored as a function of the end-to-end distance, essentially following the formation of a hairpin bend from an initially fully extended state. While elongated conformations are among the most stable state in both gas and solvated states, the influence of the solvent is to favor more compact states. At the most compact conformation, the hairpin bend, the solvent imparts a stabilizing free-energy contribution of 6 W/mol. The torsional strain of bending the molecule into this hairpin conformation is partly relieved by the reduction of entropically unfavorable water molecules associated with the reduction in the alkane's surface exposed area. The resultant solvent contribution to the free-energy change is characterized by an unfavorable enthalpic and a competing favorable entropic (-TAS) contribution. The conformational change associated with the transition from an extended state to a compact state for n-dodecane is accompanied by a relatively small, non-uniform change in surface exposed area. A simplified model of the free energy as a function of this area can account for the gross effects but is unable pick up the molecular details inherent in the system.
1. Introduction
The solvation and aggregation of nonpolar molecules in aqueous media are different manifestations of the hydrophobic effect. As apolar molecules are relatively insoluble in water they tend to avoid aqueous environments, and if the concentration is sufficient, they create their own hydrophobic phase by forming larger aggregates. The ubiquitous presence of these hydrophobic forces in virtually all biochemical systems has been appreciated for a long time,' though their magnitude and importance have been debated.2 Hydrophobicity is thus thought of as the driving force to collapse initially extended protein assemblies into a nonnative molten g l o b ~ l e . ~Hydrophilic attractions, on the other hand, are often thought to be of a more specific nature, e.g. the specific matching between hydrogen bonding groups that occurs among DNA base pairs and within the interior of globular proteins. Yet the role of surface hydrophobicities in providing identification sites for molecular recognition among protein-protein complexes is ~ t r i k i n gThe .~ molecular nature of the hydrophobic effect is intimately connected with both the solute and solvent interactions in the presence of each other. The large gain in free energy, upon removal of portions of solute molecules that do not interact particularly strongly with water molecules, is due partially to the more favorable water-water interactions that can develop as a consequence. It is, however, not a purely enthalpic effect; strong entropic contributions are also postulated to arise from the release of water molecules from a hydrophobic surface. Computer simulation of molecularly detailed systems has proven to be a particular useful means of probing the nature of hydrophobic interactions in terms of explicit solvent and solute molecules. Thus initial efforts at showing the hydrophobic effect between rare-gas or methane molecule^,^-^^ revealed the importance of including explicit solvent molecules in order to obtain a quantitative description of the system's free energy. The resultant water-cage structures around the solutes 'Abstract published in Advance ACS Abstrucrs, July 1, 1995.
0022-365419512099-13118$09.0010
were found to actually stabilize the solvent-separated solute pairs, a counter-intuitive result impossible to predict from empirical models of surface exposure. Although these watercage structures are reminiscent of clathrate structure^,'^ they are not rigid; they exchange solvent molecules with the bulk and possess properties that are more similar to the liquid phase than to the ice The conformational equilibrium of a chain molecule is highly dependent on its surrounding environment. Consequently the population of compact versus extended states for a single chain molecule will depend on whether the molecule is in the gas phase, a liquid hydrocarbon phase, or an aqueous phase. The central dogma of the hydrophobic effect is that hydrophobic entities are being induced to aggregate in an aqueous phase, primarily through the influence of an entropic reduction of the surrounding water molecules. The same forces should also be operative for a chain molecule and should influence the conformational equilibrium. Theoretical investigations of nbutane in various phases have shown inconclusive result^;'^-^^ the energetic consequences of changes toward a compact gauche state have been found to be negligible or only slightly energetically favorable. The entropic contribution to the hydrophobic interaction for this system have also been disp ~ t e d . ~ ~The . ' ~variability of these results stems largely from the difference in simulation and interactions parameters used for the systems investigated. The current report does indeed show that the hydrophobic effect is present and of great importance in the compaction of hydrocarbon chains. Studying a comparatively large system, a 12-monomer alkane chain versus n-butane, permits additional assessment of the hydrophobic interaction. In our system the hydrophobic effect is not large enough to actually shift the equilibrium states of extended chains to a compact chain, but the contributing solvent free-energy component is unequivocally present and is compatible with a large negative entropic (-TAS) and a competing positive enthalpic component. The solvent contribution per water molecule surrounding the alkane molecule 0 1995 American Chemical Society
Free Energy of n-Dodecane in Water is only -0.1 kJ/mol-a value that is hard to accurately evaluate for smaller systems where the number of affected water molecules are fewer, as is the case for n-butane. The solvent contribution to the free energy of bending the n-dodecane molecule has been explored using measures of ~ ~ .find ~ ~ that the concept surface area of the chain m o l e c ~ l e . We of parametrizing free-energy changes as proportional to a surface area is only of a qualitative nature. The nonquantitative nature of this relationship stems from cooperative effects among the waters associated with the alkane chain. Effectively the freeenergy penalty of maintaining one water molecule in the water shell of an extended chain is different from that of a compact chain. This is not inconsistent with the empirical observation that transfer free energies of n-alkanes between water and the gas phase are proportional to the surface area of the extended alkane chain. The equilibrium conformations of low molecular weight alkane molecules correspond to an essentially extended state, and since the free-energy penalty per water molecule for maintaining shell waters around an extended state is constant, the free energy of transfer is proportional to the number of methylene units in the chain. As we show below, the freeenergy penalty is however dependent on the conformation of the chain molecule. Thus, even though the concept of area is useful for some cases, it is not a general result and should not be used indiscriminately between different systems. The results presented here also show that the bending motion of the 12-monomer model hydrocarbon chain is composed of three well-defined parts. First, there is a large-scale overall bending motion that essentially leaves the hydration shell intact. This is commensurate with an extended state and is characteristic of the set of equilibrium conformations present in the aqueous solution for our model system. Further contraction of the chain is hindered by an energetic barrier which must be overcome in order to proceed to a compact state. For n-dodecane this barrier is approximately 7 kJ/mol and is a result of the internal motion of the alkane chain. Once a semicompact state has been achieved, this conformation is stabilized by the solvent. Maintaining a water molecule at an alkane/water interface is associated with a free-energy penalty; water molecules would rather be solvated in a bulk liquid phase. The removal of such energetically unfavorable water molecules via a lowered surface area in the compact state, as compared to the extended state, is the driving force behind this stabilization.
2. Computational Procedures
J. Phys. Chem., Vol. 99, No. 35, 1995 13119 verify that we can reproduce the distribution, and hence the free energy, of the gas-phase system by calculating the free energy via a perturbation path along this “reaction” coordinate. Then we repeat the free-energy calculation in the presence of the solvent. The equations of motions were integrated using the Rattle version2’ of the velocity Verlet algorithm2*in order to maintain the internal bond lengths as well as the bond angle of the water molecule. Similarly, the constrained end-to-end distance of the hydrocarbon chain was also maintained via this algorithm. The time step was set to 2.0 fs. Temperatures were maintained at room temperature by periodically rescaling the velocities. Both the translational and rotational temperatures were monitored so as to avoid a temperature imbalance between these degrees of freedom. All interactions were spherically truncated at 9.3 A, which has been shown to be adequate to account for the longrange interactions in most water proper tie^.^^ 2.1. Model Interactions. For water-water interactions we have employed a recently developed effective, painvise additive potential RER,30 whose molecule-molecule potential energy is given by
where {R,}denotes the coordinates of water molecule i, roo is the distance between oxygen atoms on water i and j, and rG is the radial distance between two atoms on the ith and jth water molecule. This effective pair-potential model has a permanent dipole moment of 2.60 D. The total potential energy for this nonpolarizable model needs to be adjusted for the energy cost of creating the effective dipole moments, Vself.3i-30 The parameter values for the model are given in Table la. The hydrocarbon chain has intramolecular energy contributions arising from bond-angle bending and torsional motion. Bond vibrations are suppressed by keeping the bond length fixed to 1.53 A. Bond angles were introduce to properly account for the energy flow between bond-bending and the dihedral angle rotations. A harmonic bond-bending term between three connected methylene units forming an angle 8 is thus used,
In this investigation we have used 717 water molecules enclosed in a periodically replicated cubic box whose volume (3) is 21 718 A3, together with a single n-dodecane molecule. Although the water model explicitly includes both oxygen and where ke is the harmonic force constant and 8 0 is the equilibrium hydrogen atoms, the hydrocarbon chain was modeled as a set bond angle. Parameter values are those of Wiener et al.33and of connected CH3- or -CH2- units. The exclusion of hyare given in Table lb. The dihedral potential was taken from drogen atoms and the subsequent reparametrization of the the work of J o r g e n ~ e nand ~ ~is represented as a truncated Fourier potential energies can account for variations of up to 10% in the equilibrium between gauche and trans states in n - b ~ t a n e . ~ ~ series. Simulations were carried out both in the gas phase and in the V&ri,rj,rk,r,) = V, Vi cos 4 V2 cos 24 V3cos 34 (4) aqueous phase in order to determine the free energy as a function of conformation. While a n-dodecane molecule has 9 dihedral angles, a complete characterization of the free-energy surface where C#I is the dihedral angle between the connected methylene of all possible conformations is untractable. Instead we have units, r;, rj, rk, and rl. The coefficients are summarized in Table used the end-to-end distance of the chain molecule, 1b. Intramolecular interactions between methylene units that are separated by four or more bonds utilizes a Lennard-Jones potential with the parameters given in Table 1b.34.35 The water oxygen atoms interacts with the methyl and methylene units on the n-dodecane molecule via a Lennardwhere f1.N denotes the coordinates of the end groups on the Jones interaction between atom pairs i and j as dodecane molecule, as our “reaction” coordinate. We initially
+
+
+
Wallqvist and Cove11
13120 J. Phys. Chem., Vol. 99, No. 35, 1995 I7
TABLE 1
Gas-Phase Perturbation +
(a) Water Interaction Parameters" 40, e
qn, e CIZ,kJl(mol&'?) c6. W/(mol Ab)
-0.920 0.460 3 500 000 -3100
C4,kJ/(mol C,, Wlmol Wr,
A?)
15.0 - 1.000 1.5 4.5
,A-2
rr, A
8
(b) Hydrocarbon Interaction Parameters Dihedral Potential
VO
9.74 kJ/mol 2.95 kJ/mol
VI
-0.57 kJ/mol 6.58 kJ/mol
V2 v3
Bending Potential ks
527.2 kJ/mol/rad?
-A
00
1.962 rad
Intramolecular LJ Potential A,,'?, 107
C,,'?, 103
atom i/atom j
(A1?kJ)/mol
(A6kJ)/mol
CH,/CH' CH?/CH3 CH?/CH?
3.682 2.483 3.023
10.385 7.004 8.527
(c) WaterkIydrocarbon Interaction Parameters: LJ Potential atom i/atom j
,A,,~z,107 (A'?kJ)/mol
O(HzO)/CH, O(H*O)/CH?
0.96 1 0.789
C,'?, IO' kJ)/moI
(A6
5.146 4.226
(7 The monomer geometry of a water molecule is given by a bond length ro = 0.96 8, and a bond angle 0, = 104.52'. In order to recover the potential energy of the water-water interaction in kilojoules per mole when using eq 2, I / ( ~ X E ~ )should be set to 1389.0. The correction term to be added to the RER potential arising from the self-energy of the dipole moment, V5e~f, is 13.2 kJ/mol at 300 K.
where the ry denotes the interatomic distance and the A'* and c6 coefficients are given in Table lc.34 2.2. Free-Energy Calculations. Thermodynamic and statistical perturbation theory originally due to Z ~ a n z i g ~ was ~.~' used to calculate relative free energies for different end-to-end separations, ree, of the hydrocarbon chain. The change in Helmholtz free energy between two systems ro and rl containing N particles at a given temperature, T, and volume, V, is given by
AA
= A ( r , ) - A(ro) = -kTln(e-[Li(rl)-u(rn)llkT >o ( 6 )
where k is the Boltzmann constant and the brackets indicate an ensemble average of the exponential of perturbation in the reference system. The entropy contribution to the free-energy change in the system can be evaluated from the temperature derivative of the Helmholtz free energy to ~ i e l d ' ~ . ~ ~
where 5G and SU;; are the perturbed and the reference Hamiltonians, respectively. Even though this estimator is associated with a large uncertainty as calculated from molecular simulations, it is still useful for gauging the overall entropic response of the system. A set of simulations with the end-to-end separation, re,, between 14.0 and 4.0 8,was carried out in order to characterize
Figure 1. Helmholtz free energy of n-dodecane in the gas phase as a function of the end-to-end distance, re,,calculated via the perturbation approach (+). The potential of mean force calculated from the actual gas-phase distribution of re,-distances, P(ree),from an unconstrained simulation in the gas phase (0).The overall change in free energy is calculated with respect to the fully extended, all-trans chain.
the free-energy change according to eq 6 and eq 7. A separation of 14.0 8,characterizes the fully extended all-trans conformation and serves as the overall reference state, while 4.0 8, is the turning point of the LJ interaction between two end-group units. The distance, re,, was changed in steps of 0.1 8, for the reference states, and the free-energy changes were accumulated in intervals of f0.02 8, for each reference state. Thermodynamic averages were then calculated from block averaging five simulations of 5.0 ps each. The error of the free energy for the perturbed systems was estimated as twice the standard deviation obtained from block averaging. As the perturbation technique estimates the change in free energy relative to each reference system simulated, the total free energy change is reconstructed by connecting the individual free energy curves from each reference system. The total simulation time was approximately 1000 CPU hours using a DEC-alpha Model 3000/800 workstation.
3. Results and Discussions 3.1. n-Dodecane in Gas Phase. In order to characterize a reference state and investigate the feasibility of using the endto-end distance as a perturbation path to describe the conformations of the dodecane molecule, we performed both unconstrained simulations and perturbation calculations of n-dodecane in the gas phase. We calculated the probability distribution of end-to-end distances, P(ree),for the free, unconstrained molecule and converted this probability to a potential of mean force, W(ree),or the Helmholtz free energy, via39
+
"(re,) = kT In P(ree) C
(8)
where k is the Boltzmann constant, Tis the absolute temperature, and C is an arbitrary constant. In Figure 1 we can gauge the correspondence between the calculated free energies using the actual distribution, AW(ree),and the free-energy change estimated from the perturbation path for slowly decreasing the constrained end-to-end distance, AAgas(ree).The largest deviation between the perturbation path and the gas-phase simulations lies in the transition from the fully extended state at re, = 14.0 A, assigned a free energy of 0.0, and states compatible with 10 < ree/8, < 13. The artificial fixed coupling between the end units introduces a correlated motion along the chain that is not
Free Energy of n-Dodecane in Water
J. Phys. Chem., Vol. 99, No. 35, 1995 13121
6
5
T4 z
4
3
2 3
6
9
12
15
rw/A
Figure 3. Snapshots of the n-dodecane conformation for different constraiaed ree’s. From right to left ree= 14.0, 12.0, 10.0, 8.0,6.0, and 4.0 A. The initial compaction of the chain is via a dihedral rotation at the chain end, whereas below 10 %, a bend roughly located at the midpoint of the chain is developed.
Figure 2. Distribution of rg and reevalues in the gas phase indicated by the contour diagram. The perturbation path in the gas phase is indicated by ( 0 )and closely follows the most probable distribution from the unconstrained simujation. This path is characterized by the relation r, = 0.17ree 2.0 A. The liquid-phase perturbation path is marked by asterisks and similarly follows the most probable path.
inn, Aqueous-Phase (aq) u Gas-Phase (gas) * Influence of Solvent (sol) +e-
+
present for the free chain. This appears to create an effectively larger barrier in the first rotation of a dihedral angle at the chain end. Yet the weights of the states at 10 -= reeiA < 13 are sufficient to ensure that both of these regions are the most probable states in the two systems-accordingly this is what we refer to as an extended chain. The subsequent collapse to a compact state is well reproduced by the free-energy perturbation scheme. This allows us to estimate that the fully extended state is about 8 kJ/mol more stable than the compact state at ree = 4.0 A. Another way of characterizing the conformational equilibrium of the chain molecule is to define a pseudo “equilibrium con~tant”,~ K, between configurations that are compact (ree < 4 A) or extended (ree > 9 A) as -W(r)lkT 2
r dr
-W(r)lkT 2
50
c
E. 00
. . 4 r,
Irl
a
-5.0
L$
-10.0 4
6
8
rc,JA
10
12
14
Figure 4. Calculated free-energy change in the gas phase and in the aqueous phase as a function of ree. Solvent contributions to the freeenergy change in the aqueous environment AA,,I,, is defined as the difference between the gas- and aqueous free-energy curves. The extended state at ree = 14 A is arbitrarily assigned a free energy of zero.
(9)
r dr
where we have somewhat arbitrarily chosen the definition of compact versus extended conformations. The yield of Kga,as 19 indicates that the number of observed compact state in the gas phase is indeed small. 3.2. n-Dodecane in Aqueous Phase. Although the freeenergy change could be reproduced in the gas phase by the perturbation method, this does not a priori need to he true in the aqueous phase. As a further check on the validity of the procedure, we monitored the radius of gyration as a function of end-to-end distance for the various conformations generated. Thus, a radius of gyration, rg, was calculated from the massweighted average distance from the center of mass, rcom, as
where mi is the mass of the united atom and M is the total mass of the molecule. The distribution of rg and re, in the gas phase is given as a contour diagram in Figure 2. There is a roughly linear relation between the end-to-end distance and the compactness of the conformations, as indicated by smaller r,-values at smaller ree. The perturbation path in the gas phase is identical
to the maxima of this distribution; similarly the free-energy calculation in the aqueous phase closely resembles this most probable path. On the basis of this agreement, we believe that this approach represents a reasonable means of calculating the free-energy of the conformational collapse along the most probable path, even for the solvated case. This result is specific to the system; as the number of monomer units grows in the alkane chain, the end-to-end distance is no longer expected to be able to adequately characterize the full set of conformational variations of the molecule. A set of configurational snapshots taken from the simulation of the solvated n-dodecane molecules in Figure 3 shows the process of bending the chain. From an initially fully extended molecule the perturbation path takes the chain into its most compact form, a hairpin bend characterized by re, = 4.0 8, and (r,) = 2.7 A. Water molecules are expelled from the interstitial volume between the two strands of the bent conformation around re, = 7 A. 3.2.1. Free Energy. In Figure 4 we compare the free-energy change upon contracting the hydrocarbon chain ends from 14 t o 4 A in both gas and aqueous phases. The solvent contribution is indicated as a difference between the calculated free-energy changes,
13122 J. Phys. Chem., Vol. 99, No. 35, 1995
-
Wallqvist and Cove11 Free Energy Entropy Enthalpy
-
1SO
*
Extended
Compact
+
%
1.25
10
8. 24 ri
s
g o
-
a 0
A
42
1 .oo
0
s
0.75
51
+
*)
0.50
0.25 0
*”
4
6
IO
8
12
14
T d A
n.nn __
0
2
6
4
8
10
rcoJ8,
Figure 5. Solvent contribution AAsO, subdivided into AH,,, and -TASso1 components from the differences between gas and aqueous phases. The competing nature of opposing enthalpic and entropic
contributions is evident. We can see that the general features of the extended states resemble each other regardless of the n-dodecane environment. Even though this qualitative similarity in behavior exists, the influence of the solvent is to decrease the stability of the extended state. Thus, the values for AAsOl,in the range 9 < re,lA < 14 are positive. As the molecule bends further, as in Figure 3, we closely follow the gas-phase free-energy curve. In the range 7 < ree/A 9 there is an almost vanishingly small contribution from the solvent. Further below this end-to-end distance any interstitial waters between the hydrocarbon chain ends are expelled and a significant stabilizing force arising from solvent contributions is developed. In this range, AAaq becomes increasingly less repulsive for the most compact state at ree 4 A. The “equilibrium constant”, Kaq,defined in eq 9, reduces to 14, indicating that the solvent is favoring the compact state, but for this system the hydrophobic effect is not large enough to drive the system into a compact equilibrium conformation. Given the large solvent contribution to AA favoring compact conformations, it can be speculated that AAsol. is indeed the principal driving force of the collapse of extended hydrophobic aggregates in liquid water. At a sufficiently close end-to-end distance the two terminal monomer units will experience a repulsion, as their electron clouds will be close enough to experience an exchange repulsion force, and the free energy of the system will be dominated by this repulsive interaction. A further subdivision of the solvent Contribution into enthalpic AHsOl.and entropic -TASsol, components can be done by comparing the estimates of the entropy via eq 7 for the two phases considered. In Figure 5 these quantities are given as a function of the end-to-end distance. Even though the error estimates are larger than for the free energy, two significant points can be made. First of all the individual contributions to AAsol, are larger in magnitude than AAS01.itself, and they are of opposite signs. Second, for the extended states (ree > 9 A) the enthalpic contribution is negative while that of -TASs01, is positive; both are relatively constant and essentially compensate each other. At lower end-to-end distances, where the conformations are more compact, the enthalpic losses are compensated by entropic gains, yielding a net attractive contribution to the free energy. 3.2.2. Molecular Properties of the Solvent. As the solvent envelope around the dodecane molecule consists of individual water entities, it is of interest to determine any changes in properties of these molecules as the chain bends. As water
-
Figure 6. Radial distribution function gco(r) for reC= 13 A, denoted as an “extended”conformation,and a corresponding “compact” gco(r) for ret 4 A. Both distributions display a clearly exclud$d region between 0.0 and 2.9 A, a shell region in the range 3.0-5.5 A, beyond which there is a bulk water phase.
molecules, in contrast to e.g. CC4, have a strong local anisotropy, small changes in orientational or positional properties may influence energetic properties. Even though such changes in properties are expected to be small, the sheer number of water molecules affected by a macromolecule can determine the energetic balance of the system. Although n-dodecane is not a macromolecule, it is associated with a substantially larger number of water molecules than comparable simulations employing methane dimers or butane molecules. In order to investigate structural changes among solvent molecules near the solute, we calculated the radial distribution function of water molecules with respect to the methylene units on the alkane chain, gco(r). The pair-distribution function was calculated as
where eo is the number density of water, and we have performed an ensemble average over all configurations for a given constrained end-to-end distance. In Figure 6 this distribution is given for an extended, as well as for a compact, chain conformation. For both conformational states a well developed shell exists with a first maximum at -4 A and and minimum at 5.5 A, beyond which there is little structure. This indicates that the alkane molecule is influencing only water molecules that are at the alkane/water interface. Long-range effects are not observed in the water density beyond the first water shell. Even though the shell is present for both conformations, water density is depleted for the compact states. As the alkane chain folds on itself, shell waters are being replaced by other methylene units from the chain. The number of actual water molecules within the first solvation shell of the dodecane molecule, nn,, was also monitored as
where the integration goes to the first minimum of the methylene-oxygen pair-distribution function, and eo is the number density of water molecules. In Figure 7 the number of nearest-neighbor water molecules, nn,, is given as a function of end-to-end distance. There is an approximate average
J. Phys. Chem., Vol. 99, No. 35, 1995 13123
Free Energy of n-Dodecane in Water in ‘ V
I
I
65
60
6 55
50
45
6
4
8
10
12
14
Tm/A
Figure 7. Number of nearest-neighbor water molecules nn, as a function of the chain end-to-end distance ree. Although characterized by large fluctuations,there is a significant decrease of shell waters as the chain becomes more compact.
0’4
0.2
l----l 1
-I
b c B
:0
8. A h
5
$ -0.2
> 6 a h I
-0.4
-0.6
’
4
6
Tc,d.A
10
Figure 4,the solvent contribution to the overall free energy is -6 kJ/mol, or, in terms of a contribution per water molecule associated with the solvent shell around the chain, -0.1 (kJ/ mol)/nn,. For noncompact states this number is even smaller, indicating a nonlinear dependence of the solvent contribution to the free energy of the system. The changing number of water molecules at the interface also affects energetic properties of shell waters. Thus, shell waters interacting with all other water molecules exhibit a 0.40 (kJ/ mol)/water more repulsive potential energy for compact states versus the fully hydrated case of an extended chain. This additional strain on the surface-exposed water molecules is to some degree compensated by additional potential energy developed between methylene units on touching strands of the compact chain. A further cooperative effect is that for compact states below re, < 9 8, the average water-alkane interaction is 0.03 (kJ/mol)/water more attractive than for the extended states. In order for such small changes among molecular properties to be noticeable and significant for the overall system of interest, a large number of molecules need to be affected. 3.2.3. Changes in Solvent Exposed “Area” During Hairpin Formation. It is of interest to investigate the correspondence between the solvation free energy, AAS01.,and the amount of surface exposure of the alkane chain. The solvent exposed surface area has been very successful in empirically accounting for many thermodynamic properties of both hydrophobic and hydrophilic molecules, as well as for larger biomolecules.40-42 While there is no theoretical justification for such an appr~ach,”~ it is highly desirable to quantify such a relationship from an “exact” model. Computer simulations yield an ensemble of configurations that are representative of the equilibrium system and can thus be used to calculate mechanical properties, such as an average area for the system. A popular view used to relate solvent effects to properties of the solute is to examine the amount of either buried solute area or solvent-exposed area44,45 to the free energy of the system,
I
12
14
Figure 8. Entropic contribution of the solvent to the free-energy change of bending the n-dodecane molecules on a per water molecule basis. The values of -ThS,,~/nn, versus reeexhibit large differences between compact (Tee < 9 A) and extended (ree> 9 A) states.
decrease of 10 water molecules associated with the collapse of the hydrocarbon molecules from an extended state to the compact state at 4 8,. The error bars indicate that this is a dynamic system where fluctuations of two to three water molecules in the shell are common. The number of nearestneighbor water molecules correlates with the enthalpic and entropic changes of AA,,I, in Figure 5 . Thus, for an end-toend distance larger than 9 A there is very little change in the number of water molecules closely associated with the hydrocarbon chain. Similarly the entropic penalty of maintaining this solvent shell is mainly constant and does not change significantly until the re,-values drop below 9 8,. However, below these distances the free energy is no longer strictly proportional to the number of water molecules in the shell. Once we go from a fully hydrated methylene unit in the alltrans conformation to a partially hydrated one, the entropic contribution per water molecule in the shell is no longer constant. As shown in Figure 8 the quantity -TAS,,hn, depends strongly on the end-to-end distance, especially for the compact states below re, 9 A where there is a large variation. This result indicates that free-energy changes depend on the solvent environment around a particular solute conformation. At the solvent-induced local minima, located at re, 4.2 8, in N
-
where the summation runs over all surface-area patches -4 that are associated with an energy/area coefficient ui. It is empirically found that this relationship can account for many experimental observations and can be useful for understanding complex phenomena where a fully molecular or atomistic treatment of the problem is not possible. Computer simulations of model systems designed to probe the relationship between macroscopic properties of water and the actual area change have indicated that the local behavior of the solvent near molecular solutes is critical and basically unaccounted for by simple area
consideration^.^^.^' In this work we are modeling a chemically realistic hydrocarbon chain in order to gauge the free-energy behavior with respect to changes in solvent/solute properties. The bending of the n-dodecane molecule into a hairpin conformation is one of the largest possible area deformations the chain can undertake without breaking apart. Given the numerous definitions of the appropriate area available, we have investigated two different methods to characterize the solvent surrounding the solute. The molecular surface, f,,,l, of the dodecane molecule wasocalculated by assigning a van der Waals (vdW) radius of 1.7 A to all methylene units and calculating re-entrant surfaces via a 1.4-8, probe.45 An excluded area, -f,,,,, was derived from the region of zero water density around any carbon atom (see below); Le. the excluded area is equivalent to a molecular surface with a carbon radius of 2.9 8, and zero probe radius-neglecting reI
Wallqvist and Cove11
13124 J. Phys. Chem., Vol. 99, No. 35, I995 400
380
35
360
-i $
0
4
4
e
01
e
340
320 185
200
I95
190
205
A3?,,,i/A2
Figure 9. Close relationshipbetween the investigated areas, displayed
as a scatter diagram of the excluded area {,,Iversus the molecular area Go,.Both of these areas have considerable fluctuations associated with them, as indicated by the box errors for a few data points. ~
~
4. Conclusions
6 4 0
2 w
E o .. Y rr
> + -2 a
e 0
-4
e 0
-6
appears to be no theoretical ground for such an a s s e s ~ m e n t . 4 ~ . ~ - ~ ~ Our calculated proportionality constants do, however, span both the experimental macroscopic value of 0.3 (kJ/mol)/A2as well as the two proposed “microscopic” tensions of 0.1 and 0.2 (kJ/ mol)/A2.40 Distributed among the first shell waters, the penalty for exposing any one water molecule to the average hydrocarbon interface is similarly 0.69 (kJ/mol)/nn,. As discussed in the previous section, there is an even stronger correlation between the free-energy penalty of maintaining an aqueous/alkane interface and the character of that interface, rendering linear relationships like that of eq 14 inoperative. The variation in surface area upon bending the chain molecule is relatively small compared to the area of the entire system. The volume occupied by the n-dodecane molecule compared to the entire system volume is likewise not changed significantly; consequently the pressure-volume work performed by the bending motion of chain is negligible. Variations in pressure during the simulation are less than 100 atm, and thus the results presented should be identical to those performed in an NFT ensemble.
v
I
e
0
.Q
185
195
190
200
205
A,,,,JA2
Figure 10. Solvent contribution to the free-energy change, as a function of the molecular area . o,{~. Even though there is a correlation between these quantities, there is a large amount of scatter in
the data. Illustrative errors have been indicated as boxes for selected data. entrant surfaces. There is also a close correspondence between the area measures and the number of water molecules associated with the chain, nn,, given in Figure 7. The excluded area -fexl and the molecular area -fmolare strongly correlated, as shown in Figure 9, where fexl 4.1 wG,ot. Likewise the excluded area is proportional to the number of water molecules in the first shell, -fexl 4.7 nn,. The simulations and the natural fluctuations of the hydrocarbodwater system do not allow discrimination between these measures. Consequently there is a rather broad range of area definitions that are not mutually exclusive and that allow for a wide latitude of representations and interpretations. As an example, the solvent contribution to the free energy, AAsol, is given as a function of the molecular surface area, fmol,of the dodecane molecule in Figure 10. There is a roughly linear relationship between the free-energy change and the area, but given the fluctuations of the system, the relationship is not quantitative and cannot reliably distinguish between states that are separated by less than -2.5 kJ/mol. The proportionality constant in eq 14 is 0.60 (kJ/mo1)/A2; if we instead employ the excluded area, -lex], the coefficient changes to 0.15 (kJ/mol)/A2. Although such numbers bear a superficial resemblance to interfacial surface tension coefficients, there I
-
-
In this work we have studied the detailed contributions of solvent to the overall free-energy landscape of a flexible hydrocarbon chain. The solvent component of the free energy is small when considered on a per solvent molecule basis. Thus, it is necessary that a substantial number of water molecules be involved before hydrophobic effects become the determining factor in establishing the equilibrium properties of the system. The solvent contribution is also dependent not only on the number of solvent molecules directly associated with the water/ alkane interface but also on the state of the alkane molecules-in our case a large difference exists between extended and compact chain conformations. Such dependencies indicate that a general formulation of free energies as being solely proportional to an exposed surface area is not feasible. Yet it cannot be denied that the empirical formulation of thermodynamic properties in terms of areas remains quantitative for comparing a select set of systems. For the short hydrocarbon chain studied here the most probable conformation is always an extended one, whether in the gas phase or in aqueous solution. The shape of the freeenergy curve as a function of end-to-end distance in Figure 1 shows a slight abatement in the free-energy penalty for compact states, a stabilization due to the increased number of methylenemethylene contacts between the chain ends. Given enough such favorable contacts, the entropic penalty of confining hydrocarbon strands to interlock will eventually be overcome for a large enough chain. Simulation of a hydrocarbon chain of 60 monomers shows that on average this system can support two to three hairpin bends in the gas phase. In this case additional enthalpic gains occur by forming interchain contacts among more than two strands. In the aqueous phase the solvent contribution to the stability of compact structures will be substantial, dominating the behavior of such chains in solution. In such a system the alkane “strain” at maintaining connected strands in such a fashion is largely due to entropy. For a larger chain this entropy penalty will grow, counteracted by the solvent gain in free energy from the release of entropy upon a minimization of the surface exposed by the hydrocarbon molecule. This delicate balance will determine the most probable conformational state of the hydrocarbon chain.
Acknowledgment. A.W. thanks Dr. R. Batik for insightful discussions on hydrophobic interactions.
Free Energy of n-Dodecane in Water
References and Notes (1) Kauzmann, W. Some factors in the interpretation of protein denaturation. Adv. Prorein Chem. 1959, 14, 1-63. (2) Blokzijl, W.; Engberts, J. B. F. N. Hydrophobic effects. Opinions and facts. Angew. Chem., Int. Ed. Engl. 1993, 32, 1545-1579. (3) Dill, K. A. Dominant forces in protein folding. Biochemistp 1990, 29, 7133-7155. (4) Young, L.; Jemigan, R. L.; Covell, D. G. A role for surface hydrophobicity in protein-protein recognition. Prorein Sci. 1994,3, 717729. ( 5 ) Pangali, C.; Rao, M.; Beme, B. J. A Monte Carlo simulation of the hydrophobic interaction. J. Chem. Phys. 1979, 71, 2975-2981. (6) Pangali, C.; Rao, M.; Berne, B. J. Hydrophobic hydration around a pair of apolar species in water. J . Chem. Phys. 1979, 71, 2982-2990. (7) Watanabe, K.: Andersen, H. C. Molecular dynamics study of the hydrophobic interaction in an aqueous solution of Krypton. J. Phys. Chem. 1986, 90, 795-802. (8) Tobias, D. J.; Brooks, C. L., 111. Calculation of free energy surfaces using the methods of thermodynamic perturbation theory. Chem. Phys. Leti. 1987, 142, 472-476. (9) Ravishanker, G.; Mezei, M.; Beveridge, D. L. Monte Carlo computer simulation study of the hydrophobic effect. Faraday Symp. Chem. SOC.1982, 17, 79-91. (IO) Guillot, B.; Guissani, Y.; Bratos, S. A computer-simulation study of hydrophobic hydration of rare gases and of methane. I. Thermodynamics and structural properties. J . Chem. Phys. 1991, 95, 3643-3648. (1 1) Smith, D. E.; Zhang, L.; Haymet, A. D. J. Entropy of association of methane in water: A new molecular dynamics computer simulations. J . Am. Chem. Soc. 1992, 114, 5875-5876. (12) Smith, D. E.; Haymet, A. D. J. Free energy, entropy, and internal energy of hydrophobic interactions: Computer simulations. J . Chem. Phys. 1993, 98, 6445-6454. (13) van Belle, D.; Wodak, S. J. Molecular dynamics study of methane hydration and methane association in a polarizable water phase. J. Am. Chem. Soc. 1993, 115, 647-652. (14) Dang. L. X. Potential of mean force for the methane-methane pair in water. J . Chem. Phys. 1994, 100, 9032-9034. (15) Sloan, E. D. Claihrafe Hydrates of Natural Gases; Dekker: New York, 1990. (16) Zichi, D. A.: Rossky, P. J. The equilibrium solvation structure for the solvent separated hydrophobic bond. J . Chem. Phys. 1985, 83, 797808. (17) Zichi, D. A.; Rossky, P. J. Solvent molecular dynamics in regions of hydrophobic hydration. J. Chem. Phys. 1986, 84, 2814-2822. (18) Pratt. L. R.; Chandler, D. Theory of the hydrophobic effect. J . Chem. Phys. 1977, 67, 3683-3704. (19) Rebertus, D. W.; Berne, B. J.; Chandler, D. A molecular dynamics and Monte Carlo study of solvent effect on the conformational equilibrium of n-butane in cc14. J. Chem. Phys. 1979, 70, 3395-3400. (20) Pratt, L. R.; Rosenberg, R. 0.;Berne, B. J.; Chandler, D. Comment on the structure of a simple liquid solvent near a n-butane solute molecule. J . Chem. Phys. 1980, 73, 1002-1003. (21) Rosenberg. R. 0.;Mikkilineni, R.; Berne, B. J. Hydrophobic effect on chain folding. The trans to gauche isomerization of n-butane in water. J. Am. Chem. Soc. 1982, 104, 7647-7649. (22) Jorgensen, W. L. Monte Carlo simulation of n-butane in water. Conformational evidence for the hydrophobic effect. J . Chem. Phys. 1982, 77, 5757-5765. (23) Zichi, D. A.; Rossky, P. J. Molecular conformational equilibria in liquids. J . Chem. Phys. 1986, 84, 1712-1723. (24) Tobias, D. J.; Brooks, C. L., 111. The thermodynamics of solvophobic effects: A molecular dynamics study of n-butane in carbon tetrachloride and water. J. Chem. Phys. 1990, 92, 2582-2592. (25) Hermann, R. B. Theory of hydrophobic bonding. I. The solubility of hydrocarbons in water, within the context of the significant structure theory of liquids. J. Phys. Chem. 1971, 75, 363-368.
J. Phys. Chem., Vol. 99, No. 35, 1995 13125 (26) Hermann, R. B. Theory of hydrophobic bonding. 11. The correlation of hydrocarbon solubility in water with solvent cavity surface area. J . Phys. Chem. 1972, 76, 2754-2759. (27) Andersen, H. C. Rattle: A ‘velocity’ version of the shake algorithm for molecular dynamics calculations. J. Comput. Phys. 1983, 52, 24-34. (28) Swope, W. C.; Andersen, H. C.; Berens, P. H.; Wilson, K. R. A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules: Applications to small water clusters. J . Chem. Phys. 1982, 76, 637-649. (29) Wallqvist, A,; Teleman, 0. Properties of flexible water models. Mol. Phys. 1991, 74, 515-533. (30) Wallqvist, A.; Beme, B. J. Effective potentials for liquid water using polarizable and nonpolarizable models. J. Phys. Chem. 1993, 97, 1384113851. (31) Berendsen, H. J. C.; Grigera, J. R.: Straatsma. T. P. The missing term in effective pair potentials. J. Phys. Chem. 1987, 91, 6269-6271. (32) Watanabe, K.; Klein, M. L. Effective pair potentials and the properties of water. Chem. Phys. 1989, 131, 157-167. (33) Weiner, S. J.; Kollman, P. A.; Case, D. A.; Singh, U. C.; Ghio, C.; Alagone. G.; Profeta, S.; Weiner, P. A new force field for molecular mechanical simulation of nucleic acids and proteins. J. Am. Chem. Soc. 1984, 106, 765-784. (34) Jorgensen, W. L.; Madura, J . D.; Swenson, C. J. Optimized intermolecular potential functions for liquid hydrocarbons. J. Am. Chem. Soc. 1984, 106, 6638-6646. (35) Hams, J. G. Liquid-vapor interface of alkane oligomers. Structure and thermodynamics from molecular dynamics simulations of chemically realistic models. J. Phys. Chem. 1992, 96, 5077-5086. (36) Zwanzig, R. W. High-temperature equation of state by a perturbation method. I. Nonpolar gases. J . Chem. Phvs. 1954, 22, 1420-1426. (37) Beveridge, D. L.; DiCapua, F. M. Free energy via molecular simulations: Applications to chemical and biomolecular systems. Annu. Rev. Biophys. Biophys. Chem. 1989, 18, 431-492. (38) Postma, J. P. M.; Berendsen, H. J. C.; Haak, J. R. Thermodynamics of cavity formation in water. Faraday Symp. Chem. Soc. 1982, 17, 5567. (39) McQuarrie. D. A. Siatistical Mechanics; Harper Row: New York, 1976. (40) Nicholls, A.; Sharp, K. A.: Honig, B. Protein folding and association: Insights from the interfacial and thermodynamic properties of hydrocarbons. Proteins: Struct., Fund., Gen. 1991, I I , 281-296. (41) Makhatadze, G. I.; Privalov, P. L. Contribution of hydration to protein folding thermodynamics. I. The enthalpy of hydration. J . Mol. B i d . 1993, 232, 639-659. (42) Privalov, P. L.; Makhatadze, G. I. Contribution of hydration to protein folding thermodynamics. 11. The entropy and Gibbs energy of hydration. J . Mol. Biol. 1993, 232, 660-679. (43) Ben-Naim, A,; Mazo, R. M. Size dependence of the solvation free energies of large solutes. J. Phys. Chem. 1993, 97, 10829- 10834. (44) Richards, F. M. Areas, volumes, packing and protein structure. Annu. Rev. Biophys. Bioeng. 1977, 6, 151-176. (45) Connolly, M. L. Solvent-accessible surfaces of proteins and nucleic acids. Science 1983, 221, 709-713. (46) Wallqvist, A,: Beme, B. J. Molecular dynamics study of the dependence of water solvation free energy on solute curvature and surface area. J . Phys. Chem. 1995, 99, 2885-2892. (47) Wallqvist, A,; Beme, B. J. Computer simulation of hydrophobic hydration forces on stacked plates as short range. J . Phys. Chem. 1995, 99, 2893-2899. (48) Ben-Naim, A. Solvation: from small to macro molecules. Curr. Opin. Sfruct. Biol. 1994, 4 , 264-268. JP95092 I +
