Hydrophobic Interactions from Surface Areas, Curvature, and

Chapter 23

Hydrophobic Interactions from Surface Areas, Curvature, and Molecular Dynamics Use of the Kirkwood Superposition Approximation To Assemble Solvent Distribution Functions from Fragments 1

Downloaded by MONASH UNIV on October 26, 2012 | http://pubs.acs.org Publication Date: September 29, 1994 | doi: 10.1021/bk-1994-0568.ch023

Robert B. Hermann

Lilly Research Laboratories, E l i Lilly and Company, Indianapolis, IN 46285

A previously developed methodology for the calculation of hydrophobic interactions based on resolving solution-free energies into solvent cavity potentials and solute-solvent interaction energies is investigated further. The aromatic systems benzene and toluene now are well fitted along with the aliphatic compounds. The inclusion of the local mean curvature of the accessible surface is examined. The solvent effect on cyclohexane dimerization is found to be unfavorable, in accord with previous estimates. The calculations rely on a molecular dynamics determined solute-solvent interaction energy. However, it is possible to build the solute-solvent distribution function for a molecule or dimer from molecular fragment distribution functions via the Kirkwood superposition approximation. Distribution functions built from such transferable solvent distribution functions, in the case of methane dimer and ethane, give encouraging results for the solute -solvent interaction energy.

Approximate intermolecular potentials and semi-empirical methods are widely used today in the calculation of molecular and bulk thermodynamic properties. Solvent effects models including hydrophobic interaction models are useful in conjunction with such methods for the semi-empirical calculation of free energies and dissociation constants in solution (1-11). In this regard, the development of a working model for the calculation of hydrophobic interactions is continued here. In a previous paper (4), a method was developed which allows one to calculate the hydrophobic interactions between small hydrocarbon molecules. In that method, two important quantities are solvent cavity potentials and solute-solvent interaction energies. Solute-solvent interaction energies are found from molecular dynamics simulations, and cavity potentials are found by calibrating a parameterized version of scaled particle theory with experimental hydrocarbon solvation energies. 1

Current address: Department of Chemistry, Indiana University/Purdue University at Indianapolis, 402 North Blackford Street, LD3326, Indianapolis, IN 46202-3274

0097-6156/94/0568-0335S08.72/0 © 1994 American Chemical Society

In Structure and Reactivity in Aqueous Solution; Cramer, C., et al.; ACS Symposium Series; American Chemical Society: Washington, DC, 1994.


33é

STRUCTURE AND REACTIVITY IN AQUEOUS SOLUTION

Here this methodology is applied to several other interesting eases. In addition, the use of local curvature (12) is investigated rather than just the use of inherent curvature of spherical cavities as was done before. Some current methods for continuum solvation effects treat different regions of the accessible surfaces (13) with different parameters. Local curvature considerations are important for treating portions of accessible surfaces (12), e.g. in the vicinity of hydrophobic groups, since these have their own associated curvature. In addition to the development of a means to obtain cavity energies, a simpler method currendy under development to obtain interaction energies within the scope of continuum methods is outiined; namely, an attempt is made to construct a solvent distribution function for a larger solute from the distribution functions of smaller molecules, or fragments. While the solvent distribution function of the fragments may be obtained from molecular dynamics calculations, the solvent distribution function associated with the larger system is obtained by putting these fragment distribution functions together from a set of rules. Theory The solubility of hydrocarbons in water has been treated by breaking down the solvation-free energies into two main contributions-the cavity potential and solutesolvent interaction energy (14-20). The following equation for dilute solutions relates the solvation potential to Henry's law constant: k T l n K = μ* + kT In NrjkT/V

(1)

H

where K H is Henry's law constant, k is Boltzman's constant, Τ is the temperature, No is Avogadro's number and V is the molar volume of the solvent. The method as developed in a previous paper (4) and outlined here is based on the idea that the solvation potential μ* can be expressed as μ* = μο + βΐ

(2)

where ei is defined as the following average: ei=

(3)

0

where N is the number of water molecules and the solute subscript is 0. μο, the cavity potential, is then defined by the difference w

μ = μ*-βί

(4)

€

By this definition, contains all entropie contributions to μ* and ei contains none. The quantity ei may be found then by averaging molecular dynamics runs using a hydrocarbon whose coordinates are fixed and ΊΤΡ3 water (21). Fixed conformations of hydrocarbons are used and each must be calculated individually. The energy was calculated out to 15 Â so that a correction 2

Ecoix = 4πρι Juo r * W

Rc


(5)

23. HERMANN

337

Hydrophobic Interactions from Surface Areas

where Rc = 15 Â must be applied. The results of the molecular dynamics calculations are shown in Table I. The column labeled ei is the sum of the dynamics result plus EcorrThe calculation of the cavity potential is done by first measuring the accessible surface area (13) of the molecule or system. The solute-solvent thermal radii (22) are used which simplifies the parameter choice problem somewhat. After carrying out such calculations for 13 hydrocarbons, the resulting calculated cavity potentials were fitted to the following equation (4) suggested by scaled particle theory (77, 18, 23, 24):

μ = aA(l + b/r + c/r ) 2

(6)

0


2

where r is found from the accessible surface area A = 4itr , and a is approximated by the surface tension of water (25), 103.5 cal Â , b is 3.09 Â and c is 3.25 Â . The parameters a, b and c are, in general, temperature dependent; however, they have been found by calibrating against solubility data at a particular temperature (4). To calculate hydrophobic interactions between pairs of hydrocarbon molecules, the molecular dynamics calculations are performed on the (associated) pair to get ei and the cavity potential μο for the pair is found from the area of the pair and Equation 6. The curvature is taken into account through the term b/r. Table II lists the areas and Table III gives the resulting cavity potentials. Figure 1 shows the cavity potentials of the monomers plotted against the accessible surface areas. The cavity potential as given by calibrating Equation 6 is dependent of the nature of the molecules occupying the cavity. The above calibration should represent cavity potentials for aliphatic systems generally. 2

2

Combined Treatment of Aliphatic and Aromatic Systems It would be desirable to treat aliphatic and aromatic systems the same way, i.e. with the same set of cavity parameters. Benzene using OPLS parameters (26) gives an interaction energy of -14.132 kcal m o l . This value, together with the accessible surface area of 230.5 and an experimental solvation energy of -883 cal mo\ (27) produces a cavity potential from Equation 6 much too high to fit on the cavity potential vs. area curve of Figure 1. Presumably because of its larger quadrupole moment, benzene interacts more strongly with water than an aliphatic compound of similar molecular weight. Because of this extra electrostatic interaction, the cavity formed by aromatic systems is perturbed relative to the cavity formed by aliphatic molecules. The resulting cavity energy would be significantly higher in energy per unit area from the aliphatic cavities. Therefore, die parameters for Equation 6, defining the cavity potential for a solute, would not apply to aromatic systems. A solution to this problem applied here was to carry out the molecular dynamics simulations on benzene using aliphatic charges rather than aromatic charges but otherwise using aromatic OPLS parameters. This insures that the cavity surface has a similar structure to that formed by the aliphatic molecules, so that the cavity energy per unit area is comparable. When evaluating the interaction energy ei afterward by averaging over all the configurations, the benzene molecule is then given the aromatic charges. 1

l

The energy is found from N

W

ei = < Z Uqj >ai

(7)


338


Table I Molecular Dynamics Interaction Energies and Related Data System

(MD)* (kcal mol )


1

Methane Ethane Propane trans «-Butane gauche n-Butane Isobutane trans n-pentane gauche n-pentane Neopentane Cyclopentane Dimethylbutane Cyclohexane Dimethylpentane Cycloheptane Isooctane Benzene Toluene

b

0

0

C H 4 - C H 4 0.0 A C H 4 - C H 4 1.54 A CH4 - C H 3.45 A CH4 - C H 4.0 A CH4 - CH4 5.0 A CH4 - C H 6.0 A CH4-CH47.I6A CH4-CH48.OA 4

D

4

4

CH4-C H C H -C H /-QHio-i-QHîo /i-C Hi2-n-C Hi2 CôH^-QHn 2

2

6

6

2

5

6

5

6

QH12- Q H i 2 CÔHÔ - CéHg (C H )3 2

6

(C H ) 2

6

4

f

-3.172 -5.186 -7.086 -8.670 -8.702 -8.654 -10.630 -10.242 -9.873 -9.786 -11.578 -11.120 -12.914 -12.618 -14.151 -11.664 -13.325 -7.070 -6.380 -5.766 -5.850 -5.866 -6.081 -6.383 -6.363 -7.800 -9.340 -14.989 -17.757 -18.900 -19.959 -20.911 -13.270 -15.864

ei (kcal mol ) 1

-3.216 -5.264 -7.198 -8.815 -8.847 -8.799 -10.808 -10.420 -10.052 -9.954 -11.790 -11.321 -13.159 -12.852 -14.430 -11.841 -13.536 -7.159 -6.469 -5.855 -5.939 -5.955 -6.170 -6.472 -6.452 -7.922 -9.496 -15.279 -18.114 -19.303 -20.361 -21.265 -13.504 -16.176

Run Waters Length (ps) 180 180 300 210 180 180 180 180 180 180 180 180 180 180 300 300 270 300 240 210 248 240 300 270 270 210 180 180 180 300 270 270 300 330

926 960 1025 1074 1136 1111 1118 1167 1171 1100 1223 1165 1247 1207 1282 1045 1137 963 1014 1082 1082, 335 1136 1150 1160 1239 1096 1096 1294 1330 1358 1376 1356 1220 1321

a

Aliphatic monomers from Ref. 4. Conformation with one gauche interaction. Interaction energy calculated from Equation 7. 48 ps with 1082 waters and 200 ps with 335 waters. Parallel or stacked configuration. See Figure 2b. ^Perpendicular or "T" configuration. See Figure 2c. c

d

c


23.

HERMANN


339

Table Π. Geometric Features of the Cavitv Surfaces Average Mean Average Squared Accessible Mean Curvature^* Curvature * Surface Area^

System*

0

2

K (Â-2)


(A*)

Methane Ethane Propane trans n-Butane gauche Λ-Butane Isobutane trans n-Pentane gauche n-Pentane Neopentane Cyclopentane Dimethylbutane Cyclohexane Dimethylpentane Cycloheptane Isooctane Benzene Toluene Œ4-CH4O.OÂ

f

CH4-CH4 1 . 5 4 Â CH4-CH4345A

CH4-CH44.OÂ CH4-CH45.OÂ CH4-CH46.OÂ CH4-CH47.I6A GH4-CH48.OÂ CH4-C2H6

C H 2

C H

6

2

6

i-QHio - I-QHJO

n-C Hi2-n-C Hi2 5

5

QH12- QH128 C6H -C6H h

12

12

CôHg - CeHô (C H )3 (C H ) 2

2

6

6

4

M

135.5 173.2 204.4 235.3 231.1 229.6 266.1 262.0 250.7 236.8 273.2 258.4 306.2 278.5 321.5 230.5 259.4 148.1 167.1 207.5 222.5 244.0 270.9 271.7 271.0 244.1 266.8 358.5 393.4

0.312 0.278 0.259 0.244 0.243 0.243 0.232 0.230 0.232 0.240 0.223 0.230 0.215 0.221 0.208 0.243 0.232 0.295 0.281 0.259 0.253 0.245 0.239 0.285 0.304 0.242 0.232 0.199 0.191

407.0

0.186

0.099 0.081 0.071 0.064 0.064 0.063 0.060 0.059 0.059 0.063 0.055 0.057 0.053 0.053 0.050 0.064 0.060 0.088 0.081 0.072 0.071 0.071 0.076 0.088 0.096 0.066 0.061 0.050 0.045 0.047

430.7 379.3 357.2 383.7

0.185 0.198 0.205 0.194

0.045 0.052 0.055 0.048

a

See text for dimer configurations. bCalculated using M O L A R E A , Quantum Chemistry Program Exchange 225, Indiana University, Bloomington, IN. Mean curvatures at a point were measured by the method of Nicholls et al.; Ref. 12. A 2.8 radius probe was used. dThe mean curvature at a point was squared to get the squared mean curvature at that point. Averages are for the entire surface. tonly one gauche interaction. ^Parallel or stacked configuration. See Figure 2b. ^Perpendicular or "T" configuration. See Figure 2c. c

e



340


Figure 1. Cavity potentials from equation 4 for aliphatic hydrocarbons, toluene and benzene.


23. HERMANN

Hydrophobic InteractionsfromSurface Areas

341

where the subscript al indicates that the configurational average is over the aliphatically charged cavity and not die aromatic charged cavity. Uqj is the complete aromatic molecule-water molecule interaction potential. Thus the difference between the aromatic and aliphatic charge-solvent interaction is treated as a first order perturbation of a zeroth order dynamics distribution function. With the ei calculated in this way, benzene and toluene are fitted reasonably well, as can be seen from Figure 1. For toluene, OPLS aromatic parameters are used for the phenyl ring and OPLS aliphatic parameters (4) are used for the methyl group. The charges on the methyl group carbon and the phenyl group carbon to which it is joined are then adjusted to give a dipole moment of 0.37 Debye (28) with the positive end of the dipole in the direction of the methyl group. Downloaded by MONASH UNIV on October 26, 2012 | http://pubs.acs.org Publication Date: September 29, 1994 | doi: 10.1021/bk-1994-0568.ch023

Inclusion of Local Curvature In the calculation of He in Equation 6, the curvature of the cavity is taken into account implicitly since it is found from the accessible surface area of the molecule as the inverse of the radius of the sphere which has the same area as the accessible surface. This corresponds to the second term in Equation 6 involving the parameter b. The assumption is made that the energy of the cavity formed by the molecule is the same as that for a spherical cavity of equal surface area (4,20,29). This equivalent sphere method of representing the cavity is necessary i f the liquid state theory used to calculate the cavity energy is applicable only to spherical cavities. In the equivalent sphere approximation, it is assumed that when the sphere is deformed to the accessible surface of the molecule, regions that were deformed such that curvature increased were somewhat compensated for by corresponding regions where curvature is thereby necessarily decreased or negative; since the cavity potential depends on curvature, such free energy changes tended to cancel. The good correlation with the data indicates that this assumption is reasonable at least for die accessible surfaces of small hydrocarbons. This equivalent sphere approximation was used previously in connection with Barker Henderson perturbation theory based on a Perçus-Yevick zeroth order hard sphere representation and produced results in qualitative agreement with this paper (20,29). Nicholls et al. (12), in their treatment of hydrophobic interactions through interfacial tension, have demonstrated the effect of local curvature on the interfacial free energy. They have used the accessible surface area and interfacial tension together with a modification due to local curvature. In this paper, the detailed curvature of the accessible surface cavity is similarly considered. The two descriptions of curvature of a surface are the mean curvature and the Gaussian, or total curvature (30). It is this mean curvature that is of interest here. The mean curvature K at a point Ρ on the smooth cavity surface is given by m

Κ

Λ

= (1*ι + 1Λΐ)/2

(8)

where Γχ and τ% are the principal radii at a point Ρ on the surface. For a sphere of radius r, this reduces to K = 1/r. The average value of the mean curvature is the value of the mean curvature at each point on the surface, averaged over the entire surface. This will be designated K . Another quantity which will be of interest is the average value of the square of the mean curvature K . The curvature at each point on the surface is squared to get the mean curvature squared at that point. This is then averaged over the surface. The mean curvature is measured at a point on the surface in the manner of Nicholls et al. (12). They calculate the local curvature at a point on the accessible surface of a molecule by determining the accessible surface of a (spherical) water molecule at that point For a planar accessible surface, one half of the water molecule m

m

2

m


342


is accessible but varies more or less for curved surfaces. Calculations of this accessibility determines the approximate curvature at that point on the surface. Table II gives the necessary geometrical parameters to be usai in subsequent calculations on molecules and complexes. Five thousand points per spherical surface were used in the calculation in this paper. Two alternatives to Equation 6 are now considered, as possible methods for the treatment of local curvature. These possibilities are given by Equations 9 and 10. In both cases in the curvature term, 1/r is replaced by , the average mean curvature. Hc = aA(l+b*ic; c7r2)

(9)

+


^ic = a A ( l + b * ^ + c " i C ^ )

(10)

In Equation 10, the last term involves the average of the square of the mean curvature instead of 1/r . The results for the inclusion of local curvature, based on Equations 9 and 10 is given in Table III along with the results of Equation 6. The agreement with experimental solubilities is only slightly better than the equivalent sphere method of Equation 6. The aromatic systems, given the special treatment mentioned above, can in all cases also be accommodated with good agreement. The parameter a is the surface tension of water and is the same for Equations 6, 9 and 10. The parameters b' and c' for Equation 9 were found to be 2.77 and 2.39, respectively while b" and c" for Equation 10 were found to be 3.09 and 3.25, respectively. The parameters for Equation 10 are the same as for Equation 6. 2

K

2 m

and the Dissociation of Dimers

In some cases, e.g. the dissociation of dimers, the solute system can consist of more than one molecule, thereby requiring two or more distinct cavities, i.e. two or more separate surfaces. In order to treat such composite solutes involving multiple surfaces, Equation 10, rather than Equation 6 or 9 is necessary. This is easily demonstrated for a collection of spherical cavities, as follows: First of all, for a cavity having a single spherical surface, the contribution to the potential due to the last term of Equation 10 is n

A a c " K ^ = 47iac .

(11)

This is independent of the cavity radius and is related to the work of introducing a single point particle into the solvent. In scaled particle theory, the cavity potential is derived from the process of introducing a point particle into the solvent, and then allowing it to grow up to the size of the final desired (spherical) particle. Introducing a point particle increases the chemical potential, and such a process should be associated with each distinct spherical surface. Equations 6 and 9, as in scaled particle theory, can account for the introduction of only one point particle. The important difference is that in addition to being able to describe a single surface, Equation 10, unlike Equations 6 or 9 can also describe several spherical surfaces. If the total area of the system is A = 4πτ and the system is made up of η molecules, then 2

2

A = 4π(η + r

2 2

2

+ ·· ·r ) n

(12)

where the molecules 1, 2, · · · η have the spherical curvatures K = l / n , l/r , · · · l/r . m

2

n

(13)


23.

HERMANN


Table III. Solvation Potentials and Cavity Potentials

System^

μα Médiane Ethane Propane trans n-Butane* gauche n-Butane Isobutane trans w-Pentane gauche n-Pentane Neopentane Cyclopentane Dimethylbutane Cyclohexane Dimethylpentane Cycloheptane Isooctane Benzene Toluene f


f

f

CH4 C H 4 O . 0 A C H 4 - C H 4 1.54 A C H - C H 4 3.45 A CH4-CH44.0A CH4-CH45.0A CH4-CH46.OA CH4-Œ47.I6A C H 4 - C H | 8.0 A 4

CH -C H C H6 - C H6 4

2

2

6

2

I-QHIO-I-QHÎO

n-C H -n-C H C6H -C6H g C6H -C6H h C^HÔ - CoHé (C H ) (C H ) 5

a

1 2

5

12

12

1 2

1 2

2

6

3

2

6

4

1 2

Λ

μ*(αχρ)

5190 1974 7069 1805 9143 1946 10968 2051 10691 2051 11029 2230 13117 2346 12840 2346 12693 2641 11169 1216 14245 2456 12554 1234 16008 2849 13655 803 17354 2924 10957 -883 12771 -765

Equation 6 /ic 5053 7231 9171 11189 10912 10811 13278 12997 12221 11288 13767 12748 16083 14129 17177 10869 12812 5760 6868 9372 10344 11769 13604 13663 13612 11779 13323 19865 22457 23482 25279 21406 19774 21730

1837 1968 1973 2373 2064 2013 2469 2576 2170 1334 1978 1427 2923 1276 2747 -971 -724 -1398 399 3516 4405 5814 7433 7191 7160 3857 3827 4585 4343 4179 4918 141 6270 5554

Equation 9 μ*« με 5029 7237 9098 11032 10934 10902 12939 12932 12347 11319 13909 12829 15923 14267 17233 10924 12702 5897 6937 9205 9974 11252 12588 9057 7518 11418 12984 19747 22321 23530 24903 20857 19095 21445

1813 1973 1900 2216 2087 2103 2131 2511 2296 1366 2119 1509 2764 1415 2803 -916 -834 -1262 469 3350 4035 5298 6418 2586 1066 3496 3488 4468 4207 4227 4542 -408 5591 5269

343

3

Equation 10 με 5042 7258 9150 11116 10932 10852 13141 13003 12249 11322 13801 12723 16100 14122 17268 10896 12801 5747 6834 9313 10282 12008 14271 11455 10399 11796 13339 20277 22718 24290 25732 21918 20142 22106

1

1826 1994 1953 2300 2084 2053 2332 2583 2198 1368 2012 1402 2941 1270 2839 -944 -729 -1412 365 3458 4343 6053 8101 4983 3947 3874 3842 4998 4604 4987 5371 653 6638 5930

In cal m o l . For dimer configuration, see text. Except for n-butane and n-pentane, this is found from Uç = μ*- ei, where μ* is the experimental solvation energy (column 2) and ei is from Table I. Aliphatic entries previously listed in Ref. 4; Benzene and Toluene values from Ref.27. Calculated using from the preceding column and e( from Table I. ΐ ο τ treatment when more than one conformation is present, see Ref. 4. ^Parallel or stacked configuration. See Figure 2b. Perpendicular or "T* configuration. See Figure 2c.

b

c

d

e

h


344


The average squared mean curvature consists of these η terms, each weighted with their corresponding areas and divided by the total area A : 2

= {4πτ 2(ΐ/η) + 4πτ 2(ΐ/ )2 + . . . + 4πτ 2(1/ )2 }A1

2

Γ2

η

1

Γιι

(14)

so that


Aac'lu

5

= 47cac"n

(15)

It can be seen that the quantity 4rcac" enters in once for each distinct spherical surface in the system. In this manner, Equation 10 can represent η spherical surfaces, or cavities. The expectation is that Equation 10 will represent any number of distinct solute surfaces. This is critical, for example, in following dimer dissociation such as methane-methane (see below). In order to get the correct description of a dimer in the limit of large separation, Equation 10 is necessary. Applications to Selected Molecular Pairs and Aggregates The numerical results for the molecular systems below are shown in Tables ΙΠ and IV. The choice of dimer configuration shown in Figure 2 is somewhat arbitrary and usually depends in part on maximum contact. The remaining dimer configurations were defined in ref. (4). Once the general dimer configuration was chosen, the energy was minimized using M M 2 or M M P 2 (31 ) on the isolated system, with the constraint that the chosen relative orientation be maintained. The solvation energies calculated for these dimers are for these fixed configurations and are given in Table III. The binding energies in Table IV are found from these fixed configuration solvation energies and the solvation energies of the monomers. Comparing the solvation energies of these selected dimers, while incomplete in the sense that all possible configurations were not considered, gives some insight as to whether the solvation effect augments or hinders association in solution. Cyclohexane Dimer. Two configurations were treated, as shown in Figures 2b and 2c. The relative positions were minimized using M M 2 , giving a 4.55 Â between ring centers for parallel and 4.65 A between centers for the perpendicular configuration. The parallel configuration of the isolated dimer had a binding energy of -3.29 kcal m o l , while the perpendicular configuration has a binding energy of -1.93 kcal m o l . The calculated solvation free energies associated with each indicate less association in water than in the gas phase. The values for both the equivalent sphere method and the explicit local curvature methods are qualitatively similar. The M M 2 energy and solvation energy together were then minimized, varying only the intermolecular distance (between ring centers). For the parallel configuration, a new minimum was found at 4.45 A, favoring a greater association by 190 cal m o l over the results of Table IV. In the case of the perpendicular configuration, a similar minimization placed the molecules .2 A closer and decreased the energy by 172 cal m o l . Subtracting this from the results in Table VI, solvent effects still favor dissociation of the dimer. Based on the results for the one parallel conformation only, the equivalent sphere method predicts the dimer is favored 7 to 1 in the gas phase over solution, i.e. the solvent effect favors dissociation. Second virial coefficient data (32) and osmotic coefficients calculated for cyclohexane has suggested a dimer concentration three times greater in the gas phase than in solution (33). There is a reduction in both cavity area and cavity energy upon association. Usual methods for calculations of changes in hydrophobic interactions would give a reduction in area and would, therefore, incorrectly predict a greater association in solution. 1

1

1

1


23.

HERMANN


Table IV. Calculated Hydrophobic Interactions Pair

(Equation 6) Binding potential


CH4-CH4

0.0 A

1.54 3.45 4.0 5.0 6.0 7.16 8.0 CH4 - C H G 2

C H -C H 2

6

2

6

1-C4H10 - Î-C4H10

n-C Hi -n-C H C H -C6H c C H -C6H d CeHé - Cells C H -(C H ) 5

2

6

1 2

6

1 2

2

6

5

1 2

1 2

2

6

3

1 2

-5072 (-5347) -3274 (-3549) -158 (-432) 731 (456) 2140 (1866) 3760 (3485) 3517 (3243) 3486 (3212) 52 (77) -109 (217) 559 (125) -664 (-348) 1324 (1711) 2064 (2451) 2081 (1909) -2683

a

1

(Equation 9) Binding potential * 1

-4888 -3157 -276 409 1672 2793 -1040 -2560 -289 -458 262 -306 1209 1525 1425 -2295

(-5211) (-3479) (-598) (86) (1349) (2470) (-1363) (-2883) (-283) (-122) (8) (-484) (1759) (2075) (1360)

345

3

(Equation 10) Binding tial poten b

-5064 (-5360) -3288 (-3583) -195 (-490) 690 (394) 2400 (2105) 4448 (4152) 1330 (1034) 285 (-D 53 (95) -147 (232) 891 (538) -226 (-87) 2183 (2519) 2567 (2904) 2541 (2420) -2703

Reduction in area 123.0 103.9 63.5 48.5 27.0 .1 -.7 0.0 64.6 79.6 100.0 136.1 110.0 86.3 81.7 146.7

2

Binding potentials in cal mol" . Areas in Â . ^Binding potential is dimer solvation potential minus the solvation potential of two monomers from Table ΠΙ. The first binding potential is based on the monomer calculated energies while the numbers in parentheses are based on the experimental monomer energies. Solute-solute interaction not included. From Figure 2b. From Figure 2c. c

d

Benzene Dimer. Table IV shows the results when benzene and the dimer are treated as discussed above. For both monomer and dimer, the molecular dynamics was carried out on the aliphatic charge model as outlined above and the interaction energy evaluated afterward with aromatic charges as shown in Equation 7. The dimer configuration that was suggested by Jorgensen and Severance (26) was used here and shown in Figure 2d. The distance between ring centers is 4.99 Â measured as indicated above. A s Table IV shows, the solvation energy favors dissociation. Again, this is qualitatively in accord with second virial coefficient data (32) and estimates of the osmotic coefficient (33, 34, 35). Based on this single configuration only, the solvent effect in all three cases is somewhat overestimated. The M M P 2 energy and solvation energy were minimized by varying the intermolecular distance along a line through the center of one ring and perpendicular to the plane of the other ring. The distance was thereby reduced by .2 A and the energy by 95 cal m o l . In spite of the small change, the overall effect of solvation is found to favor dissociation. 1

Ethane Dimer, Isobutane Dimer, n-Pentane Dimer and Methane-Ethane. Calculations for the equivalent sphere method were reported previously (4), but are included for comparison. The dimer configuration chosen (4) was one of the more


346


H IN / \1 -c—i—»,cNN


(a)

(b)

Is (c) 1

(f)

(e)

Figure 2. Dimer and aggregate configurations of hydrocarbons.

13 υ

-a eο

Eq.10

H3 .S β

Carbon-Carbon Separation (Â)

Figure 3. Methane-methane hydrophobic interaction according to Equations 6,9, and 10. In Structure and Reactivity in Aqueous Solution; Cramer, C., et al.; ACS Symposium Series; American Chemical Society: Washington, DC, 1994.

23.

HERMANN

347


closely packed possibilities to demonstrate an upper limit to the solvation effect. A slightly lessened contribution of the solvent effect toward dissociation was found upon the inclusion of local curvature as in Equation 9, but not in Equation 10. In all cases, a reduction in accessible surface area was observed upon association. Ethane Aggregation. To demonstrate that under certain conditions, the solvation effect can favor strong association, the interaction of an ethane molecule with an aggregate of three strategically placed ethanes was considered. Table IV gives the relevant binding energies. The configurations are shown in Figures 2e and 2f. It can be seen, in accord with previous results (29), that solution favors association by 2.32.7 kcal m o l . In this case, like the others above, association results in a reduction in area. Downloaded by MONASH UNIV on October 26, 2012 | http://pubs.acs.org Publication Date: September 29, 1994 | doi: 10.1021/bk-1994-0568.ch023

1

Methane-Methane. A range of intermolecular distances is considered for the methane-methane interaction. A plot of the association energy, based on the three methods for the calculation of the cavity energy, is shown in Figure 3. The advantage of the third term in Equation 10 involving the average squared mean curvature can now be easily seen. The solvent effect correctly goes to zero for large separations, as can be seen from the curve representing Equation 10. Equation 6, on the other hand, implies too low a curvature resulting in a high cavity energy upon separation. The last term in Equation 9 corresponds to only one point particle and gives too low an energy upon greater separation. Beginning with the M M 2 vacuum minimized distance of 3.45 A between the methane molecules, a new minimum at 3.35 A was obtained by minimizing again but including the solvent effect. This resulted in a small increase in dimer stability of 55 cal m o l over the results of Table IV. These results are for a closest approach interaction as shown in Figure 2a and neglect other methane-methane orientations at contact. 1

Modeling the Interaction Energy e|. A further simplification would be desirable so that a molecular dynamics calculation does not have to be done on each system for which ei is needed. It would be useful to model the interaction so that with a set of rules and parameters, ei for a molecule could be found, to sufficient accuracy. In order to obtain ei, the solvent distribution function must be known so that ei may be calculated from the integral over the 2

intermolecular potential weighted with the solvent distribution function (36), p:

(16)

t h

where Uk is the interatomic potential between the k solute atom and a water molecule, m is the total number of solute atoms and p = d g< is the water distribution about the solute atom, where g< ) is the solute-solvent correlation function, d$ is the solvent density and dt indicates integration over the relevant coordinates. The method under investigation is based on the idea of obtaining the solvent distribution function for a molecular "fragment", e.g. a methane molecule, and then just as these fragments may be combined to build a large molecule, the distribution functions associated with the fragments are combined to produce the large molecule 2

2)

s

2


348


n

distribution function. The discussion below is in terms of the correlation function g rather than the distribution function p< ). n

Methane Dimer-Water Correlation Functions. It is first instructive to examine the solute-solvent correlation functions of some systems for which there are molecular dynamics results. In Figures 4 and 5, the distribution of water oxygens about the methane dimer (inter-carbon distance = 7.16 A) as determined by a molecular dynamics calculation is plotted. These figures are different views of a cylindrical plot in which the methane-methane axis is taken as the cylinder axis x. The water molecule oxygen positions obtained from 270 ps dynamics run sampling 13,500 configurations using Amber 3.0A (37) were allocated to 27irx.25x.25 A bins over three dimensional space. Adding the waters around the cylinder axis in this manner gives a representation of the correlation function in two dimensional space. The .25 A grid spacing in the χ and r directions can be seen from the figures. The intensity of the correlation function is then plotted in the ζ (upward) direction. Referring to Figure 4, the position of the methanes is easily seen to be where the correlation function is zero. The methanes are oriented as shown in Figure 2a and a slight egg-shaped asymmetry may be noticed due to the axial hydrogens of the methanes pointing away from the center. The first peak in the correlation function is clearly indicated as a wall of oxygen atom density around the methanes. Since methane is not exactly cylindrically symmetrical, the first peak may be slightly low due to being slightiy smoothed out around the three non-axial hydrogens of each methane. A slight trough lies right behind the wall. After a very slight peak after this, the density beyond is fairly uniform. Examination of Figures 4 and 5 indicate interesting secondary features. First, a build-up of water between the methanes can be seen. The presence of a water molecule between two hydrophobic solutes has been pointed out previously (38, 39, 40, 41, 42). Second, a slight deepening of the trough of the correlation function behind the peak, can be seen equidistant from the methanes. This deepening can be seen in the center of Figure 5. Finally, there are two slight dips in the wall on either side of the peak, approximately at the position where the trough would intersect the crest. These features may be interpreted as an interference pattern between two methane centered methane-water correlation functions if a separate methane water correlation function is assumed to be present around each methane. However, the two assumed functions multiply rather than add. Figure 6 shows a simple methane-water correlation function, from a 510 ps molecular dynamics run on methane. Figures 7 to 9 show results of other dynamics runs. In general, there is an increase of density of water between the methanes where possible. In Figure 8 at 3.45 A, a slight increase in density can still be seen. Figure 9 shows the molecules coalesced to the unnatural distance of 1.54 A.


3

Kirkwood Superposition Approximation. In order to produce similar results from a synthetic approach, it is necessary to make use of the Kirkwood superposition approximation (43, 36). The superposition approximation gives the three body correlation function g< )(123), in this case two methanes, designated 1 and 2, and a water molecule, designated 3, as a product of three two body correlation functions: 3

2

g0>(123) = g< >(12)g(2)(13)g(2)(23)

(17)

If the methanes are held fixed, then the superposition approximation may be written in unsymmetrical form (43,36): g[31(123) = g(2)(13)g(2)(23)


(18)


23. HERMANN


Figure 4. Methane-methane-water correlation function for carbon-carbon distance of 7.16 A from a 270 ps molecular dynamics run.

Figure 5. Front view of Figure 4. Carbon-carbon axis horizontal.


349


350


Figure 6. Methane-water correlation function from a 510 ps molecular dynamics run. Blip in front center indicates position of carbon atom.

Figure 7. Methane-methane-water correlation function for carbon-carbon distance of 6 Â from a 300 ps molecular dynamics run.



HERMANN


Figure 8. Methane-methane-water correlation function for carbon-carbon distance of 3.45 Â from a 210 ps molecular dynamics run.

Figure 9. Methane-methane-water correlation function for carbon-carbon distance of 1.54 Â from a 240 ps molecular dynamics run.


352


2

Figure 6 is regarded as g< )(13) and the symmetry operation χ —» -x is performed on this to get g< >(23). If these correlation functions, namely g(13) and g< (23), are translated relative to each other so that the carbons are 7.16 A apart, and then the product g( >( 13) g< >(23) is formed, the result is shown in Figure 10. The product g< >(13) $ >(23) reproduces the true function (Figure 4) rather well, except for possibly an extra large build-up between the methanes. Of course, sampling on or very near the χ axis is much less adequate than at larger r values. The same procedure of obtaining a correlation function from the product of two methane water correlation functions was applied to the other methane-methane separations. The results are shown in Figures 11, 12, and 13. In all cases, the superposition of the peaks of the correlation functions lead to some build-up between the methanes, even in the 1.54 A and 3.45 A cases. In Figure 13 for the 1.54 A case, the crest-trough cancellation of the product g< )(13)g< )(23) in the region where the water molecule is in line with and outside of the two methanes is apparent. This is intuitively incorrect; it is unlikely that the presence of one methane can change the water distribution on the far side of the other methane. This failure does not occur at the large methane-methane separations because there the crest only overlaps the other g where it is equal to 1. A t the smaller R values, it becomes an important difficulty. The failure of the superposition approximation for overlapping cavities has been discussed by Pratt and Chandler (44). In building the function gl l(123) from the superposition approximation, the shape of the function between the methanes is qualitatively correct, although in some cases too large. In the outer regions, it can lead to unphysical crest-trough canceling. Because of these considerations, the treatment is modified in the following way. In the region between the methanes, the superposition approximation is applied. On the other side of the methanes, the correlation function pertaining solely to that methane is retained: 2

2

2)

2


2

2

2

2

2

2

3

(Domain 1) 2

2

g[3](123)= g< >(13)-g< >(23) for χι < χ < x

2

(19)

(Domain 2) 2

g[31(123)= g(13) for x
x (21) Two cases were considered: I) x% and x are given by the methane carbon positions, and II) they are each .75 A inside the C-C distance. It was found that placing the domain boundary directly at the carbon atoms as in method I still gives too much of the field to the superposition approximation. Therefore, in method II, the boundary is moved so that it is about .75 A (one half a van der Waals radius) inside the space between the two fragments. The result for the 1.54 A case is shown in Figure 14. In addition to a reasonable appearance of the correlation function, it is necessary that reasonable interaction energies may be calculated. In Table V , the energies for the several approximations are compared to the dynamics results. Method II gives good results for all practical distances. Even in the case of 1.54 A, which is important for building a larger molecule from fragments, the agreement is satisfactory. A t the unimportant distance of 0 A, there is an ambiguity in how a correlation function can be constructed, due to the lack of symmetry. 2

2

2



HERMANN


Figure 10. Methane-methane-water correlation function for carbon-carbon distance of 7.16 Â built from the superposition approximation.

Figure 11. Methane-methane-water correlation function for carbon-carbon distance of 6 A from the superposition approximation.


353


354


Figure 12. Methane-methane-water correlation function for carbon-carbon distance of 3.45 Â from the superposition approximation.

Figure 13. Methane-methane-water correlation function for carbon-carbon distance of 1.54 Â from the superposition approximation.



23. HERMANN


Figure 14. Methane-methane-water correlation function for carbon-carbon distance of 1.54 Â using Method II from Equations 19,20 and 21.


355

356


3

The energies in Table V were calculated from the tabular form of gt l as shown in the figures. The energy for the two methanes at infinite separation calculated in this manner agrees well with the direct dynamic results in Table I. Table V . Methane Dimer-Water Interaction Energies Using Fragment Constructed Distribution Functions 8


Carbon-Carbon Interatomic Distance 0.0 Â 1.54 3.45 5.0 6.0 7.16 8.0 oo a

Results from MD -7.070 -6.380 -5.766 -5.866 -6.081 -6.383 -6.355 -6.344

Method I

Superposition Approximation

-5.61 -6.73 -5.79 -5.73 -6.11 -6.61 -6.35 -6.37

-8.70 -7.03 -5.61 -5.80 -6.15 -6.70 -6.38 -6.37

Method II -5.61 -6.19 -5.90 -5.72 -6.10 -6.59 -6.34 -6.37

1

Distances in Angstroms. Energies in kcal m o l .

Application to Ethane. To form the function for ethane, two methane correlation functions are brought together from the other direction than in the cases above, to a distance of 1.54 A . The result of the superposition approximation produced crest trough cancellations in the wrong regions. The energy was -4.87 kcal m o l in fair agreement with the result in Table I (-5.19 kcal m o l ) . Method II produced a correlation function in good agreement with the authentic function, similar to Figure 14. The energy for Method II was -4.81 kcal m o l . 1

1

1

Discussion The methodology as developed so far has several positive attributes. It is in agreement with current calculations that show that the solvent effect favors methane association (40, 41, 42). It is also in agreement with current estimates that benzene and cyclohexane dimerize to a greater extent in the gas phase than in solution (33, 35). The treatment is not confined to aliphatic systems. Aromatic hydrocarbons may be treated along with aliphatic hydrocarbons with the same cavity potential parameters. Molecules studied here produce a decrease in surface area and usually a decrease in cavity potential on dimerization. Only in the case of Equation 10 for the high energy contact dimer of cyclohexane and the benzene dimer is the cavity potential slightly higher than the monomer pair. Thus, the free energy of association is not directly related to reduction in area or cavity potential. While the direct surface area methods will always predict a solvent effect favoring dimerization, the method of this paper can predict dimer dissociation. However, it is possible by burying enough surface area as shown in Figures 2e and 2f to get significant binding. It has been demonstrated that the equivalent sphere method of determining the cavity curvature gives results similar to the use of the explicit treatment of average mean curvature using a 2.8 Â radius sphere as a probe (12) in determining the curvature. For the simple hydrocarbon series solubilities, the inclusion of local curvature results in only a small improvement over the equivalent sphere method. In Structure and Reactivity in Aqueous Solution; Cramer, C., et al.; ACS Symposium Series; American Chemical Society: Washington, DC, 1994.

23. HERMANN


357

As a molecular complex separates into monomers, the free energy of the system goes to the correct separation limit only if Equation 10 is used. In the case of dimers, it may be seen that the hydrophobic interaction results for Equation 9 are mostly lower, and the Equation 10 results mostly higher than the results for the equivalent sphere method. The average squared mean curvature is much more dependent on the detailed shape of a surface than is the average mean curvature. In the method above of treating the curvature of cavity surfaces in which measurements of K of the accessible surface are made with a 2.8 Â spherical probe, it should be noted that the calculated value of K is sensitive to the size of the probe. Compared to potential of mean force calculations (40, 41, 42), the barrier for pulling apart two molecules from a dimer configuration to separate monomers is too high. A second minimum is present in some (40, 41, 42) but not all (45) P M F calculations. One possible reason for the higher barrier is that the accessible surface is too convoluted in the barrier region and does not represent very well the first solvation shell. In an effort to calculate ei generally, it appears that it may be possible to build distribution functions for larger systems to a good approximation from transferable solvent distribution functions of fragments, or small molecules. A method is developed here in which novel features such as build-up and interference effects appear automatically between fragments and in the vicinity of their union. Table V shows that the molecular dynamics results for ei are reasonably well reproduced by Equations 19-21. Such functions built in this manner could be incorporated into continuum solvation models. The representation of the correlation functions by analytical functions would be a simplification of the treatment. 2

m

2


m

Acknowledgments The author wishes to acknowledge discussions with Dr. James Metz and Dr. Bo Saxberg of Lilly Research Laboratories. The author also wishes to thank Professor W. Jorgensen for the use of the OPLS aliphatic parameters. Literature Cited 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16.

Cramer, C. J.; Truhlar, D. G. J. Am. Chem. Soc. 1991, 113, p. 8305. Cramer, C. J.; Truhlar, D. G. Science 1992, 256, p. 213. Hermann, R. B. J. Phys. Chem. 1972, 76, p. 2754. Hermann, R. B. J. Comp. Chem. 1993, 14, p. 741. Honig, B.; Sharp, K.; Yang, A. J. Phys. Chem. 1993, 97, p. 1101. Ooi, T.; Oobatake, M.; Némethy, G.; Scheraga, H. A. Proc. Natl. Acad. Sci. U.S.A. 1987, 84, p. 3086. Still, W. C.; Tempczyk, Α.; Hawley, R. C.; Hendrickson, T. J. Am. Chem. Soc. 1990, 112, p. 6127. Kang, Y. K.; Gibson, K. D.; Némethy, G.; Scheraga, H. A. J. Phys. Chem. 1988, 92, p. 4739. von Freyberg, B; Braun, W. J. Comp. Chem. 1993, 14, p. 510. Warshel A. J. Phys. Chem. 1979, 83, p. 1640. Wesson, L.; Eisenberg, D. Protein Science 1992, 1, p. 227. Nicholls, Α.; Sharp, Κ. Α.; Honig, B. Proteins: Struct. Funct. Genet. 1991, 11 p. 281. Lee, B.; Richards, F. M . J. Mol. Biol. 1971, 55, p. 379. Eley, D. D. Trans. Faraday Soc. 1939, 35, p. 1281. Ulig, H. H. J. Phys. Chem. 1937, 41, p. 1215. Choi, D. S.; Jhon, M. S.; Eyring, H. J. Chem. Phys. 1970, 53, p. 2608. In Structure and Reactivity in Aqueous Solution; Cramer, C., et al.; ACS Symposium Series; American Chemical Society: Washington, DC, 1994.

358

17. 18. 19. 20. 21. 22.


23. 24. 25. 26. 27. 28. 29. 30. 31. 32. 33. 34. 35. 36. 37. 38. 39. 40. 41. 42. 43. 44. 45.


Pierotti, R. A. J. Phys. Chem. 1963, 67, p. 1840. Pierotti, R. A. J. Phys. Chem. 1965, 69, p. 281. Sinanoglu, O. In Molecular Associations in Biology; Pullman, B., Ed.; Academic Press: New York, N.Y., 1968; p. 427. Hermann, R. B. J. Phys. Chem. 1975, 79, p. 163. Erratum,J.Phys. Chem. 1975, 79, p. 3080. Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M . L. J. Chem. Phys. 1983, 79, p. 926. Potsma, J. P. M.; Berendsen, J. C.; Haak, J. R. Faraday Symp. Chem. Soc. 1982, 17, p. 55. Reiss, H.; Frisch, H. L.; Lebowitz, J. L. J. Chem. Phys. 1959, 31, p. 369. Reiss, H. Adv. Chem. Phys. 1965, 9, p. 1. CRC Handbook of Chemistry and Physics; Lide, D. R., Ed.; Chemical Rubber Publishing Co.: Cleveland, OH, 1991. Jorgensen, W. L.; Severance, D.L. J. Am. Chem. Soc. 1990, 112, p. 4768. Ben Naim, Α.; Marcus, Y. J. Chem. Phys. 1984, 81, p. 2016. The Chemists Companion: A Handbook of Practical Data, Techniques and References; Gordon, A. J.; Ford, R. Α., Eds.; John Wiley & Sons: New York, N.Y., 1972. Hermann, R. B. In Seventh Jerusalem Symposium on Molecular and Quantum Pharmacology; Bergman, E; Pullman, E., Eds.; D. Reidel: Dordrecht, Holland, 1974; p. 441. DoCarmo, M. P. Differential Geometry of Curves and Surfaces; Prentice Hall: Englewood Cliffs, N.J., 1976. Burkert, U.; Allinger, N. L. Molecular Mechanics: ACS Monograph 177; American Chemical Society, Washington, D.C., 1982. Dymond, J. H.; Smith, Ε. B. Virial Coefficients of Gases: A Critical Compilation; Oxford University Press: New York, N.Y., 1980; 2nd Ed. Wood, R. H.; Thompson, P. T. Proc.Natl.Acad. Sci. U.S.A. 1990, 87, p. 946. Tucker, Ε. E.; Lane, Ε. H.; Christian, S. D. J. Solution Chem. 1981, 10, p. 1. Rossky, P. J.; Friedman, H. L. J. Phys. Chem. 1980, 84, p. 587. Hill, T. Statistical Mechanics; McGraw-Hill: New York, N.Y., 1956. Weiner, S. J.; Kolmann, P. Α.; Nguyen, D. T.; Case, D.A. J. Comp. Chem. 1986, 7, p. 230. Geiger, Α.; Rahman, Α.; Stillinger, F. H. J. Chem. Phys. 1978, 70, p. 263. Pangali, C.; Rao, M.; Berne, B. J. J. Chem. Phys. 1979, 71, p. 2982. Jorgensen, W. L.; Buckner, J. K.; Boudon, S.; Tirado-Rives, J. J. Chem. Phys. 1988, 89, p. 3742. Ravishanker, G; Mezei, M.; Beveridge, D. L. Faraday Symp. Chem. Soc. 1982, 17, p. 79. Smith, D. E.; Haymet, A. D. J. J. Chem. Phys. 1993, 98, p. 6445. Kirkwood, J. G. J. Chem. Phys. 1935, 3, p. 300. Pratt, L. R.; Chandler, D. J. Chem. Phys. 1977, 67, p. 3683. van Belle, D.; Wodak, S. J. J. Am. Chem. Soc. 1993, 115, p. 647.

RECEIVED April 25, 1994


Hydrophobic Interactions from Surface Areas, Curvature, and

Recommend Documents