Can Hydrophobic Interactions Be Correctly Reproduced by the

J. Phys. Chem. 1996, 100, 9955-9959

9955

Can Hydrophobic Interactions Be Correctly Reproduced by the Continuum Models? Jesu´ s Pitarch, Vicente Moliner,† Juan-Luis Pascual-Ahuir,* Estanislao Silla,* and In˜ aki Tun˜ o´ n Departamento de Quı´mica Fı´sica, UniVersidad de Valencia, 46100 Burjasot (Valencia), Spain ReceiVed: February 8, 1996X

The ability of the continuum models to describe hydrophobic interactions is investigated. In this work we have studied the interactions between two methane molecules in aqueous solution by means of a continuum model. The resulting potential of mean force is in good agreement with those obtained using Monte Carlo and molecular dynamics techniques. The three energy contributions appearing in the continuum energy partition (electrostatic, dispersion-repulsion, and cavitation) have been analyzed. The cavitation free energy plays the most important role of the three, determining the existence of an energy barrier between the contact minimum and the separated methane monomers. This barrier, which also appears in previous molecular simulations, can only be reproduced using the solvent-excluding surface and not with the solvent-accessible surface.

1. Introduction The hydrophobic interaction term is used to describe the tendency of nonpolar groups or molecules to aggregate in water solution.1 Hydrophobic interactions are believed to play a very important role in a variety of processes, especially in the behavior of proteins in aqueous media. The origin of these solvent-induced interactions is still unclear. In 1945 Frank and Evans2 proposed the so-called iceberg model where emphasis is made on the enhanced local structure of water around the nonpolar solute. However, computational studies and experimental advances have yielded increasing evidence against the traditional interpretation,3 and other alternative explanations, such as the reduced freedom of water molecules in the solvation shell,4 have emerged. To understand hydrophobic interactions at the microscopic level molecular simulations of nonpolar compounds in water have been carried out. The simplest model for investigating these interactions is that of two methane molecules surrounded by water molecules. Different Monte Carlo and molecular dynamics simulations of this system have led to similar results.5-11 The potential of mean force for the two methane molecules in water shows a contact minimum with a carboncarbon distance ranging from 3.7 to 4.2 Å and with an energy, with respect to the separated methanes of between 0.0 and -1.0 kcal/mol. The energy barrier to separate the two methane molecules is 0.3-1.5 kcal/mol, and the top of the barrier is placed around 4.6-5.7 Å. Computer simulations also usually predict the existence of a second solvent-separated minimum placed at a methane-methane separation between 6 and 7 Å. The nature12,13 and the existence7 of this second minimum is more controversial, and it has long been a subject of debate. The existence of the solvent-separated minimum is not essential to the hydrophobic interactions3 and quite different results have been obtained when solvent polarization is included in the simulations.7 Although molecular simulations provide valuable microscopic information on solvation problems and specifically on hydrophobic interactions, they are computationally very expensive, † Permanent address: Departamento Ciencias Experimentales, Universidad Jaume I. Castello´n, Spain. * To whom correspondence should be addressed. X Abstract published in AdVance ACS Abstracts, May 15, 1996.

S0022-3654(96)00418-2 CCC: $12.00

specially for large systems, and normally make use of oversimplified potentials. Thus, it would be interesting to investigate if other methodologies much less demanding on computing are able to reproduce the main characteristics of the hydrophobic interactions. Continuum models correctly reproduce solvation free energies of numerous compounds.14 However, to our knowledge, they have not been used to investigate hydrophobic interactions, perhaps because of the lack of information in these models of the solvent distribution around the solute and because nonelectrostatic energetic terms needs to be accurately reproduced. In order to test the ability of the continuum models to evaluate the hydrophobic interaction we have decided to carry out a study of the methane-methane interaction in water by means of a continuum model. Besides the computational advantages, the continuum model can offer a view of the hydrophobic interactions which is complementary to that given by molecular simulations; for example, in the continuum models solute and solvent polarization can easily be accounted for. Moreover, the energy partition scheme used (electrostatic, dispersion-repulsion, and cavitation free energies) can be useful to analyze the nature of the hydrophobicity problem. We will show that, when a proper solute surface model is used, the continuum model predicts the existence of a contact minimum for the methanemethane dimer which can be broken into two monomers after passing an energy barrier, in good agreement with most of the molecular simulations cited above. The continuum model used here is limited to the study of a single cavity in solution and thus does not provide information about the solvent-separated minimum where the solute cavity breaks into two. For example, no theory exits about the cavitation energy needed to create a cavity in the close vicinity of another one. This illustrates a fundamental limitation of the continuum model as opposed to Molecular Dynamics or Monte Carlo fluid simulations where the entire potential of mean force can be obtained. However, this problem could be overcome in the near future either by explicit incorporation of some solvent molecules in the calculations or by methodology advances of the continuum models. 2. Methodology To calculate the potential of mean force for two approaching methane molecules as a function of the intermolecular distance, © 1996 American Chemical Society

9956 J. Phys. Chem., Vol. 100, No. 23, 1996

Pitarch et al.

we have used the following thermodynamic cycle:12

Figure 1. Methane dimer at the contact distance.

where ∆E(R) is the interaction energy of the two methane molecules at a distance R in the gas phase, ∆Gsol(R) is the solvation free energy of the two methane molecules at a fixed distance (R), and ∆Gsol(CH4) is the solvation free energy of a methane molecule. δGsol(R) is the solvent-induced interaction term. The gas phase internal energy has been obtained at the MP2/ 6-31G** level.15 This computational level takes into account the dispersion energy between the two methane monomers and gives results in good agreement with more sophisticated calculations.16,17 Geometries have been optimized for different values of the carbon-carbon distance (R) at the same computational level using the Berny algorithm.18 The path has been followed for the relative orientation shown in Figure 1, which corresponds to the energy minimum at contact distances. At larger distances there are several orientations with very similar energy.16,17 Calculations in the gas phase have been carried out with the Gaussian92 package of programs.19 Solvation free energies are calculated as the sum of three different terms, electrostatic, dispersion-repulsion, and cavitation free energies:

∆Gsol ) ∆Gele + ∆Gd-r + ∆Gcav and thus also the solvent induced term can be partitioned into three terms:

δGsol ) δGele + δGd-r + δGcav The electrostatic free energies have been obtained by means of the PCM model20-22 at the MP2/6-31G** level using the gas phase optimized geometries and a dielectric constant of 78.4. The direct PCM-BEM version22 of the model as implemented in the HONDO program23 has been employed for these calculations. The cavity used to calculate this contribution to the solvation free energy is the solvent-excluding surface (SES).24,25 As usual, the van der Waals radii26 (Rc ) 1.7 Å and RH ) 1.2 Å) have been augmented by 20% for electrostatic calculations. The solvent sphere radius used was 1.4 Å. Dispersion-repulsion contributions have been obtained from a surface integration using the uniform distribution approximation for the solvent.27,28 This term was obtained using a 6-12 Lennard-Jones potential. The methane parameters for this potential were taken from Jorgensen et al.,29 and for water we used TIP4P dispersion-repulsion coefficients.30 Cavitation free energy was calculated with the formula proposed in ref 31 which makes use of the solute’s area of the solvent-excluding surface and of the solvent’s surface tension (104 cal/mol Å2 for water):

∆Gcav ) γS + Wo

Figure 2. Solvent-accessible surface (upper) and solvent-excluding surface (lower) of a particular solute.

TABLE 1: Free Energy Contributions to the Solvation Free Energy of a Methane Molecule in Water. Calculated and Experimental Solvation Free Energies Refer to the 1M f 1M Process (in kcal/mol) ∆Gele

∆Gd-r

∆Gcav

∆Gsol

∆Gexp

-0.19

-2.91

5.26

2.16

1.93a 2.01b

a

Reference 33. b Reference 34.

where the independent term (287 cal/mol) comes from the Scaled Particle Theory.32 In this energy partition scheme the molecular surface model chosen to represent each term is of great importance. Thus, for the dispersion-repulsion term, the solvent-accessible surface is the proper model to be used since it is obtained by means of a center-to-center potential. As can be seen in Figure 2, the solvent-accessible surface is defined by the center of the solvent sphere when rolling over the solute van der Waals surface. In the case of the electrostatic term, produced by interactions between solute and solvent charge distributions and where the electronic cloud polarization plays an important role, the solventexcluding surface (see Figure 2) is the more adequate model.22 The correct modelization of the cavitation term is still unclear. This term plays a very important role in water solutions because of the strong water-water interactions (as reflected in the large surface tension), and one goal of this paper is to shown that it can be correctly estimated using the solvent-excluding surface as proposed in reference 31. In order to provide a test of the goodness of the methodology used in this work on investigating hydrophobic interactions, we have first calculated the solvation free energy of methane in water. The different free energy contributions to the solvation free energy of a methane molecule in water at 298 K are given

Continuum Models for Hydrophobic Interactions

J. Phys. Chem., Vol. 100, No. 23, 1996 9957

Figure 3. Methane-methane potential of mean force in the gas phase.

Figure 5. Methane-methane potential of mean force in the gas phase and in aqueous solution.

Figure 4. Solvent-induced contributions to the methane-methane potential of mean force.

in Table 1. The calculated solvation free energy of methane in water compares well with the experimental estimations33,34 and with the results of molecular simulations.5 3. Results and Discussion Figure 3 shows the MP2/6-31G** energy profile as a function of the carbon-carbon distance in the gas phase. The methanemethane minimum in the gas phase is found at a carbon-carbon distance of 3.89 Å, the interaction energy relative to the separated methanes being -0.26 kcal/mol. This value is in agreement with others obtained using larger basis sets.16,17,35 The best estimation found in the literature (MP3/6-311G(3d,3p) + BSSE)17 gives an interaction energy of -0.39 kcal/mol and an equilibrium distance of about 3.8 Å. For comparison with molecular mechanics potentials used in molecular simulations, the values found using the methane-methane interaction potential of Jorgensen et al.5 are 4.2 Å and -0.29 kcal/mol. The electrostatic, dispersion-repulsion, and cavitation contributions to the potential of mean force are presented in Figure 4. Comparing the three different solvent-induced contributions some important conclusions can be made. As can be expected from its small contribution to the methane solvation, the electrostatic free energy does not produce an important contribution to the potential of mean force. In fact, in the range between 3 and 6 Å, the solvent-induced electrostatic term is nearly constant and zero. The dispersion-repulsion term is also relatively constant in the range between 3 and 6 Å although with a small positive value (about 0.5 kcal/mol). More

Figure 6. Variation of the areas of the solvent-accessible surface and solvent-excluding surface of the methane dimer as a function of the intermolecular distance.

interesting than this is the behavior of the cavitation free energy contribution. This term shows a maximum around 4.8 Å and then decreases with the intermolecular distance R. We will analyze this term more carefully later. The resulting potential of mean force is shown in Figure 5. The gas phase interaction energy is also shown for comparative purposes. The potential of mean force in solution shows a minimum placed at a carbon-carbon distance of about 3.5 Å, the net binding energy being -0.23 kcal/mol relative to the infinitely separated monomers. The maximum of the free energy curve is placed at 4.85 Å, and the free energy barrier needed to separate the methanes contact minimum in aqueous solution is 1.03 kcal/mol. These findings are in reasonable agreement with the results found by means of molecular simulations.5-11 As explained before, the investigation of the solvent-separated minimum is outside the scope of the methodology used in its present form. Returning to Figure 4 it is clear that the shape of the calculated potential of mean force in solution comes from the cavitation energy term. In turn, this is due to the behavior of the solvent-excluding surface since we have used a linear relationship between both quantities. It has been previously noted36,37 that this surface shows an interesting dependence on the intermolecular distance when a dimer dissociates. In Figure 6 the area of the more widely used solvent-accessible surface (SAS)38 is shown together with the area of the solvent-excluding surface (SES).39 The SES reaches a maximum and then diminishes up to the constant value corresponding to sum of

9958 J. Phys. Chem., Vol. 100, No. 23, 1996

Pitarch et al.

Figure 7. The solvent-excluding surface of a particular solute is the envelope of the volume excluded to the solvent considered as a whole sphere that represents its charge density.

the two separated cavities. As can be seen in Figure 4, in the continuum model partition scheme it is the cavitation term, and thus the nature of the SES, which determines the solvent-induced interaction between the two methanes. However, the SAS grows continuously up to the moment in which the cavity breaks into two (approximately at 6 Å). Of course, a parallel evolution would be observed for the van der Waals surface, because both models result from a set of spheres centered on the solute atoms. Since the cavitation term, originated by the disruption of solvent-solvent interactions, is a function which is growing with the solute’s surface area, using the solvent-accessible surface this energy term would diminish continuously from separated methanes to the contact minimum without any energy barrier. Combining the cavitation energy obtained in this way with all the other energy contributions (dispersion-repulsion and electrostatic) would result in a potential of mean force presenting a single minimum and without an energy barrier between the dimer and the separated methanes. An example of potential of mean force obtained using the solvent-accessible surface is shown below. This different behavior when using both surface models could be explained taking into account that in water solutions the electrostatic term is the main interaction between solvent molecules and consequently the larger contribution to the cavitation free energy. This electrostatic term cannot be correctly represented by interactions among only one center by each water molecule, such as the construction of the solventaccessible surface implies. The solvent-excluding surface, which gives the area of the cavity not accessible to the solvent whole sphere and which should be close to the true envelope of the volume inaccessible to the solvent charge distribution, would be a more appropriate model (see Figure 7). As seen in Figure 6, this area first increases and then diminishes when the methanes are approached. An illustration of this behavior can be found in Figure 8 in which the evolution of the solvent-excluding surface of the methane dimer in the range studied is shown. We have represented a distance for which the solvent molecule passes through the two methanes (A) and a contact distance in which a small excluded volume appears among the spheres corresponding to the atoms (C). Finally in B the distance between the carbon atoms is close to the energy barrier top, and it can be seen that the solvent-excluding surface is notably larger and consequently also the cavitation free energy. All these considerations lead us to choose the solventexcluding surface as the proper model to use to take into account these effects in continuum calculations. However, most of the continuum models use the solvent-accessible surface to parametrize the nonelectrostatic interactions.40-45 This approach offers some advantages such as the computation speed and the assignment of all the spheres to particular atoms and conse-

Figure 8. Solvent-excluding surface of the methane dimer as obtained with the GEPOL program at three different distances: (A) > 6 Å, (B) 4.75 Å, (C) 3 Å.

Figure 9. Methane-methane total potential of mean force in aqueous solution. Electrostatic and nonelectrostatic solvent-induced contributions obtained using a solvent-accessible surface based continuum model.

quently an easier evaluation of the contribution of the spheres to the solvation free energy. In many cases these models are very useful in predicting solvent effects on chemical processes.40-45 Nevertheless, the choice of the solute surface model can be decisive in reproducing some solvent-induced interactions. We have carried out the same study of methanemethane interaction in solution by means of a continuumparametrized model recently developed by some of us.45 It uses the ellipsoidal cavity46,47 to obtain the electrostatic interactions at the HF/6-31G* level while the nonelectrostatic terms are calculated by means of a parametrization based on the solventaccessible surface.45 This model has been used successfully to interpret the features of a Diels-Alder reaction in solution.45 In this paper, we have used it to reproduce the potential of mean force of the two methane molecules in aqueous solution. The total potential of mean force in solution, together with the electrostatic and nonelectrostatic solvent-induced interactions, is shown in Figure 9. As in our previous calculation, the net electrostatic contribution is close to zero but now the nonelectrostatic term has a linear dependence with the distance reflecting the behavior of the solvent-accessible surface in Figure 6. The resulting potential of mean force shows a single minimum at

Continuum Models for Hydrophobic Interactions 4.2 Å, close to previous findings, but no energy barrier is found when the methanes are approaching each other. 4. Conclusions By means of a continuum model the methane-methane interaction in aqueous solution has been studied as an example of hydrophobic interactions. In agreement with more expensive molecular simulations, we have found that the continuum model correctly predicts the existence of a contact minimum which can be reached after passing an energy barrier. In the partition energy scheme used in our model, the cavitation free energy is the most important contribution to the potential of mean force, being responsible for the energy barrier that separates the contact minimum. The electrostatic contribution to the potential of mean force for two methane molecules in water is close to zero and the dispersion-repulsion term is approximately constant in the interval studied. It should be pointed out that in our approach the cavitation free energy only depends on the solute surface area and on the solvent physical properties (surface tension and density). So, although the probe sphere radii would be the same for two solvents, differences arise from their physical properties. What makes water a unique solvent, in our continuum model, is its large surface tension due to the strong attractive water-water interactions. Thus, this cavitation term presents the same qualitative behavior for the dissociation of the methane dimer as for the dissociation of a molecule into two charged fragments (i.e. NaCl f Na+ + Cl-). What distinguishes these two problems, from the solvent point of view, is that, in the methane-methane dissociation, the net electrostatic contribution is close to zero while in the dissociation of two charged species the electrostatic contribution grows continuously (in absolute value) with the intermolecular distance. In other words, it seems that in the continuum solvation model the lack of a net solvent-induced electrostatic interaction and consequently the prevalence of the cavitation term over the other contributions is the distinctive feature of the hydrophobic interactions. Finally, we have shown that the ability of the continuum models to reproduce hydrophobic interactions critically depends on the molecular surface model selected to compute the cavitation free energy. Only the solvent-excluding surface has the correct dependence with the intermolecular distance. When using a solvent-accessible based parametrization the resulting potential of mean force does not agree with previous findings of molecular simulations. The solvent-excluding surface is able to reproduce the existence of an energy barrier and consequently is the more correct surface model for describing separation/ aggregation processes in solution and particularly hydrophobic interactions. This energy barrier could also play an important role during bond breaking and forming in many reactive processes. Acknowledgment. Calculations were carried out on a cluster of RISC/6000 of the Departamento de Quı´mica-Fı´sica (financial support of the Conselleria de Educacio´ i Ciencia de Valencia and of the Ministerio de Educacio´n y Ciencia de Espan˜a is acknowledged). The authors thank W. Dı´az for technical assistance and Ce´sar Lo´pez for helping to prepare Figure 8. J.P. acknowledges a doctoral fellowship of MEC (Spain). I.T. thanks MEC for a postdoctoral contract. Authors wish to acknowledge Prof. Tomasi for providing them with the solvation subroutines for the HONDO program. This work was supported by DGICYT Project number PB93-0699. References and Notes (1) Kauzmann, W. AdV. Protein Chem. 1959, 14, 1. (2) Frank, H. S.; Evans, M. W. J. Chem. Phys. 1945, 13, 507.

J. Phys. Chem., Vol. 100, No. 23, 1996 9959 (3) Blokzilj, W.; Engberts, J. B. F. N. Angew. Chem., Int. Ed. Engl. 1993, 32, 1545. (4) (a) Lee, B. Biopolymers 1985, 24, 813. (b) Lee, B. Biopolymers 1991, 31, 993. (5) Jorgensen, W. L.; Buckner, J. K.; Boudon, S.; Tirado-Rives, J. J. Chem. Phys. 1988, 89, 3742. (6) Pangali, C.; Rao, M.; Berne, B. J. J. Chem. Phys. 1979, 71, 2975. (7) (a) van Belle, D.; Wodak, S. J. J. Am. Chem. Soc. 1993, 115, 647. (b) New, M. H.; Berne, B. J. J. Am. Chem. Soc. 1995, 117, 7172. (8) Smith, D. E.; Zhang, L.; Haymet, A. D. J. J. Am. Chem. Soc. 1992, 114, 5875. (9) Smith, D. E.; Haymet, A. D. J. J. Chem. Phys. 1993, 98, 6445. (10) Dang, L. X. J. Chem. Phys. 1994, 100, 9032. (11) Head-Gordon, T. Chem. Phys. Lett. 1994, 227, 215. (12) Ben-Naim, A. J. Chem. Phys. 1989, 90, 7412. (13) Zichi, D. A.; Rossky, P. J. J. Chem. Phys. 1985, 83, 797. (14) Tomasi, J.; Persico, M. Chem. ReV. 1994, 94, 2027. (15) Clark, T.; Chandrasekhar, J.; Spitznagel, G. W.; Schleyer, P. v. R. J. Comput. Chem. 1983, 4, 294. (16) Tsuzuki, T.; Tanabe, K. J. Phys. Chem. 1991, 95, 2272. (17) Tsuzuki, T.; Uchimaru, T.; Tanabe, K.; Kuwajima, S. J. Phys. Chem. 1994, 98, 1830. (18) Schlegel, H. B. J. Comput. Chem. 1982, 3, 214. (19) Frisch, M. J.; Trucks, G. W.; Head-Gordon, M.; Gill, P. M. W.; Wrong, M. W.; Foresman, J. B.; Johnson, B. G.; Schlegel, H. B.; Robb, M. A.; Reprogle, E. S.; Gomperts, R.; Andres, J. L.; Raghavachari, K.; Binkley, J. S.; Gonzalez, C.; Martin, R. L.; Fox, D. J.; Defrees, D. J.; Baker, J.; Stewart, J. J. P.; Pople, J. A. Gaussian 92. Carnegie-Mellon Quantum Chemistry Publishing Unit, Pittsburgh, PA, 1992. (20) Miertus, S.; Scrocco, E.; Tomasi, J. Chem. Phys. 1981, 55, 117. (21) (a) Bonaccorsi, R.; Cimiraglia, R.; Tomasi, J. J. Comput. Chem. 1983, 4, 567. (b) Pascual-Ahuir, J. L.; Silla, E.; Bonaccorsi, R.; Tomasi, J. J. Comput. Chem. 1987, 8, 778. (22) Cammi, R.; Tomasi, J. J. Comput. Chem. 1995, 16, 1449. (23) Dupuis, M.; Farazdel, A.; Karna, S. P.; Maluenes, S. A., HONDO8, IBM Corporation, Scientific and Engineering Computations, Kingston, NY. (24) Pascual-Ahuir, J. L.; Silla, E.; Tun˜o´n, I. J. Comput. Chem. 1994, 15, 1127. (25) Pascual-Ahuir, J. L.; Silla, E.; Tun˜o´n, I. QCPE Program No. 554 (1988, 1992, 1994). (26) Bondi, A. J. Phys. Chem. 1964, 68, 441. (27) (a) Floris, F. M.; Tomasi, J. J. Comput. Chem. 1990, 10, 616. (b) Floris, F. M.; Tomasi, J.; Pascual-Ahuir, J. L. J. Comput. Chem. 1991, 12, 784. (28) Floris, F. M.; Tani, A.; Tomasi, J. Chem. Phys. 1993, 169, 11. (29) Jorgensen, W. L.; Madura, J. D.; Swenson, C. J. J. Am. Chem. Soc. 1984, 106, 6638. (30) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (31) Tun˜o´n, I.; Silla, E.; Pascual-Ahuir, J. L. Chem. Phys. Lett. 1993, 203, 289. (32) Pierotti, R. A. Chem. ReV. 1976, 76, 717. (33) Ben-Naim, A.; Marcus, Y. J. Chem. Phys. 1984, 81, 2016. (34) Abraham, M. H. J. Chem. Soc., Faraday Trans. 1 1984, 80, 153. (35) Tsuzuki, S.; Uchimaru, T.; Tanabe, K. J. Mol. Struct. (THEOCHEM) 1994, 307, 107. (36) Pascual-Ahuir, J. L.; Silla, E. J. Comput. Chem. 1990, 11, 1047. (37) Jackson, R. M.; Sternberg, M. J. E. Protein Eng. 1994, 7, 383. (38) Lee, B.; Richards, F. M. J. Mol. Biol. 1971, 55, 379. (39) Richards, F. M. Annu. ReV. Biophys. Bioeng. 1977, 6, 151. (40) Cramer, C. J.; Truhlar, D. G. J. Am. Chem. Soc. 1991, 113, 8305. (41) Cramer, C. J.; Truhlar, D. G. Science 1992, 256, 213. (42) Cramer, C. J.; Truhlar, D. G. J. Comput. Aided Mol. Des. 1992, 6, 629. (43) Luque, F. J.; Bachs, M.; Orozco, M. J. Comput. Chem. 1994, 15, 847. (44) Bachs, M.; Luque, F. J.; Orozco, M. J. Comput. Chem. 1994, 15, 446. (45) Tun˜o´n, I.; Bertra´n, J.; Ruiz-Lo´pez, M. F.; Rinaldi, D. J. Comput. Chem. 1996, 17, 148. (46) Rivail, J. L.; Rinaldi, D. Chem. Phys. 1976, 18, 233. (47) Rinaldi, D.; Rivail, J. L.; Rguini, N. J. Comput. Chem. 1992, 13, 675.

JP960418C