Article pubs.acs.org/JPCB
Pairwise Hydrophobicity at Low Temperature: Appearance of a Stable Second Solvent-Separated Minimum with Possible Implication in Cold Denaturation Sridip Parui and Biman Jana* Department of Physical Chemistry, Indian Association for the Cultivation of Science, Jadavpur, Kolkata 700032, India S Supporting Information *
ABSTRACT: The hydrophobic effect appears to be a key driving force for many chemical and biological processes, such as protein folding, protein− protein interactions, membrane bilayer self-assembly, and so forth. In this study, we calculated the potential of mean force (PMF) using umbrella sampling technique between different model hydrophobes (methane− methane, cyclobutane−cyclobutane, and between two rodlike hydrophobes) at lower than ambient temperatures (300, 260, and 240 K). We find the appearance of a second solvent-separated minimum at ∼1.0 nm apart from the usual contact and first solvent-separated minimum in the PMF profile of the methane pair at low temperature. In the PMF between both cyclobutane and the rodlike hydrophobe pairs, the second solventseparated pair (SSSP) becomes even more stable than the first solventseparated pair (FSSP) at 240 K. Analysis of the water structure shows that, at 240 K, the core water of SSSP for the rodlike hydrophobe pair is more strongly hydrogen bonded and more tetrahedrally oriented than that of the FSSP. Strongly hydrogen-bonded ordered water molecules implicate strong water−water interactions, which are responsible for stabilization of SSSP at low temperature. This weakening of hydrophobic interactions through stabilization of SSSP may play a key role in the cold denaturation of protein.
1. INTRODUCTION In general, the unfavorable water-mediated interaction with nonpolar solutes is termed the hydrophobic effect.1−5 Its immediate consequence is the association of nonpolar solutes such that the effective water-accessible surface area is reduced. Because water is an integral part of the cellular environment,6−11 hydrophobic interactions play an important role in several biological processes, like protein folding, biomolecular recognition, and so forth.12−20 In particular, the spatial collapse of hydrophobic side-chain is believed to be the crucial step in protein folding.12−15 External perturbations (temperature, pressure, etc.) modify hydrophobic interactions,20−38 and such modifications are believed to be important for processes like cold39−41 and pressure-induced denaturation42−44 of protein. Frequently, the hydrophobic interaction between two hydrophobic solutes is determined by calculating the potential of mean force (PMF) between them, which measures the free energy cost of bringing two hydrated hydrophobic solutes from infinity to a particular distance. The characteristic features in a conventional PMF profile between two hydrophobic solutes are as follows: (i) a deep minimum when two solutes are in close contact with each other, which is referred to as the contact minimum (CM), (ii) a moderately stable second minimum that is generally referred to as the solvent-separated minimum (SSM), and corresponds to a configuration where a solute pair is separated by one layer of water molecules, and (iii) a barrier © 2017 American Chemical Society
in the PMF between the CM and SSM, which is referred to as the desolvation barrier. At room temperature, generally, CM is the global minimum in the free energy profile. Changes in the external conditions affect the relative stability of these two minima. Destabilization of CM and/or subsequent stabilization of SSM are the measure of weakening of the hydrophobic interaction. In particular, the temperature dependence of hydrophobic interactions has been studied in several contexts.23−26 Theoretically, pairwise hydrophobic interactions have been studied using integral equation (IE) theory in the liquid state. Pratt and Chandler derived the expression for the hydrophobic interaction of two spherical solutes dissolved in water using IE, and the main feature of their expression is that it does not depend on the solute−solute or solute−solvent interaction parameters.1,45 IE theory has also been used successfully to extract the enthalpic and entropic contributions to the PMF between two hydrophobic solutes.46 It has been shown that the CM is stabilized by the entropic term and SSM is stabilized by the enthalpic term using both IE and other approaches. Another fundamental yet debatable issue in the context of hydrophobic interactions is the structural ordering of water near the hydrophobic solutes.47−50 According to the famous Received: March 21, 2017 Revised: July 1, 2017 Published: July 4, 2017 7016
DOI: 10.1021/acs.jpcb.7b02676 J. Phys. Chem. B 2017, 121, 7016−7026
Article
The Journal of Physical Chemistry B “iceberg” model of Frank and Evans, water around the hydrophobes is more ordered compared to bulk water at the same conditions.51 Kauzmann et al. proposed that hydrophobic association is driven by an increase in entropy as the ordered water gets released from the surface of the hydrophobic solute to the bulk upon formation of the contact pair (CP).17 Consistent with this idea, recent simulations47,48 and experiments49,50 have supported the presence of ordered water around hydrophobic solutes, although no ordering around the hydrophobe was observed in neutron scattering experiments of methane solutions.52 Recently, the weakening of hydrophobic interactions at low temperature has received enormous attention, particularly due to its implication in the cold denaturation of protein.21,25,41 Here, we calculated the PMF between two methane molecules at various temperatures starting from 300 to 240 K. To investigate the effect of the size and shape of the hydrophobic solutes, we also calculated the PMF for another two sets of hydrophobes (cyclobutane−cyclobutane and a pair of rodlike hydrophobes). We extracted the contributions of enthalpy and entropy from the PMF profile in each case. Enthalpy and entropy contributions are discussed as a function of the structural ordering of water molecules at different configurations of the two hydrophobic solutes corresponding to those minima in the PMF. We then framed the importance of our results in the context of cold denaturation.
collect the probabilities to calculate the PMF. Sampling was performed at four different temperatures (320, 300, 260, and 240 K). Twelve windows with a space of 0.1 nm were used. Harmonic potential with a force constant k = 1000 kJ/(mol nm2) was used to harmonically restrain the distance to the desired value between the two methanes. In each window, first, a 2 ns equilibration was performed with the NPT ensemble. Production runs were carried out at NPT for another 6 ns. Nose-Hoover thermostat was used to keep the temperature constant and Parinello-Rahaman barostat was applied for pressure coupling. The weighted histogram analysis method (WHAM)61 was used to calculate the PMF. There will be an entropic contribution to the PMF due to rotation of the solutes. Entropic corrections were done by adding the term 2kBT ln(r) to the free energy obtained from WHAM.61 Uncertainty in the PMF was computed by bootstrap analysis.62 The resulting PMF was scaled to zero at large distance. 2.2. PMF between Two Cyclobutane-Type Molecules and between a Pair of Rodlike Hydrophobes. To check the shape and size effects, the PMF between two cyclobutane molecules and between two rodlike hydrophobes were also calculated. Cyclobutane and the rodlike hydrophobe are represented by four carbon beads in a square and 10 carbon beads in a linear shape, respectively, with a bond distance of 0.153 nm in the united atom description. The Lennard-Jones diameter of each particle is 0.37 nm and ε is 1.23 kJ/mol. Recently, Narayanan et al.63 computed the PMF between two rodlike infinite homopeptides using the umbrella sampling technique. The ends of the homopeptides were made to be infinite through a periodic boundary in the z-direction. The PMF calculations between the two cyclobutane-type molecules and between the two rodlike hydrophobes followed the same procedure.63 However, in our study, the solutes have a finite chain length. A harmonic force with a force constant of 2000 kJ/(mol nm2) was applied in the z-direction, and the x and y coordinates of the atoms were held fixed with a restraining force to keep the solutes parallel to each other.64,65 We performed umbrella sampling at different windows of desired distances, ranging from 0.35 to 1.5 nm. Windows were separated by a space of 0.025 nm. In each window, we performed 2 ns of equilibration and data were collected for 6 ns. 2.3. Decomposition of the PMF into Enthalpic and Entropic Contributions. The entropy of association was computed from the temperature derivative of free energy with the help of the finite differences method.66,67 Therefore, entropy could be calculated from the PMF at each solute-pair separation r (distance between center of mass of two solutes).
2. METHODS Simulations were performed in the GROMACS53 platform. We used the all-atom optimized potentials for liquid simulations (OPLS-AA) force field54,55 and TIP4P water model,56,57 which has a melting temperature of 232 K.58 It has been shown previously that the TIP4P water model56,57 is appropriate for simulation at low temperature with the OPLS-AA force field.54,55 The molecular dynamics equation of motions was integrated using 2 fs time steps using the leap-frog algorithm. Long range electrostatic interactions were calculated by the particle mesh Ewald59 method with a grid spacing of 0.16 nm and a cutoff of 1.0 nm. A cutoff of 1.0 nm was applied to evaluate van der Waals interactions. Simulations were performed at four different temperatures (320, 300, 260, and 240 K). Because our interest was to understand the hydrophobicity at low temperature, temperatures below 300 K were therefore chosen. Before execution of umbrella sampling, long 30 ns simulations were performed at each temperature at the isobaric−isothermal ensembles to equilibrate the systems. 2.1. PMF between Two Methane Molecules. The system was composed of 2 methane molecules and 1642 water molecules in a cubic box (3.67 nm × 3.67 nm × 3.67 nm). All atom representation of methane was employed. The parameter for methane is shown in Table 1. To avoid poor sampling, the umbrella sampling technique60 was used to
− S(r ) =
atomic center
σ (nm)
ε (kJ/mol)
q (e)
methane
C H carbon bead carbon bead
0.35 0.25 0.37 0.37
0.276 0.125 1.23 1.23
−0.24 0.06 0.00 0.00
cyclobutane rod
(1)
Free energy is written as a function of inter solute distance and temperature. In the present study, entropy was calculated at 300 and at 240 K and ΔT was chosen to be 20 K. The enthalpy contribution to free energy at temperature T, H(r), is obtained by
Table 1. Parameters for Methane, Cyclobutane-Type Hydrophobe, and Rodlike Hydrophobe solute
G(r , T + ΔT ) − G(r , T ) ΔT
H(r ) = G(r ) + TS(r )
(2)
The solvent contribution, Wsolv(r), to the PMF, G(r), is given by Wsolv(r ) = G(r ) − Usol − sol(r ) 7017
(3) DOI: 10.1021/acs.jpcb.7b02676 J. Phys. Chem. B 2017, 121, 7016−7026
Article
The Journal of Physical Chemistry B
Figure 1. Analysis of structure of water. (a) Radial distribution of water oxygen from carbon atom of methane at 300 K. (b) Water−water radial distribution of bulk water at 300 K. (c) Illustration of the hydrogen bond angle.
where Usol−sol(r) is the solute−solute interaction energy in vacuum.27,67−69 Because the solute−solute interaction energy is independent of temperature Ssolv(r ) = S(r )
vicinity of the solute do not have four nearest neighbors. Thus, we calculated conditional tetrahedrality (th),72 where the tagged water molecules have 2−4 neighbors. th is given by N−1
(4)
t hi = 1 −
Ssolv(r) is the entropic component of the solvent contribution.27 The enthalpy of association, H(r), has two parts: the direct solute−solute interaction in vacuum, Usol−sol(r), and the enthalpy from the solvent contribution, Hsolv(r). Hsolv(r) is given by Hsolv(r ) = H(r ) − Usol − sol(r )
N
2 ⎛ 1⎞ cos ψjik + ⎟ ⎝ 3⎠ k=j+1
∑
⎜
(7)
with 2 ≤ N ≤ 4. To calculate the hydrogen bond angle, the same cutoff for the water pair was used. Hydrogen bond angle (θ) is illustrated in Figure 1c. For each water pair, there are four possible O−O− H angles, the minimum of which is chosen. To obtain probability distributions of the O−O−O and hydrogen bond angles, a bin width of 1° was taken.73 The hydrogen bond angle distribution provides information about how strong or how weak the hydrogen bond is, and the O−O−O angle distribution is the measure of tetrahedral ordering of water.
(5)
Hsolv(r) can be divided into two components: one is the solute−solvent interaction denoted by Hsol−solv(r) and the remaining part is Hrem(r), which primarily captures the change in solvent−solvent interactions as well as the mechanical pressure−volume (PΔV) work term. Direct calculation of Hrem(r) from simulation is associated with large statistical uncertainty as fluctuation in the total solvent−solvent interactions is large. However, Hsol−solv(r) can be calculated directly from simulation. So, Hrem(r) can be calculated by subtracting Hsol−solv(r) from Hsolv(r)27,69 Hrem(r ) = Hsolv(r ) − Hsol − solv(r )
⎛9⎞ 1 ⎜ ⎟∑ N (N − 1) ⎝ 8 ⎠ j = 1
3. RESULTS AND DISCUSSION 3.1. Effect of Cooling on PMF between Two Methanes. To investigate the pairwise hydrophobic interaction at low temperature between two methane molecules, we calculated the PMF at 300, 260, and 240 K (Figure 2a). We did not choose a temperature below 240 K as the freezing temperature of the TIP4P water model is 232 K. Consistent with a previous report, we found a minimum at ∼0.4 nm, which is the CM, and another minimum at ∼0.7 nm, which is the SSM, with a desolvation barrier between them. At 300 K, the well depths of the CM (∼2.8 kJ/mol) and desolvation barrier (∼3.7 kJ/mol) are in good agreement with previous studies.23 In addition to these two minima, we also found a very shallow minimum at ∼1.0 nm at 300 K. Therefore, we called the conventional SSM (∼0.7 nm) the first solvent-separated minimum (FSSM) and the minimum at a larger distance (∼1.0 nm), the second solvent-separated minimum (SSSM). On cooling, the well depth of the CM decreases and those of the SSMs (both FSSM and SSSM) increase. Therefore, with decreasing temperature, the CP configuration is destabilized and the solvent-separated pairs (SSPs) are stabilized, relative to 300 K. These features agree with previous results using different water models.23−25 It is worthwhile noting that the second solvent-separated pair (SSSP), which is marginally stable at 300 K, becomes increasingly stable with decreasing temperature and the depth of the SSSM is comparable with that of FSSM at 240 K. The representative configurations at the minima are shown in Figure 2b. There is no water layer at the CP, however, one layer of water molecules is at the first solvent-separated pair (FSSP) and two layers of water molecules are at SSSP in between the two methane molecules. Extensive discussion on
(6)
2.4. Analysis of Water Structure. In this study, the first solvation shell (FSS) was defined by a cutoff distance between solute and solvent, here water. The cutoff distance (here 0.55 nm) was demarcated from the first minimum of the solute− water oxygen radial distribution function (Figure 1a). Here, the core water molecules of the solvent-separated pair that are between the two hydrophobic solutes were identified as follows. To detect core water molecules, two criteria were set. First, the water must be inside the FSS of the solute. A second criterion is that the minimum distance of water oxygen atom from both solutes must be within a cutoff which is the distance between two solutes. If a water molecule fulfills the first criterion, that is, it is inside the FSS of the solute, and the minimum distance from one solute is within the cutoff but from the other solute, the distance exceeds the cutoff, the water molecule will not be called core water. Using the definitions of FSS and core water, the structure of water was analyzed by calculating O−O−O angle and hydrogen bond angle distributions. For O−O−O angle calculation, the first coordination shell of central water was considered. The first coordination shell was defined as the O−O cutoff distance of 0.35 nm, which is the minimum of the O−O radial distribution of water70 (Figure 1b). We also calculated the tetrahedral order parameter. In this study, usual tetrahedrality71 cannot be used, as water molecules inside the core or in the 7018
DOI: 10.1021/acs.jpcb.7b02676 J. Phys. Chem. B 2017, 121, 7016−7026
Article
The Journal of Physical Chemistry B
Negative solvent contribution at the CM and SSM is because of the favorable entropy and enthalpy contributions from the solvent, respectively. At 240 K, a significant reduction in solvent contribution at the CM is due to a decrease in the entropic component. However, at 240 K, solvent contribution at the SSMs becomes more negative compared to that at 300 K and magnitudes of entropy and enthalpy contribution from the solvent are also increased (Figure 3d). The enthalpic component of solvent contribution arises from the methane−water interaction, water−water interaction, and mechanical pressure−volume work term (PΔV). At different methane−methane separations, we computed the interaction energy between methane and water. The value of this term at large methane−methane separation was shifted to zero. The difference in methane−water interaction energy at a particular separation distance from its value at the fully dissociated state is defined as methane−water enthalpy contribution (HM−solv). Subtraction of H M−solv from the enthalpy of solvent contribution (Hsolv) provides the remaining enthalpy of solvent contribution (Hrem), which includes mainly the water−water interaction as well as the mechanical pressure−volume work term (PΔV). Hsolv and its component contributions (HM−solv, Hrem) at 300 and 240 K are shown in Figure 3e,f, respectively. At 300 K, when methane−methane separation is less than the FSSM distance, the HM−solv term becomes positive and reaches ∼1.9 kJ/mol at the CM. This means that contact formation is disfavored by HM−solv. At the CM, 50% Hsolv comes from HM−solv. At separations greater than FSSM distance, HM−solv is zero. At 240 K, HM−solv remains nearly unaltered, indicating that cooling has little effect on HM−solv. The Hrem term, which consists of contributions primarily from changes in the water− water interaction and pressure−volume work term (PΔV), predominantly contributes to the enthalpy of solvent contribution, as shown in Figure 3e,f. Especially at separations larger than the FSSM distance, Hsolv is entirely due to Hrem. At 300 K, Hrem at the CM is positive and contributes almost 50% to the enthalpy of solvent contribution. However, at the FSSM and SSSM, Hrem becomes favorable and the respective values are ∼−2.3 and ∼−1.3 kJ/mol. Unlike HM−solv, Hrem is immensely sensitive to temperature. At 240 K, almost the entire Hrem plot shifts below the zero line, indicating that most of the contributions come from solvent−solvent organization. 3.3. Analysis of Water Structure around Methane. Figure 4a shows the hydrogen bond (H-bond) angle distribution of water in the FSS of methane compared to that of bulk water. The distribution consists of two peaks. In other words, there are two types of H-bond present in water, more linear (∼10°) and more bent (∼50°). First coordination shell (O−O distance
