Theory of Hydrophobic Bonding
163
Theory of Hydrophobic Bonding. 111. A Method for the Calculation of the Hydrophobic Interaction Based on Liquid State Perturbation Theory and a Simple Liquid Model' Robert B. Hermann The Lilly Research Laboratorieq. Eli Lilly and Company, Indianapolis, Indiana 46206 (Received March 22, 1974; Revised Manuscript Received October 7, 1974) Publication costs assisted by Eli Lilly and Company
The solubility of a pair of hydrocarbon molecules is calculated by a method incorporating first-order liquid-state perturbation theory. The method is first applied to the calculation of the solubilities of hydrocarbons in water. Water is treated as a simple liquid and only the free energy of solution at a specific temperature is calculated. It is not possible by this method to properly account for the entropy of solution in a fundamental way. A feature introduced into the theory, necessary for the calculation of hydrophobic interactions, is that the nonspherical character of the solute molecules is taken into account. The energy of the nonspherical cavity that accommodates the hydrocarbon molecule is assumed to be equal to the energy of a spherical cavity of equal area, so that liquid perturbation theory in its usual form can be applied. The hydrocarbon-solvent interaction energy is calculated with a Lennard-Jones potential between each hydrocarbon atom and using an approximate hard-sphere solvent distribution function. The fit with hydrocarbon solubility data serves to determine the one adjustable parameter in the method. The hydrophobic interaction energy between two ethane molecules is obtained for a range of distances. When the two molecules are brought together from infinite separation, a barrier is crossed at a distance where the two molecules are separated by less than the diameter of a water molecule. At a close packed distance there is about -450 cal/ mol of binding free energy at 298'. It is possible by this method to calculate the hydrophobic binding energy between two molecules in different relative orientations as well as at different distances.
Introduction In a previous paper2 the solubility of liquid hydrocarbons at a pressure of 1 atm in water was correlated with the surface area of the solvent cavity that accommodates the hydrocarbon molecule. In this paper we calculate the solubility of gaseous hydrocarbons in water and develop a method to calculate hydrophobic interactions between hydrocarbon molecules and between hydrocarbon side chains of molecules. The method is based on the first-order perturbation theory of liquid mixtures, and uses a computercalculated solvent cavity surface area to define solute size. NBmethy and Scheraga3 have developed a method for calculating hydrophobic interactions between hydrocarbon side chains, especially those of amino acid residues in polypeptides. That method was based on the assignment of a fixed quantity of free energy per water molecule in the first solvation shell around the side chains, then computing the energy associated with the exclusion of water molecules from this solvation shell when the side chains approach each other and touch. Friedman and Krishnan4 calculate the magnitude of hydrophobic interactions by a consideration of the degree of overlap of Gurney solvation cospheres about the solutes. Ben-Naim5s6has also estimated hydrophobic interactions between two hydrocarbon molecules in a novel way. This was done by comparing the solubilities of the two hydrocarbon molecules with the solubility of the molecule formed when the two molecules are moved together to the distance equal to the length of a carbon-carbon bond. This corresponds to the hydrophobic interaction energy change involved in bringing two molecules to a distance of 1.54 8, (and neglecting the two extra hydrogen atoms) between carbon atoms. While such a close ap-
proach is never observed, it does set an upper limit to the magnitude of the hydrophobic interaction. Built in to some of these methods of calculation is the idea that the solvation energy changes which are involved in the hydrophobic interaction are proportional to the size of the solvation shell. This is undoubtedly only approximately true since the energy of the water molecules in the solvation shell depends in part on the attraction between the water molecules and the total hydrocarbon molecule and thus on the density of hydrocarbon atoms present inside the solvation shell. For cyclic hydrocarbons the density of atoms in the cavity is higher than the corresponding straight chain analog, while for the hydrophobic bond itself the density of atoms would be lower. Previous estimates of hydrophobic attractions are thus too large, since calculational methods are usually based on parameterization not with aggregates of solute molecules but with single hydrocarbon solute molecules, where the atom density is higher. Within the last several years the solubility of gasses in simple liquids has been accounted for with the BarkerHenderson perturbation theory of liquid mixture^.^-^ Perturbation theory applies rigorously also to water as well as simple liquids, however, the intermolecular potential for two water molecules is complicated by the fact that it is dependent on the relative orientations of the two water molecules. In this paper we show that first-order perturbation theory can be used to calculate hydrocarbon solubilities in water using a spherically symmetric potential for water and determining the potential parameters from the surface tension of water. The equations in this form can then be used to calculate hydrophobic interaction free-energy curves for the interaction between arbitrarily oriented hydrocarbon molecules.1° The Journal of Physical Chemistry, Vol. 79, No. 2, 7975
Robert B. Hermann
164
Theory For a solute dissolved in a solvent, the Barker-Henderson perturbation theory gives, to first order, for the Helmholtz free energy A of the m i ~ t u r e : ~ A = AoHs- 4.irpikTN2diz2gi2(dl,)[diz - 6i2]
+
and the solvent cavity energy is therefore, for a spherical cavity pc = pHS + p c O R R - Pin (10) Rewriting eq 6 in terms of the cavity energy and the interaction energy16-lQ kT In
where AoHs is the Helmholtz free energy of the hard-sphere mixture reference fluid, p i is the number density of component i, k is Boltzmann's constant, T is the temperature, Ni is the number of i molecules, dij is the hard-sphere diameter between spherical molecules i and j , gzJ(r)is the twobody hard-sphere radial distribution functionl1-l3 between i and j molecules, g12(d12)is the particular value of this function a t the contact distance dlz between solvent and solute,14 u;j(r) is the spherically symmetric (perturbation) potential function between particles i and j , uu is the intermolecular distance where uij = 0. The quantities 812, dll, dzz, and dlz are found from
X,
= pc
+
pin
Qii
exp[- u,,(r)/kTI d r
(3)
= (4, + d,2)/2 (4) The chemical potential of dissolving a solute in a solvent is
4 2
Pz =
(WBNZ)V,T,N*
(5)
Neff and McQuarriel5 showed that if it assumed that the gas phase is ideal kT l n h =
XZ
pHS
+
+
pCoRR kT I n m V
(6)
where p 2 is the partial pressure of the solute at the given temperature and X 2 is the solubility in mole fractions pHs = kT{- In ( 1 - y )
+ yXR3 +
'/z [ 3 y / ( l - y)31R2 + [ 3 ~ / ( 1 Y ) I W + R ) ) (7)
Y = 7i.pid113/6
x
R = d,,/dii = (1 + y -k y 2 ) / ( 1 - y)3
which is the hard-sphere chemical potential, and pcoRR= - 47i.~ikTdi2~(di, - 6,2)gi2(di2) +
4 ~ p , k T S* u,,(r)gi2(r)r2 dr
(8)
"12
is that portion of the chemical potential due to the perturbation. It is convenient at this point to define the cavity potential. The portion of the chemical potential due to the solvent-solute interaction is, for a spherically symmetric potential pin = 4.irpikTJ * uiz(r)gi2(r)r2 dr "12
The Journal of Physical Chemistry, Vol. 79, No. 2 , 1975
(9)
kT In
Nk T V
~
ill)
In applying the above ideas to the solution of hydrocarbons and hydrocarbon pairs or multiplets in water, several assumptions are made. First, although water is not a simple liquid, we assume a Lennard-Jones potential for u 1 1 (r ) between water molecules
and we calibrate this potential, determining both parameters €11and u l l from a knowledge of the surface tension of water, for a very large cavity at a given temperature. This is done by obtaining the expression for the cavity energy y, from eq 3 and 7-10 and noticing that for a given value of d22, fiC is obtained as a function of €11 and all. The cavity energy should approach yn(d11 d22)2, the surface tension of water times the cavity area, as the cavity area goes to infinity. The value y plus this limiting condition are enough to determine the two parameters, €11 and ull. The equations above were solved numerically for several large values of d22. The only parameters found to be compatible with the above conditions were q 1 l k = 302'K and 811 = 2.67 &. for cavity diameters over the range of interest. The quantity 61 was 0.03331320and the integral involving the derivative of gll(r) with respect to Nz was solved numerically and with an approximation due to Neff and McQuarrie.15 From eq 3 dll is found to be 2.60 A. Because the angular dependence of the intermolecular potential for water has been neglected, our model cannot give a completely adequate thermodynamic treatment of hydrocarbon solutions. Our aim is the more modest one of calculating the magnitude of hydrophobic interaction free energies at a particular temperature in terms of a theory which can account for hydrocarbon solubilities in water. The treatment in this paper gives the cavity energy as a function of cavity area for a particular temperature. A more detailed treatment of water would again give (a) the cavity energy as a function of cavity area, and presumably (b) a more nearly correct temperature dependence. The more important concern of the two is that the assumption of eq 1 2 and the use of first-order perturbation theory leads to a cavity energy as a function of cavity area which is a good approximation to a more detailed treatment, for a particular temperature. To obtain the chemical potential of spherical cavity formation at a different temperature, it would be necessary to reparameterize the potential in eq 12, using the corresponding value of the surface tension for an infinite sized cavity a t the appropriate temperature. A second assumption must be made because perturbation theory in the form given can treat only spherical cavities. The hydrocarbon molecules of interest are in many cases far from spherical, and in any case the parameter u22 is not generally available. In particular, hydrophobic interactions are dependent also on the size of a pair of hydrocarbon moieties and u22 for such systems would be impossible to obtain by conventional methods. We must necessarily treat the solute particle shapes in some detail.
+
dii = uii -
+
Theory of Hydrophobic Bonding
165
In order to treat nonspherical solutes, it is therefore assumed that solvent cavities of equal areas have equal energies wc. In practice this means that we must compute the molecular cavity area A R W ( T )for a given solute and then find the correct choice of dl2 from A,w(T) = 4?7d,22
(13)
This assumption will not necessarily be true for seriously deformed spherical cavities2 but should be satisfactory for all cavities of interest here. Surface energies depend on the amount of curvature and for a given spherical cavity a deformation will both increase and decrease curvature at different points or in different directions. The integrated gaussian curvature is a topological invariant. In addition we assume d12 = 612. Since the differences between these quantities is already spall, this would seem to be a justifiable assumption in view of other approximations made. Finally, it is necessary to neglect the effect of the solution process on the internal degrees of freedom of the solute. In order to calculate the cavity energy pc associated with a given solute, it is first necessary to find the area A&T) of the solute cavity in water. This may be done by either assigning a hard-sphere diameter to the individual atoms or to the individual CH2 groups. In this paper we have chosen to calculate the cavity energy by treatifig the molecule as a collection of spherical CH2 groups. First the hard-sphere diameter is determined for methane using eq 3 and the value for ~ 2 and 2 622 for methane, which are 3.82 A and 137'K, respectively.21 The value for d(CH4) is found to be 3.65 8, and is used for all CH, groups. The value for d 11 for water above yields 3.125 b for the local cavity diameter of a CH, group. Using the area program22 described in an earlier paper2 which computes the exposed area of a collection of intersecting spheres, the area of the solute cavity is then determined as the sum of the exposed areas of the CH, groups. For the calculation of the interaction potential pin, a more general formula is used than the one given in eq 9, since a spherical potential is assumed in arriving at that equation. The integral pin becomes, with a slight change in notation
sin
4 d r dB d$
(14)
where u WM(r, 0,$) is the solvent--solute interaction potential, pw is the solvent molecule number density, and g WM(r,0,$) is an approximation to the hard-sphere distribution function about the solute M, and the integration limit uo is a point on a surface SO(r,e,$) parallel to the cavity surface and chosen to minimize the integral pin. This integral is considered to be the sum of the integrals for each constituent solute atom interacting with the solvent. Then the constituent atom-solvent molecule interaction potential u W A ( r ) is a function of r only. Letting the subscript A represent the constituent carbon and hydrogen atoms in the solute
sin $5 d r dB d+ with the following atom-water interaction potentials:
(15)
Figure 1. Graphical illustration of some quantities involved in the calculation of the solute-solvent interaction energy pin.
where CYAis the polarizability of atom A23 and PWis the dipole moment of water, 1.84 D,23 and
