|Home | About | Journals | Submit | Contact Us | Français|
By means of molecular dynamics simulations of a pair of methane molecules in a TIP3P periodic water box with the NVT scheme at six temperatures and, additionally, the NPT scheme at three temperatures ranging from T = 283 to 373 K, we determined the potential of mean force (PMF) of pairs of interacting methane molecules in water as functions of distance between the methane molecules. The PMFs converge to a single baseline only for r> 11 Å at all temperatures. The curves of the dimensionless PMF obtained at different temperatures with the NVT scheme overlap almost perfectly in the region of the contact minimum and still very well in the regions of the desolvation maximum and the solvent-separated minimum, which suggests that the temperature-dependent hydrophobic interaction potentials at constant volume in united-residue force fields can be obtained by scaling the respective dimensionless potentials by RT, R being the universal gas constant. For the dimensionless potentials of mean force obtained with the NPT scheme, the depth of the contact minimum increases, whereas the height of the desolvation maximum and the depth of the solvent-separated minimum decrease with temperature, in agreement with results reported in the literature.
Hydrophobic interactions can be defined as attractive solvent-driven interactions of non-polar solute molecules in water (Ben-Naim, 1980). Pairwise hydrophobic interactions can be expressed quantitatively by a pair potential of mean force (PMF), which is characterized by a contact minimum, a solvent-separated minimum and a desolvation maximum between them (Lüdemann et al., 1996; Czaplewski et al., 2000, 2002; Shimizu and Chan, 2000; Paschek, 2004a, 2004b). Hydrophobic interactions are one of the driving forces of protein folding (Kauzmann, 1959; Dill, 1990), as well as of the formation of micelles, lipid membranes and other aggregates (Tanford, 1980). Specifically, every force field designed to study protein structure and folding must account for the hydrophobic effect, by including explicit water, implicit solvation models or implicit consideration of solvent-mediated interactions in the site–site interaction potentials (Mackerell, 2004).
The main component of hydrophobic interactions is the free energy of cavity formation and the entropy of solvent reorganization; consequently, the strength of these interactions depends on temperature. At constant pressure, aggregation is known to be more pronounced as temperature increases (Némethy and Scheraga, 1962). Several simulation studies carried out at constant volume (NVT) (Lüdemann et al., 1996) and at constant pressure (NPT) (Shimizu and Chan, 2000; Paschek, 2004a, 2004b) confirmed that the depth of the contact minimum increases with temperature, whereas the effect of temperature on the height of the desolvation barrier and the depth of the solvent-separated minimum is smaller. These researchers used diverse methodologies for PMF determination and water models, namely simulations of eight solute molecules in a periodic box with TIP3P, TIP4P, TIP5P, SPC or SPCE water and direct calculation of solute pair correlation function (Paschek et al., 2004a, 2004b), thermodynamic integration and the SPC water model (Lüdemann et al., 1996), and the particle-insertion method and the TIP4P water model (Shimizu and Chan, 2000).
Coarse-grained potentials (in which a given interaction site consists of a number of atoms), which are implemented in large-scale simulations of protein folding and dynamics (Koliński and Skolnick, 2004; Tozzini, 2007; Clementi, 2008), usually include hydrophobic interactions in the form of pairwise side chain–side chain potentials, which have the sense of potentials of mean force. The currently implemented potentials are independent of temperature, even though they are used in temperature-dependent simulations of protein folding and unfolding. Recently (Liwo et al., 2007), we introduced the temperature dependence of the multibody terms of our physics-based UNRES united-residue force field (Liwo et al., 2008); to the best of our knowledge, this is, to date, the only coarse-grained force field with temperature dependence. However, at present, UNRES includes statistical potentials of side-chain interaction determined from the PDB (Liwo et al., 1997), which have been implemented as independent of temperature. To replace the statistical potentials in UNRES with physics-based potentials, we are conducting a series of molecular simulation studies of models of pairs of amino acid side chains in water, taking into account the dependence of the PMF on orientation (Makowski, 2007a, 2007b, 2007c, 2008), which are a follow-up of our earlier more general studies of hydrophobic interactions (Czaplewski et al., 2000, 2002, 2003, 2005). The dependence of the PMF of selected pairs of amino acid side chains on their orientation in bulk water and in cylindrical hydrophobic nanopores was also studied recently by Vaitheeswaran and Thirumalai (2008).
The present study is our first attempt to introduce temperature dependence into the potentials of side-chain interactions. We have carried out an MD study of the temperature dependence of hydrophobic interactions of a pair of methane molecules at constant volume and supplementary calculations at constant pressure with a much larger periodic water box than that used in previously reported studies (Lüdemann et al., 1996; Shimizu and Chan, 2000; Paschek, 2004a, 2004b). On the basis of the results, we suggest a simple way to design temperature-dependent hydrophobic potentials of side chain–side chain interactions.
MD simulations were carried out with the AMBER suite of programs (Pearlman et al., 1995; Case et al., 2002), using the AMBER 7.0 force field at temperatures T = 283, 293, 303, 313, 323 and 373 K in the NVT ensemble (constant number of particles, volume and temperature). Constant temperature was maintained by using Berendsen's thermostat (Berendsen et al., 1984) with the coupling parameter τ = 1 ps. Additional calculations were carried out in the NPT ensemble (constant number of particles, pressure and temperature) at T = 283, 323 and 373 K and atmospheric pressure (P = 1 atm). Berendsen's barostat (Berendsen et al., 1984) with a coupling parameter τ = 0.5 ps was used to maintain constant pressure in these calculations. The integration step was 2 fs at each temperature. The bond lengths and bond angles of water molecules were maintained constant by means of the SETTLE algorithm (Miyamoto and Kollman, 1992). The pair of methane molecules was placed in a cubic periodic box containing 1983 TIP3P water molecules (Jorgensen et al., 1983). The box size in NVT calculations was 39.2 Å; consequently, the bulk water density was 0.983 g/cm3. Each methane molecule was represented by a united atom with Lennard-Jones parameters taken from Jorgensen et al. (1984), as in our previous work (Czaplewski et al., 2000, 2002, 2003, 2005). Each MD simulation consisted of a 200 ps equilibration period followed by 10 ns production calculations. A 9.0 Å cutoff was imposed on Lennard-Jones interactions, whereas particle-mesh Ewald summation with a 9.0 Å cutoff was used to compute electrostatic interactions. For each system, a series of 24 windows was run with different harmonic restraint potentials imposed on the distance between the atoms closest to the center of mass of each of the molecules, as given by Eq. (1).
where k is a force constant [k = 2 kcal/(mol Å2) was used], r the interparticle distance and ri° the center of the restraint in the ith simulation (window); the values of r° were 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0, 8.5, 9.0, 9.5, 10.0, 10.5, 11.0, 11.5, 12.0, 12.5, 13.0, 13.5, 14.0, 14.5 and 15.0 Å. A total number of 50 000 configurations were collected for each window at each temperature. The PMF was computed from the raw simulation data by using the weighted histogram analysis method (WHAM) (Kumar et al., 1992). In our earlier work (Czaplewski et al., 2005), we demonstrated that the use of restrained MD simulations with WHAM provides much better accuracy of the PMF of particles of the size of methane compared with the particle-insertion method (Widom, 1963, 1982), although the advantage of the latter method is that it provides absolute free energies, as opposed to restrained MD and WHAM. Because we used a very large water box in our study, we can safely assume that the PMF is zero at the largest distances encountered in a simulation.
where W(r; T) and are the PMF and the dimensionless PMF, respectively, at distance r and absolute temperature T, g(r; T) is the pair correlation function (determined from MD simulations) at distance r and absolute temperature T, Fcav(r; T) and are the cavity potential (Czaplewski et al., 2000, 2002) and the dimensionless cavity potential, respectively, ELJ(r) is the Lennard-Jones energy of two methane molecules at distance r and R the universal gas constant.
In Fig. 1a, the PMFs [W(r; T) of Eq. (2)] for two methane particles in water calculated based on MD simulations performed in the NVT scheme at different temperatures and in Fig. 1b, the corresponding cavity potentials [Fcav(r; T) of Eq. (4)] are shown. It can be seen that all PMFs converge to a flat baseline for r > 11 Å, which enables us to draw meaningful conclusions as to the quantitative relations between the PMFs computed at various temperatures. It can be seen that the depth of the contact minimum (for the distance about r = 3.9 Å) increases with increasing temperature, this being consistent with previous works reported in the literature (Lüdemann et al., 1996; Shimizu and Chan, 2000; Paschek, 2004a, 2004b) and in agreement with the well-known experimental observation that the strength of hydrophobic interactions increases with temperature increase (Némethy and Scheraga, 1962). The increase in the depth of the contact minimum is manifested in the cavity potentials (Fig. 1b) as a shift of the left part of the curve to the right with increasing temperature. The desolvation maximum (which occurs at a distance around r = 5 Å) and the solvent-separated minimum (at a distance around r = 7 Å) are almost independent of temperature, as seen both in the total and in the cavity PMFs (Fig. 1).
A striking feature of the dimensionless PMF in Fig. 2a is that the depth of the contact minimum is virtually independent of temperature in the range of temperatures studied. The height of the desolvation maximum and the depth of the solvent-separated minimum of both and exhibit a decreasing trend with increasing temperature (Fig. 2a and b); this is consistent with the fact that these quantities are almost insensitive to a temperature change for W(r; T) and Fcav(r; T). However, the extents of the changes are almost covered in the oscillations of and .
Our results seem to be different from that obtained by Lüdemann et al. (1996), Shimizu and Chan (2000) and Paschek (2004a) who found that the first peak in the solute–solute pair correlation function increases with increasing temperature. However, Lüdemann et al. (1996) carried out simulations in a smaller water box and in a much wider temperature range (from T = 250 to 500 K), and it should be noted that the temperature dependence of the peak of the pair correlation function of Figure 4 in Lüdemann et al. (1996) is not monotonic. The studies of Shimizu and Chan (2000) and Paschek (2004a) were carried out at constant pressure, which implies a solvent density change with temperature. It was demonstrated (Pierotti, 1965, 1976; Shimizu and Chan, 2000) that the depth of the contact minimum depends on the density of the solvent. Therefore, the results of these two studies (Shimizu and Chan, 2000; Paschek, 2004a) cannot be compared directly with ours.
It should also be noted that Shimizu and Chan (2000) used a TIP4P and not a TIP3P water model and Paschek (2004a, 2004b) used xenon and not methane as a model of a hydrophobic solute; however, the VDW diameter of xenon (3.975 Å) (Paschek, 2004a) is only slightly greater than that of methane (3.730 Å) and, from Paschek's comparative study of various water models, it follows that other water models give qualitatively similar temperature dependence of the PMF (at constant pressure) as TIP3P. Therefore, the difference between our results and those of Shimizu and Chan (2000) and Paschek (2004a) seems to be caused by the difference in ensemble type. Nevertheless, because these researchers used a different protocol from ours to determine the PMF (the particle-insertion method and direct MD simulations of eight xenon particles in a box with 500 water molecules), we carried out additional constant pressure MD simulations at T = 283, 323 and 373 K. The respective PMFs and dimensionless PMFs are shown in Fig. 3a and b, respectively. It can be seen that the depth of the contact minimum in now increases with temperature, whereas the height of the desolvation maximum and the depth of the solvent-separated minimum in W(r;T) decrease with increasing temperature, as in the studies of Shimizu and Chan (2000) and Paschek et al. (2004a, 2004b).
A practical implication of our results for designing temperature-dependent potentials of interactions of hydrophobic side chains is that, in the simplest instance, we can determine an average dimensionless PMF  at constant volume, fit it to the analytical formulas proposed in our previous work (Makowski et al., 2007a, 2007b, 2007c, 2008), and scale it by RT to compute the actual potential at a given simulation temperature. However, to justify this conclusion, more complicated hydrophobic side chains, especially with non-spherical symmetry, as considered in our recent work (Makowski et al., 2007c, 2008), and more water models should be investigated. For constant pressure calculations, Paschek (2004a) showed that, for example, the TIP4P and TIP5P water models, in which the charge distribution of the water molecule is represented better than in TIP3P, give more pronounced dependence of the first peak of the correlation function, the height of the desolvation maximum and the depth of the solvent-separated minimum on temperature. Another issue is the temperature dependence of the PMF of amphiphilic side chains, such as glutamic acid, glutamine or arginine. These studies are now being carried out in our laboratory.
A more sophisticated way to derive the temperature-dependent potentials of side-chain interaction would combine the observations that the depth of the minimum of the dimensionless PMF  and the height of the desolvation maximum of the PMF [W(r, T)] are independent of temperature. Then PMF fitting could be carried out under the conditions given by Eqs. (6) and (7).
The results of our study strongly suggest that the PMF of hydrophobic interactions can, to a good approximation, be obtained by scaling the dimensionless PMF (or the negative of the pair correlation function) by temperature, in the range of temperatures of liquid water at normal pressure. This observation, if justified by more extensive studies of more complex hydrophobic solutes and with the use of a variety of water models, could open the door to the introduction of temperature-dependent coarse-grained potentials of interactions of hydrophobic side chains. The potentials implemented at present in coarse-grained protein force fields are assumed independent of temperature, while they are used in coarse-grained simulation studies of thermodynamics and kinetics of protein folding at varying temperature. It should be noted that the assumption that the depth of the minimum of the hydrophobic side-chain interaction component of coarse-grained potentials is independent of temperature neglects its actual increase with increasing temperature. The neglect of this increase is equivalent to assuming that the strength of the simulated interactions between non-polar side chains decreases with temperature at contact distance, which does not agree with the experimental finding that increasing temperature leads to hydrophobic association. Consequently, the absence of temperature dependence of the potentials of side-chain interactions can introduce artifacts into the results of simulations and lead to false conclusions. Introducing the temperature dependence of these potentials, therefore, seems necessary.
On the other hand, based on the plots shown in Fig. 3, it can be concluded that the dimensionless PMF at constant pressure depends on temperature remarkably, in agreement with the results obtained by Paschek (2004a). Therefore, more sophisticated treatment will be required to derive coarse-grained potentials describing the interactions of hydrophobic side chains at constant pressure.
Apart from the dependence on temperature and pressure, the potentials of mean force of side-chain interactions depend strongly on orientation (Makowski et al., 2007c, 2008; Vaitheeswaran and Thirumalai, 2008). This dependence is already accounted for in the knowledge-based potentials of side-chain interactions implemented in the present UNRES (Liwo et al., 1997, 2008), and we are now developing the respective physics-based orientation-dependent potentials using molecular dynamics simulations of models of side-chain pairs in water (Makowski et al., 2007a, 2007b, 2007c, 2008) to replace the knowledge-based potentials. In view of the recent work of Vaitheeswaran and Thirumalai (2008), another point to consider is the development of potentials of side-chain interactions in confined space (corresponding to chaperones and ribosomes), because the potentials of mean force for side-chain interactions are apparently different in such environments compared with bulk water.
This work was supported by grants from the U.S. National Institutes of Health (GM-14312), the U.S. National Science Foundation (MCB05-41633) and the Polish Ministry of Science and Education (N N204 152836). M.M. was supported by a grant from the ‘Homing’ program of the Foundation for Polish Science (FNP HOM/09/2007).
This research was conducted by using the resources of (i) our 818-processor Beowulf cluster at the Baker Laboratory of Chemistry and Chemical Biology, Cornell University, (ii) the National Science Foundation Terascale Computing System at the Pittsburgh Supercomputer Center, (iii) our 45-processor Beowulf cluster at the Faculty of Chemistry, University of Gdańsk, (iv) the Informatics Center of the Metropolitan Academic Network (IC MAN) in Gdańsk and (v) the Interdisciplinary Center of Mathematical and Computer Modeling (ICM) at the University of Warsaw.
Edited by Dave Thirumalai