2910
J. Phys. Chem. B 2007, 111, 2910-2916
Simple Physics-Based Analytical Formulas for the Potentials of Mean Force for the Interaction of Amino Acid Side Chains in Water. 1. Approximate Expression for the Free Energy of Hydrophobic Association Based on a Gaussian-Overlap Model Mariusz Makowski,†,‡ Adam Liwo,† and Harold A. Scheraga*,† Baker Laboratory of Chemistry and Chemical Biology, Cornell UniVersity, Ithaca, New York 14853-1301, and Faculty of Chemistry, UniVersity of Gdan´ sk, Sobieskiego 18, 80-952 Gdan´ sk, Poland ReceiVed: September 11, 2006; In Final Form: December 13, 2006
A physics-based model is proposed to derive approximate analytical expressions for the cavity component of the free energy of hydrophobic association of spherical and spheroidal solutes in water. The model is based on the difference between the number and context of the water molecules in the hydration sphere of a hydrophobic dimer and of two isolated hydrophobic solutes. It is assumed that the water molecules touching the convex part of the molecular surface of the dimer and those in the hydration spheres of the monomers contribute equally to the free energy of solvation, and those touching the saddle part of the molecular surface of the dimer result in a more pronounced increase in free energy because of their more restricted mobility (entropy loss) and fewer favorable electrostatic interactions with other water molecules. The density of water in the hydration sphere around a single solute particle is approximated by the derivative of a Gaussian centered on the solute molecule with respect to its standard deviation. On the basis of this approximation, the number of water molecules in different parts of the hydration sphere of the dimer is expressed in terms of the first and the second mixed derivatives of the two Gaussians centered on the first and second solute molecules, respectively, with respect to the standard deviations of these Gaussians, and plausible analytical expressions for the cavity component of the hydrophobic-association energy of spherical and spheroidal solutes are introduced. As opposed to earlier hydration-shell models, our expressions reproduce the desolvation maxima in the potentials of mean force of pairs of nonpolar solutes in water, and their advantage over the models based on molecular-surface area is that they have continuous gradients in the coordinates of solute centers.
Introduction The treatment of hydration in molecular simulations is of utmost importance for the correct representation of the energetics of biological molecules. Ab initio simulations of protein-folding pathways with explicit water are still out of reach, given the present computational resources. For example, it took several months on a multiple-processor machine to carry out simulations of 1 µs of a single folding pathway of a 36-residue protein starting from an extended conformation.1 Consequently, one must seek a simplified treatment in which the solvent is represented as a continuous medium.2 Using such an approach, researchers have found that ab initio folding simulations of even up to 60-residue proteins have become practical.3-5 The free energy of hydration is usually decomposed into the electrostatic and cavity components.2,6-8 A practical approach to treat the first component is use of the generalized Born model;9 the second component is usually related to the solventexposed surface area of the molecule6,9-16 or of each atom of the molecule17 (the surface models) or to the amount of water in the hydration sphere (the hydration-shell volume models).7,8,10,18,19 These models are based on the assumption that the dominant contribution to the free energy of hydrophobic hydration arises from the maintenance of a layer of water molecules (the first hydration shell) on the surface of the cavity created by inserting a solute molecule.20 This simple model is * Corresponding author. E-mail:
[email protected]. Phone: (607) 2554034. Fax: (607) 254-4700. † Baker Laboratory of Chemistry and Chemical Biology. ‡ Faculty of Chemistry.
supported by the fact that the dependence of the free energy of hydrophobic hydration on solute radius is roughly parabolic;20 however, for small-size solutes, the entropy loss related to the volume of the created cavity, which is not incorporated explicitly in the surface-area or hydration-shell volume models, also must be considered.21 For detailed discussions of the hydrophobic effect itself and the appropriateness of various simple models, the reader is referred to the literature.20-28 The surface-area models are further divided into models on the basis of molecularsurface area (the surface drawn by the van der Waals surface of a water molecule considered as a sphere that touches a solute molecule) and the solvent-accessible surface area (SASA) (the surface drawn by the center of a water molecule touching the solute). An interesting model was developed by Fukunishi and Suzuki,29 in which one component of the free energy of hydration is proportional to SASA and the other one to the volume of the empty space between the solute molecules in which no water molecules can be accommodated. The great disadvantage of the surface-area approach is that the molecularsurface area might have a discontinuous gradient in the coordinates of the atom centers.30 The hydration-shell volume models do not suffer from this disadvantage but do not represent all the features of the variation of the free energy of solvation with the distance; in particular, they are unable to reproduce the desolvation maximum. It should be noted that, when absolute energies of hydration are not important (e.g., in simulations of protein structure and dynamics in water), the reference hydration energy can be chosen as that corresponding to interacting sites at infinite
10.1021/jp065916s CCC: $37.00 © 2007 American Chemical Society Published on Web 02/27/2007
PMF Formulas for Interaction of Side Chains in Water (1) separation; consequently, expressions for free energies of association in water and not absolute hydration energies need to be parametrized. An attractive feature of such an approach is that the association energies can, in the first approximation, be decomposed into a sum of pairwise components; such an approach was implemented in the one of Augspurger and Scheraga.19 The calculation of the cavity component of the free energy of association becomes even more difficult when coarse-grain models of polypeptide chains are used in which the interacting particles are super-atomic objects and are, generally, not spherically symmetric. For example, in our UNRES model of polypeptide chains,31-37 each side chain is represented as an ellipsoid of revolution (a spheroid). The molecular-surface area model of hydrophobic hydration would be very difficult to use for spheroidal particles; in particular, the surface area of a cluster of spheroidal particles cannot be computed analytically, as opposed to that of a cluster of spherical particles. Thus far, we incorporated the free energy of hydration into the effective potential of side-chain-side-chain interaction, which has been determined32 from Protein Data Bank (PDB) statistics, but we intend to replace this potential with a physics-based one. Alternatives to volume- or surface-based solvation models have been proposed. Levitt38 used a simple polynomial expression for hydrophobic-association energy; however, this expression is not generalizable to solutes with non-spherical symmetry. Hummer and co-workers22 introduced a model of hydrophobic solvation and association on the basis of information theory. This approach is too expensive to implement in practical computation of hydration energy, because it requires numerical evaluation of multidimensional integrals over the volume occupied by the molecule. On the basis of this model, Hummer23 proposed a simple rational expression for the cavity potential of pairs of spherical hydrophobic solutes; his expression reproduces both the desolvation maximum and solvent-separated minimum. However, his expression is a straightforward Pade´ approximation to the cavity potential determined on the basis of information theory of hydrophobic hydration. Consequently, its generalization to non-spherical solutes or inclusion of multibody terms is not straightforward. Wagoner and Baker39 proposed an approach on the basis of the scaled-particle theory (SPT)40,41 model; however, their approach involves numerical integration as does the approach of Hummer.23 In this paper, we present a simple model for the energetics of hydrophobic association, which is based on the analysis of the number and the context of the water molecules present in different parts of the hydration sphere of pairs of hydrophobic solutes. To obtain analytical expressions, we approximate the density of a single solute particle and the solvent around the solute particle as a Gaussian with dimensions and symmetry as those of the solute. Using this approach, we derive the prototype of the expression of the difference (∆Fcav) between the free energy for creation of a cavity corresponding to a pair of spherical hydrophobic particles at a given distance and that for creating two cavities corresponding to separate solute molecules. On the basis of this expression, we propose a reasonable analytical expression of practical use; we demonstrate that, with appropriate choice of parameters, the expression is able to reproduce the desolvation maximum. We also extend this approach to spheroidal particles. The derived expressions have continuous derivatives and can therefore be used in molecular simulations. 2. Theory 2.1. Pairs of Nonpolar Spherical Particles in Water. We approximate the density (in the sense of the number of non-
J. Phys. Chem. B, Vol. 111, No. 11, 2007 2911
Figure 1. Illustration of the Gaussian model of solvent density in the hydration sphere of a spherical solute. The narrower Gaussian [with standard deviation (σ)] represents the density of the solute as a function of the distance from its center and the broader Gaussian (with standard deviation σ + ∆) that of the solute + solvent from the first hydration shell. Consequently, the difference between the second and the first Gaussian (dotted line) represents the density of water in the first hydration shell.
hydrogen atoms per unit volume) of a spherical solute particle (FS) and that of the solute particle plus the water of the first hydration shell (FSw) as Gaussians centered at the center of the solute particle and their standard deviations σ and σ + ∆, respectively (eqs 1 and 2).
FS(x,y,z;xo,yo,zo,σ) ) Fog(x,y,z;xo,yo,zo,σ)
(1)
FSw(x,y,z;xo,yo,zo,σ+∆) ) Fog(x,y,z;xo,yo,zo,σ+∆)
(2)
with
g(x,y,z;xo,yo,zo,σ) )
[
exp -
]
(x - xo)2 + (y - yo)2 + (z - zo)2 2σ2
(3)
where xo, yo, and zo are the coordinates of the center of the solute particle, Fo is the density at the center of the solute particle, and g is a Gaussian. It should be noted that Fo is the same for the Gaussian representing the solute density (eq 1) and that representing the solute plus the solvent of the first solvation shell because, in both cases, the density at the center of the system corresponds to that of the solute particle. The constant σ can be taken as the radius of the solute particle, and ∆ can be identified with the thickness of the first hydration shell. It should be noted that Gaussian approximations to the density of water molecules around a solute or solvation-energy density were also used in previous volume-based hydration models.7,8,18,19 The number density of water, N, in the first hydration shell of a spherical solute particle can now be expressed in terms of the difference in the density of the solute plus solvent (eq 2) and the solute systems (eq 1) by eq 4, i.e., as the difference between two Gaussians.
N(x,y,z) ) FSw(x,y,z;xo,yo,zo,σ+∆) - FS(x,y,z;xo,yo,zo,σ) ) Fo[g(x,y,z;xo,yo,zo,σ+∆) - g(x,y,z;xo,yo,zo,σ)] (4) A plot of the resulting density and the component Gaussians along the radius of the site is shown in Figure 1.
2912 J. Phys. Chem. B, Vol. 111, No. 11, 2007
Makowski et al.
Figure 3. Partition of the hydration sphere of a pair of spherical solute particles implemented in the derivation of eq 6. See text for details. Figure 2. Solid line: plot of the difference of two Gaussians with σ ) 1.9 Å + 1.4 Å ) 3.3 Å (corresponding to the sum of the van der Waals radii of the methane and the water molecules) and σ ) 1.9 Å (corresponding to the van der Waals radius of the methane molecule) representing the water density around a spherical solute particle (eq 4). Dashed line: the differential of the Gaussian with σ ) 1.9 Å and ∆ ) 1.4 Å (eq 5). Dotted line: the differential of the Gaussian with σ' ) 2.45 Å and ∆′ ) 1.29 Å; this differential fits the difference of two Gaussians (drawn as a solid line) most closely. The coordinate x is a coordinate along the diameter of the solute particle with origin in the center of the particle.
For easier handling of the equations in further considerations, we approximate the difference of two Gaussians in eq 4 by a differential expression:
N(x,y,z) ≈ No
∂g(x,y,z;xo,yo,zo,σ) ∆ ∂σ
(5)
Equation 5 is only a qualitative approximation when considering the difference of the radius of a water molecule and a typical nonpolar solute molecule. For example, for a methane molecule, the mean van der Waals radius can be taken as 1.9 Å and that of a water molecule as 1.4 Å.22 However, we note that eq 5 still captures the qualitative feature of the parent difference expression (eq 4) even for ∆ comparable to σ, as illustrated in Figure 2; moreover, Figure 2 shows that values of σ′ and ∆′ exist such that the curve computed from eq 5 with σ′ and ∆′ follows very closely the curve computed from eq 4 with σ and ∆. Moreover, the peak of the distribution of the solvent around the solute resembles the peak of the carbon-oxygen correlation function determined experimentally for the methane-water system; the first peak comes from the water molecules in the first hydration sphere of methane.40 It should be noted that No∆ should be approximately proportional to σ2, because the number of water molecules in the first hydration shell is proportional to the surface area of the shell [4π(σ + ∆)2]. As in the model of hydrophobic hydration and hydrophobic interaction of Ne´methy and Scheraga,24,25 and the analytical hydration-shell models developed later,10,13-15,42,43 we assume that the free energy of hydration of a nonpolar solute arises from the change of the pattern of interactions of the water molecules of the first hydration shell. The water molecules that are in contact with the nonpolar solute lose favorable electrostatic interactions with the neighboring water molecules as well as part of the configurational entropy, because orientation with protons pointing opposite to the solute surface are forced to maximize hydrogen-bonding interactions with the water molecules from the second hydration shell. Let us partition the first hydration shell of a nonpolar-solute dimer in water, as in Figure 3, into the exposed parts (A), each coming from a different monomer, the part corresponding to double overlap of the hydration shells of the monomers (C) excluding the cavity (D) created between the monomers, which
arises because of steric repulsion between the solvent molecules and between the solvent and solute molecules. Part B, also shown in Figure 3, corresponds to water expelled from the hydration shells of the monomers because of the overlap with the neighboring solute particle. As in the hydration-shell volume models,7-8,10,18,19 we assume that, for isolated spherical monomeric solutes, the free energy of hydration can be assumed to be proportional to the number of water molecules in its hydration shell, because each such water molecule experiences the same average surrounding. Because the water molecules in part A of the hydration shell of the dimer in Figure 3 are in the same immediate surrounding as those in the hydration shell of a monomer, we also assume that the part of the hydration energy of the dimer corresponding to that part of the hydration shell is proportional to the number of water molecules in A with the same coefficient as for the monomers. By contrast, the water molecules in part C of the hydration shell touch a saddle surface of the cavity between the solute molecules and not a convex surface and, consequently, have fewer neighboring water molecules compared to water in part A. Therefore, they lose more favorable electrostatic and nonbonded interactions with other water molecules. Moreover, the mobility of the water molecules in part C is more restricted than that of the water molecules in part A, which causes an increased entropy loss. The second observation is supported by the results of our recent molecular-simulation study,44 in which we demonstrated that the water molecules of the first hydration shell of a pair of methane molecules in water that are in the region between the two methane molecules are more ordered (and therefore, their mobility is more restricted) compared with those in the other part of the hydration shell. Consequently, a water molecule in part C has a more positive contribution to the free energy of hydration than that in part A. The presence of a cavity in part D also increases the free energy of the solute-solvent system.45 Thus, the cavity part of the free energy of formation of the dimer (∆Fcav) can be expressed as follows: (2) ∆Fcav ) F(1,2) - F(1) - F(2) ) η[NA + wNC - N(1) s - Ns ] ) (2) (1) (2) η[N(1) s + Ns - NB + (w - 2)NC - Ns - Ns ] ) -ηNB + η[(w - 2)NC + (w′ - 2)ND] (6)
because (see Figure 3) (2) NA ) N(1) s + Ns - NB - 2NC - 2ND
(7)
where F(1,2) is the cavity part of the free energy of hydration of the dimer, F(1) and F(2) are the cavity contributions to the free energies of hydration of monomers 1 and 2, respectively, (2) are the numbers of water molecules of the N(1) s and Ns hydration shells of the isolated first and the second monomers, respectively, NA and NC are the numbers of water molecules in regions A and C, NB and ND are the numbers of water molecules expelled from regions B and D of the hydration shells of the
PMF Formulas for Interaction of Side Chains in Water (1) isolated monomers, and η is a coefficient relating the number of water molecules in the hydration shell of a spherical solute particle to its free energy of hydration (and also, by the assumption about the water molecules in part C, as stated above eq 6, the part of the free energy of hydration of the dimer corresponding to part A of the dimer hydration shell, and consequently to the number of water molecules in part A); w > 1 is the ratio of the contribution of a water molecule in part C to the hydration free energy to that of a water molecule in the hydration shell of an isolated monomer, and w′ > 1 is the ratio of the free-energy contribution of a cavity (left in place of a water molecule that occupied a part of the hydration shell of a single solute particle) to that of a molecule in the hydration shell of a single solute particle in water; it should be noted that this part of the free energy of hydrophobic association has the same sense as the void-volume contribution (Evoid) to the free energy of hydration in the model of Fukunishi and Suzuki.28 Using the Gaussian approximation of the number density of water molecules given by eq 4 and of the density of the solute molecules (eq 1), we can express the number of water molecules of part B of Figure 3 as a sum of overlap integrals of N1 (the number density of the water molecules in the hydration shell of the isolated first particle) with g2 (the density of the second particle) and N2 (the number density of the water molecules in the hydration shell of the second particle) with g1 (the density of the first isolated particle). We further assume that (w - 2)NC + (w′ - 2)ND > 0 and (w - 2)NC + (w′ - 2)ND is proportional to the overlap integral between the water densities of the unperturbed hydration shells of the first and the second solute molecules (eq 4) over the range s1 + s2 < d < s1 + s2 + 2∆ (i.e., from the closest-approach distance of the solute particles to the largest distance at which their hydration spheres overlap). Without loss of generality, we can locate the center of the first solute molecule at the origin of the coordinate system and that of the second solute particle on the x-axis at a distance d from the origin. Thus, ∆Fcav can be approximated by eq 8.
∫∫∫ N1g2 dx dy dz + ∫∫∫ N2g1 dx dy dz + β ∫∫∫ N1N2 dx dy dz)
∆Fcav(d) = -R(
With the use of eq 5, this becomes
(
∆Fcav(d) ) -R ∆1 ∆2
∂g
∂g
∫∫∫ ∂σ11 g2 dx dy dz +
∫∫∫ g1 ∂σ22 dx dy dz
) (
) -R ∆1
)
∂G12 ∂G12 + ∆2 + ∂σ1 ∂σ2 2
β∆1∆2
∂ G12 (8) ∂σ1∂σ2
where R and β are coefficients introduced to encompass η, w, and w′ of eq 6 as well as No of eq 4; these coefficients can be regarded as the weight of the contribution due to region B of Figure 3 (which comes from expulsion of the water from the solvation spheres of isolated solutes) and the cumulative contribution due to regions C and D (which can be attributed to creating an additional cavity corresponding to region D and to the expulsion of part of water from the solvation spheres and change of the state of the water molecules remaining in this part of the solvation spheres, as explained earlier in this
J. Phys. Chem. B, Vol. 111, No. 11, 2007 2913 section), respectively, and
G12 ) G12(d;σ1,σ2) )
∫∫∫ g1g2 dx dy dz )
( ) [ 3
1d2 1 + 2 2 exp 2 2 σ1 σ2 2(σ1 + σ22)
]
N1 ) N(x,y,z;0,0,0,σ1,∆1) N2 ) N(x,y,z;d,0,0,σ2,∆2) g1 ) g(x,y,z;0,0,0,σ1) g2 ) g(x,y,z;d,0,0,σ2)
(9)
where the integration in eqs 8 and 9 is carried out over the entire space. Finally, the approximate expression for ∆Fcav is given by eq 10.
( )
∆Fcav ) [R(1) x4 + R(2) x2 - R(3)] exp -
x2 2
(10)
with
x)
R(1) )
R(2) )
(
d
(11)
xσ21 + σ22
(
)
5 1 - 2 βδ1δ2 1 + σ21 σ22 σ21 + σ22
) { [( 5
)
(12)
]
1 1 11 10 + 2 2 βδ1δ2 3 2 + 2 - 2 2 σ1 σ2 σ1 σ2 σ1 + σ22 R (δ + δ2) σ1σ2 1
R(3) )
( ) [( )
}
(13)
]
5 δ1δ2 1 -2 1 1 1 + 3 2 + 2 (δ1 + δ2) - 15 2 2 2 σ1 σ2 σ1 σ2 σ1 + σ22 (14) ∆1 ∆2 δ1 ) δ2 ) (15) σ1 σ2
Sample graphs of eq 10 for a pair of solute particles with the same size and for different value of the coefficient β are shown in Figure 4a. It can be seen that, for small β, the desolvation maximum is not reproduced and, for too large values of β, ∆Fcav becomes greater at distances less that σ1 + σ2, which is unphysical. There is, however, a range of β for which the graph is qualitatively similar to that of the cavity component of the potential of mean force of a pair of nonpolar solutes determined by molecular simulations43,44,46 with a well-developed desolvation maximum; the best value of β appears to be 1.0; β ) 1.5 (short-dashed lines in Figure 4) results in too high values of the approximate ∆Fcav on the left of the maximum. A comparison of the ∆Fcav curve drawn based on eq 10 with the ratio β/R ) 1.0 with that of the cavity contribution to the PMF based on the simulation data of the methane-dimer PMF in water of ref 49 is presented in Figure 4b; both R and β of Figure 4a have been scaled by 4.5 to fit the height of the maximum of ∆Fcav of eq 10 to that of the cavity contribution to the PMF obtained by simulations. It can be seen that the shape of the simulated PMF is qualitatively reproduced by ∆Fcav of eq 10 except that it does not reproduce the second solvent-separated
2914 J. Phys. Chem. B, Vol. 111, No. 11, 2007
Makowski et al. molecular-surface area and the cavity component of the potential of mean force with distance. 1
∆Fcav )
Figure 4. (a) Variation of the Gaussian-overlap-based approximate expression for ∆Fcav of two spherical solute particles (eq 10) with the distance between the centers of the particles for R ) 1 (where R is the weight of the contribution from region B of Figure 3) and different values of β (the cumulative weight of the contribution from regions C and D of Figure 3); R and β appear in eq 8. The values σ1 ) σ2) 1.5 Å and δ1 ) δ2 ) 1 (eq 13) were assumed. Solid line, β ) 0.5; longdashed line, β ) 1.0; short-dashed line, β ) 1.5; dotted line, β ) 2.0. The curve obtained for β ) 1.0 resembles the dependence of the cavity component of the PMF (∆Fcav; see, e.g., Figure 6 in ref 36) and that of the molecular-surface area31 on the distance between the centers of the two solute particles. (b) Superposition of the plot of the cavity contribution to the PMF of hydrophobic association of the methane dimer in TIP3P water drawn using the simulation data of ref 41 (solid line) and the simplified expression for ∆Fcav (eq 10) with R ) 4.5 and β ) 4.5, the other parameters of eq 10 being those in part a (longdashed line). The values of R and β have been chosen to adjust the height of the maximum of the simplified expression to that of the simulated PMF; it should be noted that, except for the scaling factor of 4.5 for a and b, the long-dashed-line plot in part b of the figure corresponds to the long-dashed-line plot in part a.
minimum (which is not handled by our Gaussian-overlap model). Also, the curve of eq 10 is only less steep than that of the cavity contribution to the simulated PMF. The maximum in the expression given by eq 10 arises from the term corresponding to the overlap of the hydration shells. The appearance of this maximum is a direct consequence of the fact that the water present in part C of the hydration shell (Figure 3), whose amount is approximated by the overlap integral of the water density in the hydration shells of isolated monomers, contributes a net positive free energy (eq 6); this is the crucial part of our model. It should be noted that the graph computed from eq 10 has a minimum at x ) 0, which is not a feature of the cavity component of the PMF. We note that, for d < σ1 + σ2, water molecules from part C start touching the convex part of the molecular surface, and consequently, the factor w in eq 7 decreases until reaching 1 for d ) 0 (i.e., complete overlap of the solute particles). This is not an issue though because the repulsive part of the effective site-site potentials will prevent the particles from reaching small distances. Nevertheless, we rectified this undesirable feature by replacing the parent expression (eq 10) with one in which the lowest power of x is 1/2 and not 2 and keeping the ratio of the powers of x at 1:2 as in eq 10; such a modification reproduces the steep ascent of the cavity contribution to the PMF of a pair of hydrophobic solutes in water at short distances compared to that arising from applying eq 10 directly (Figure 4b). We also replaced the exponential term of eq 10, which is expensive to evaluate, with a rational function decreasing with the distance. After a number of trials, we found that the expression given by eq 16 gives the best resemblance to the variation of the
R(1)x2 + R(2)x - R(3)
(16)
1 + R(4) x12
Interestingly, as shown in the second paper in this series,46 this equation can also reproduce the solvent-separated minimum, although it was not part of our model. It should be noted that models based on the Gaussian representation of the hydration-energy density have already been developed and implemented in molecular simulations.7,8,18,19 However, these models partition the hydration shell (and, consequently, the coefficient that relates hydration free energy to solvent density, η in eq 6) only according to the kind of the atom to which the respective part of the hydration shell belongs and not according to the kind of surface to which it is bordered on the solute side (Figure 3). This omission prevents the earlier Gaussian-overlap-based models from reproducing the desolvation maximum in the cavity potential of hydrophobic association. Our treatment emphasizes the free-energy increase due to the presence of water molecules near the saddle or concave part of the molecular surface of a hydrophobic cluster (part C shown in Figure 3), which is consistent with the model of hydrophobic hydration and hydrophobic interaction developed by Ne´methy and Scheraga.24,25 This feature enables us to reproduce the desolvation maximum in the cavity potential. We also note that all hydration-shell models,10,13-16,42,43 except the excluded hydration-shell molecular volume (ESSMV) model proposed in our earlier work,43 are unable to reproduce the desolvation maximum for exactly the same reason as the earlier Gaussianoverlap models.7,8,18,19 The ESSMV model43 reproduces the desolvation maximum because the hydration shell of a hydrophobic cluster is redefined as part of the space around the cluster formed by a water molecule rolled over its solvent-accessible surface and not the van der Waals surface. However, such a modification is not justified by the molecular origin of hydrophobic hydration. Because eq 16 approximates the cavity component of the relative hydrophobic-association free energy, the cavity component of the relative hydrophobic association of the whole system in the pairwise approximation is expressed by
∆Fcav )
∆Fcav;i,j ∑ i 2 sites in the standard deviations of all Gaussians. Such an approach looks promising with regard to the introduction of multibody terms in the side-chain interaction potential of UNRES; however, it should be noted that parametrization of these multibody terms would require much greater effort than the parametrization of the effective pairwise potentials of side-chain interaction, because it requires the determination of the potentials of mean force of the association of amino acid side
2916 J. Phys. Chem. B, Vol. 111, No. 11, 2007 chains of all kinds into clusters composed of at least three side chains. We therefore leave the derivation and parametrization of such multibody terms to our future work. Acknowledgment. This work was supported by grants from the U.S. National Institutes of Health (GM-14312), the U.S. National Science Foundation (MCB05-41633), the National Institutes of Health Fogarty International Center Grant TW007193, and the Polish Ministry of Science and Informatization (1 T09A 099 30). Mariusz Makowski was supported by a grant from the “Homing” program of the Foundation for Polish Science FNP. This research was conducted by using the resources of (a) our 818-processor Beowulf cluster at the Baker Laboratory of Chemistry and Chemical Biology, Cornell University, (b) the National Science Foundation Terascale Computing System at the Pittsburgh Supercomputer Center, (c) our 45-processor Beowulf cluster at the Faculty of Chemistry, University of Gdan´sk, (d) the Informatics Center of the Metropolitan Academic Network (IC MAN) in Gdan´sk, and (e) the Interdisciplinary Center of Mathematical and Computer Modeling (ICM) at the University of Warsaw. References and Notes (1) Duan, Y.; Kollman, P. A. Science 1998, 282, 740. (2) Cramer, C. J.; Truhlar, D. G. Chem. ReV. 1999, 99, 2161. (3) Vila, J. A.; Ripoll, D. R.; Scheraga, H. A. Proc. Natl. Acad. Sci. U.S.A. 2003, 100, 14812. (4) Jang, H. B.; Hall, C. K.; Zhou, Y. Q. Biophys. J. 2004, 86, 31. (5) Schug, A.; Wenzel, W. Biophys. J. 2006, 90, 4273. (6) Fraternali, F.; van Gunsteren, W. F. J. Mol. Biol. 1996, 256, 939. (7) Lazaridis, T.; Karplus, M. Proteins: Struct. Func. Genet. 1999, 35, 133. (8) Wagner, F.; Simonson, T. J. Comput. Chem. 1999, 20, 322. (9) Tomasi, J.; Mennucci, B.; Cammi, R. Chem. ReV. 2005, 105, 2999. (10) Gibson, K. D.; Scheraga, H. A. Proc. Natl. Acad. Sci. U.S.A. 1967, 58, 420. (11) Hodes, Z. I.; Ne´methy, G.; Scheraga, H. A. Biopolymers 1979, 18, 1565. (12) Hodes, Z. I.; Ne´methy, G.; Scheraga, H. A. Biopolymers 1979, 18, 1611. (13) Kang, Y. K.; Ne´methy, G.; Scheraga, H. A. J. Chem. Phys. 1987, 91, 4105. (14) Kang, Y. K.; Ne´methy, G.; Scheraga, H. A. J. Chem. Phys. 1987, 91, 4109. (15) Kang, Y. K.; Ne´methy, G.; Scheraga, H. A. J. Chem. Phys. 1987, 91, 4118. (16) Kang, Y. K.; Ne´methy, G.; Scheraga, H. A. J. Phys. Chem. 1988, 92, 4739. (17) Vila, J.; Williams, R. L.; Vasquez, M.; Scheraga, H. A. Proteins: Struct. Funct. Genet. 1991, 10, 199.
Makowski et al. (18) Stouten, P. F. W.; Fro¨mmel, C.; Nakamura, H.; Sander, C. Mol. Simul. 1993, 102, 97. (19) Augspurger, J. D.; Scheraga, H. A. J. Comput. Chem. 1996, 13, 1549. (20) Ho¨finger, S.; Zerbetto, F. Chem. Soc. ReV. 2005, 34, 1012. (21) Graziano, G. J. Chem. Soc., Faraday Trans. 1998, 94, 3345. (22) Hummer, G.; Garde, S.; Garcı´a, A. E.; Pohorille, A.; Pratt, L. R. Proc. Natl. Acad. Sci. U.S.A. 1996, 93, 8951. (23) Hummer, G. J. Am. Chem. Soc. 1999, 121, 6299. (24) Ne´methy, G.; Scheraga, H. A. J. Chem. Phys. 1962, 36, 3401. (25) Ne´methy, G.; Scheraga, H. A. J. Phys. Chem. 1962, 66, 1773. (26) Lum, K.; Chandler, D.; Weeks, J. D. J. Phys. Chem. B 1999, 103, 4570. (27) Southall, N. T.; Dill, K. A.; Haymet, D. J. J. Phys. Chem. B 2002, 106, 521. (28) Paschek, D. J. Chem. Phys. 2004, 120, 6674. (29) Fukunishi, Y.; Suzuki, M. J. Comput. Chem. 1997, 18, 1656. (30) Wawak, R. J.; Gibson, K. D.; Scheraga, H. A. J. Math. Chem. 1994, 15, 207. (31) Liwo, A.; Pincus, M. R.; Wawak, R. J.; Rackovsky, S.; Scheraga, H. A. Protein Sci. 1993, 2, 1697. (32) Liwo, A.; Ołdziej, S.; Pincus, M. R.; Wawak, R. J.; Rackovsky, S.; Scheraga, H. A. J. Comput. Chem. 1997, 18, 849. (33) Liwo, A.; Czaplewski, C.; Pillardy, J.; Scheraga H. A. J. Chem. Phys. 2001, 115, 2324. (34) Liwo, A.; Ołdziej, S.; Czaplewski, C.; Kozłowska, U.; Scheraga, H. A. J. Phys. Chem. B 2004, 108, 9421. (35) Czaplewski, C.; Liwo, A.; Ołdziej, S.; Scheraga, H. A. Polymer 2004, 45, 677. (36) Ołdziej, S.; Lagiewka, J.; Liwo, A.; Czaplewski, C.; Chinchio, M.; Nanias, M.; Scheraga, H. A. J. Phys. Chem. B 2004, 108, 16950. (37) Scheraga, H. A.; Liwo, A.; Ołdziej, S.; Czaplewski, C.; Pillardy, J.; Ripoll, D. R.; Vila, J. A.; Kaz´mierkiewicz, R.; Saunders, J. A.; Arnautova, Y. A.; Jagielska, A.; Chinchio, M.; Nanias, M. Front. Biosci. 2004, 9, 3296. (38) Levitt, M. J. Mol. Biol. 1976, 104, 59. (39) Wagoner, J. A.; Baker, N. A. Proc. Natl. Acad. Sci. U.S.A. 2006, 103, 8331. (40) De Jong, P. H. K.; Wilson, J. E.; Neilson, G. W.; Buckingham, A. D. Mol. Phys. 1997, 91, 99. (41) Czaplewski, C.; Kalinowski, S.; Liwo, A.; Scheraga, H. A. Mol. Phys. 2005, 3157. (42) Perrot, G.; Cheng, B.; Gibson, K. D.; Vila, J.; Palmer, K. A.; Nayeem, A.; Maigret, B.; Scheraga, H. A. J. Comput. Chem. 1992, 13, 1. (43) Czaplewski, C.; Ripoll, D. R.; Liwo, A.; Rodziewicz-Motowidło, S.; Wawak, R. J.; Scheraga, H. A. Int. J. Quantum Chem. 2002, 88, 41. (44) Czaplewski, C.; Liwo, A.; Ripoll, D. R.; Scheraga, H. A. J. Phys. Chem. B 2005, 109, 8108. (45) Lee, B. J. Chem. Phys. 1985, 83, 2421. (46) Makowski, M.; Liwo, A.; Maksimiak, K.; Makowska, J.; Scheraga, H. A. J. Phys. Chem. B 2007, 111, 2917. (47) Berne, B. J.; Pechukas, P. J. Chem. Phys. 1972, 56, 4213. (48) Gay, J. G.; Berne, B. J. J. Chem. Phys. 1981, 74, 3316. (49) Makowski, M.; Liwo, A.; Sobolewski, E.; Czaplewski, C.; Ołdziej, S.; Scheraga, H. A. J. Phys. Chem. B 2007, 111, 2925.
