Revisiting Hydrogen Bond Thermodynamics in Molecular Simulations

May 10, 2017 - The bond that can form between a hydroxyl group and an oxygen atom can be quantified by the interoxygen distance r and the angle θ. B)...
1 downloads 9 Views 1MB Size
Download PDF
Article pubs.acs.org/JCTC

Revisiting Hydrogen Bond Thermodynamics in Molecular Simulations Liel Sapir† and Daniel Harries* Institute of Chemistry and The Fritz Haber Research Center, The Hebrew University, Jerusalem 91904, Israel S Supporting Information *

ABSTRACT: In processes involving aqueous solutions and in almost every biomolecular interaction, hydrogen bonds play important roles. Though weak compared to the covalent bond, hydrogen bonds modify the stability and conformation of numerous small and large molecules and modulate their intermolecular interactions. We propose a simple methodology for extracting hydrogen bond strength from atomistic level simulations. The free energy associated with hydrogen bond formation is conveniently calculated as the reversible work required to reshape a completely random pair probability distribution reference state into the one found in simulations where hydrogen bonds are formed. Requiring only the probability density distribution of donor−acceptor pairs in the first solvation shell of an electronegative atom, the method uniquely defines the free energy, entropy, and enthalpy of the hydrogen bond. The method can be easily extended to molecules other than water and to multiple component mixtures. We demonstrate and apply this methodology to hydrogen bonds that form in molecular dynamics simulations between water molecules in pure water, as well as to bonds formed between different molecules in a binary mixture of a sugar (trehalose) and water. Finally, we comment on how the method should be useful in assessing the role of hydrogen bonds in different molecular mechanisms.

1. INTRODUCTION From water’s anomalous properties to biomacromolecular stability in aqueous solutions, hydrogen bonds (or “Hbonds”) are pivotal.1 Hbonds are central, for example, to the emergence of the hydrophobic interaction responsible for multiple processes in aqueous solutions,2 thus helping proteins to fold and micelles to assemble. It is often important to determine quantitatively not only the free energy of Hbond formation but also the associated entropic and enthalpic components. Simple and consistent methods to measure Hbond thermodynamic properties are therefore important because they can help to determine details of molecular-level solvation mechanisms. Numerous experimental3−9 and computational10−12 studies have approached the problem of evaluating the strength of Hbonds, with many focusing specifically on water−water interactions. The derived values of the different thermodynamic quantities from the various methods spread over a wide range (see refs 13 and 14). One challenge is the lack of a clear and unambiguous “unbound” reference state with which to compare the Hbond. In many spectroscopic studies, for example, a two-state model is assumed, in which an Hbond can be either intact or “broken”.13,15,16 The equilibrium between these states is said to represent the Hbond formation “reaction”. In other studies, Hbonds in ice are used as a reference, and the distortion of Hbonds upon heating is quantified from the changes in the angular distribution.17−19 In atomistic level simulations an additional complication often arises because empirical force fields do not uniquely define a priori the Hbond interaction. Instead, the Hbond in these © 2017 American Chemical Society

simulations is an assignable property that emerges from all calibrated forces, including electrostatic and van der Waals interactions. Thus, several approaches have been developed to define Hbonds in simulations based on empirical force fields and, in certain cases, to also quantify their strength. Motivated by the fact that Hbonds can be inferred from the respective nonrandom probability distributions of mutual positions and angles of atoms participating in the bond, many proposed methods use geometric,11,20−25 or a combination of geometric and energetic,10,13 strict criteria to identify and enumerate Hbonds. Several strategies have been proposed to ascertain Hbond strength from computer simulations. One strategy is to first enumerate Hbonds, as defined above, and then deduce the free energy of bond formation from the corresponding distribution of the number of Hbonds that form with each electronegative atom in solution.26 A similar approach uses the potential of mean force describing the nearing of donor to acceptor.20 Other approaches use the underlying force field or ab initio derived energy terms.27−30 An alternative strategy is based on empirical calibration from experimental evidence of the hydrogen bond energy with respect to one or more geometric properties of the Hbond, such as the angle or distance between electronegative atoms.6,31 We present a somewhat different, pragmatic, and straightforward methodology for extracting the free energy associated with Hbond formation in all-atom molecular dynamics (MD) simulations. This method can be easily applied also to Monte Received: March 6, 2017 Published: May 10, 2017 2851

DOI: 10.1021/acs.jctc.7b00238 J. Chem. Theory Comput. 2017, 13, 2851−2857

Article

Journal of Chemical Theory and Computation Carlo and ab initio simulations. In accordance with previous studies,17,18,31,32 we use the underlying spatial (distance and angular orientation) probability distribution function of donor− acceptor pairs as the basic quantity in our calculations. Then, the free energy change upon Hbond formation is calculated as the reversible work associated with forging this configurational distribution from a reference random (“non-Hbond-interacting”) distribution. Hence by this definition, hydrogen bonding is associated with the forces leading to the nonrandom spatial orientation of a donor−acceptor pair. Importantly, the calculation of the Hbond free energy does not require separation into distinct (e.g., bound vs unbound) populations, and the only distance cutoffs required to define the Hbond extent can be determined directly from properties of the distributions. Thus, the method does not introduce any external energetic calibration that is not already evident in the simulation itself. The same calculation can be performed for any kind of donor− acceptor Hbond pair, including moieties from molecules other than water (sugars, proteins, etc.) without resorting to advance knowledge of the requirements for Hbond formations between and within these molecules. The underlying entropic and enthalpic contributions to the free energy can then be extracted from the temperature dependence of the free energy witnessed in a series of simulations performed at different temperatures. Our methodology should be particularly helpful in determining the strength of Hbonds in mixtures that contain more than one Hbond former (unlike, e.g., pure water). In such mixtures, evaluating Hbond strength is even more involved than in pure water, as a larger number of donor−acceptor pair types are available. As a relevant demonstration, we focus on the importance of Hbond modification in aqueous solutions of osmolytes (small naturally occurring molecules used to adapt to environmental stresses33), where it has been established that the water−water Hbonds are stronger than in pure water.34−37 Of particular interest among osmolytes is the nonreducing disaccharide trehalose, α-D-glucopyranosyl-α-D-glucopyranoside.38 Not only does trehalose help stabilize proteins in their native state,39−41 but it also confers survival capabilities under extreme conditions to many forms of life through a state of suspended animation termed “anhydrobiosis”.42,43 Hbonding in trehalose solutions have been previously studied, indicating that trehalose−trehalose Hbonds in aqueous solutions lead to trehalose aggregation already at moderate concentrations.44−46 In the following sections we first describe in detail the methodology for evaluating the free energy of Hbond formation in simulations and then present an implementation for hydrogen bonding in water and in water-trehalose binary solutions. Finally, where available, we discuss and compare the findings with corresponding experiments.

Figure 1. Hydrogen bonding between molecules. A) Hydrogen bond formed between two water molecules. The bond that can form between a hydroxyl group and an oxygen atom can be quantified by the interoxygen distance r and the angle θ. B) Trehalose consists of 8 hydroxyl groups and 3 ethereal oxygens that can form Hbonds with water or with other trehalose molecules. Molecules are rendered in the CPK representation.

involving a hydroxyl group and an oxygen atom, for example Hbonds between water and sugar hydroxyl groups, Figure 1B. Similarly, these parameters can also be used to describe other Hbonds that may form between other electronegative atoms, including fluorine or nitrogen.47,48 2.2. Hbond Probability Density and the Potential of Mean Force. The fundamental quantity we use to determine Hbond stability from MD simulations is the probability density P(r,θ) for a specific Hbond configuration defined by r and θ. Figure 2A shows P(r,θ) for Hbonds formed between water molecules in simulations of pure water. The distribution shows a distinct population at small (r,θ) that can be ascribed to bona fide, “well-formed” Hbonds, to contrast with other “broken” Hbond contacts. Such delineation has previously been used to afford a qualitative description of Hbond formation, as well as its dependence on changes in environmental variables, such as the temperature and solution composition.36,45,49−51 However, while a useful distinction can be made between formed and unformed bonds based on the probability distribution, we show in the following that such a distinction is unnecessary in order to evaluate the free energy associated with the Hbond. The distribution of hydrogen bonds includes two important free energy contributions. First, there is a spatial (orientational and translational) degeneracy that favors random orientations. Second, there is a preference for linear alignment of the O−H···O link, which tends to reduce θ from its values in the random distribution. Hence, we define a natural reference state as the completely random distribution, Prand(r,θ), in the absence of any aligning tendency. The Potential of Mean Force, PMFHB(r,θ), is then defined as

2. METHODOLOGY 2.1. Definitions. We characterize an Hbond between a hydroxyl group “donor” (or a hydrogen connected to some other electronegative atom) and an oxygen (or another electronegative) “acceptor” atom by two parameters: the interoxygen distance, r, and the angle formed between the hydroxyl and the vector connecting the oxygen atoms, θ, Figure 1A. We note that this choice is not restrictive, and any other set of parameters could also be used within our methodology. For simplicity, when more than one hydrogen atom can form an Hbond between the oxygen atoms (for example, when both are hydroxyls), only the Hbond with the smallest θ is considered. The same parameters are useful both for water−water Hbonds and for any other bond

PMFHB(r , θ ) = −RT ln

P(r , θ ) Prand(r , θ )

(1)

where R is the gas constant, and T is the absolute temperature. The random distribution Prand(r,θ) is the product of the spherical-shell distance degeneracy (4πr2) and the random Hbond orientational distribution with respect to θ, Prand(θ), which in turn depends on the number of available hydrogen atoms in the donor−acceptor pair. The Hbond PMF has previously been used to describe the free energy change 2852

DOI: 10.1021/acs.jctc.7b00238 J. Chem. Theory Comput. 2017, 13, 2851−2857

Article

Journal of Chemical Theory and Computation

Figure 2. Hbond configurational landscape in pure water. A) Probability density function, P(r,θ), for water−water Hbonds at 298 K with respect to θ and r. B) Potential of mean force for Hbond formation, PMFHB, shown for the same solution as in panel (A). C) A cross section of panel (B) along r in the range 0° < θ < 3° yields PMFHB as a function of r for strictly linear Hbond orientations.

associated with the formation of an acceptor−donor configuration with respect to the random probability (see, e.g., ref 20). Figure S1 in the Supporting Information (SI) shows Prand(θ) for several Hbond donor−acceptor pairs that are frequently encountered in aqueous solutions. For example, in the simple case of a single hydroxyl interacting with an isolated electronegative atom, such as an ethereal oxygen, the distribution simply follows sinθ, representing the angular degeneracy of Prand(θ). More generally, we calculated numerically these distributions by randomly sampling the distributions of the donor and acceptor groups, given a set of constraints, i.e. the number of hydroxyls and the intramolecular angles. We note that eq 1 is exact only if both P(r,θ) and Prand(r,θ) extend to the bulk (typically above 10 Å). If the distribution is truncated before the bulk values, the expression can still be corrected by simply adding the remainder term −RT ln[g(r)∫ π0Prand(r,θ)dθ/∫ π0P(r,θ)dθ] to PMFHB in eq 1, where g(r) is the radial distribution function that does extend to bulk solution, g(r) = ∫ π0P(r,θ)dθ/Z(r), and defining Z(r) = ∫ π0Prand(r,θ)dθ. Figure 2B shows PMFHB(r,θ) for water−water Hbonds in pure water at T = 298 K. Notably, as has previously been determined,20,28 the configuration with lowest PMFHB is the completely linear arrangement between the two oxygens and the hydrogen, so that the global minimum is found at θ = 0. In fact, the corresponding P(r,θ), Figure 2A, displays a prominent peak at θ ≈ 10° only because of the angular degeneracy that favors angles away from θ = 0, Figure S1. We can now follow the free energy associated with the transition between any configuration in the bulk to the global minimum found in the PMFHB, shown for example in the cross section in Figure 2C for a transition from large r distances to close distances, at which the Hbond forms. This hypothetical process is associated with a free energy difference, ΔGPMF. The underlying entropic and enthalpic contributions, ΔSPMF and ΔHPMF, respectively, can be obtained from the tangent of ΔGPMF vs temperature. 2.3. Hbond Free Energy. To evaluate Hbond strength, we consider the reversible work required to form the specific probability distribution of neighboring electronegative atoms found in simulations that is different from the random distribution. This reversible work can then be interpreted as the free energy of Hbond formation because it takes into consideration the entire configurational ensemble of the first solvation shell, rather than a single specific configuration. Within this scheme, the first solvation shell contains all the water molecules that, because of their proximity, have the potential to form an Hbond. The idea is to consider each PMFHB(r,θ) term weighed by

Figure 3. Thermodynamics of hydrogen bonds. A) The weighted potential of mean force, P(r,θ)·PMFHB. B) The free energy ΔG associated with creating the Hbond configurational ensemble, as a function of temperature, for the different types of Hbond studied.

the probability for that particular configuration, P(r,θ)·PMFHB, see example for water−water Hbonds in Figure 3A. Integration over the configurational ensemble yields the free energy for Hbond formation π

ΔG = −RT

∫0 ∫r

rmax

min

ζP(r , θ ) ln

ζP(r , θ ) dθdr ζrandPrand(r , θ ) (2)

Here we integrate only over the volume of the first solvation shell, so that the distribution functions need to be rescaled by π rmax 1/ζ = ∫ ∫ P(r , θ )dθdr and by the corresponding expres0

rmin

sion for 1/ζrand. The range of integration over r, i.e. the range that constitutes the first solvation shell, is the only cutoff we need define in our scheme. We choose the limits of integration, [rmin, rmax], to be those for which P(r,θ)·PMFHB = 0 for θ → 0, 2853

DOI: 10.1021/acs.jctc.7b00238 J. Chem. Theory Comput. 2017, 13, 2851−2857

Article

Journal of Chemical Theory and Computation

We note in passing that the PMFHB can be also used to infer kinetic parameters of water molecules involved in Hbonding. Hbond formation and breakage will follow paths along the PMFHB surface, from which transition states and kinetic parameters can in principle be extracted. Specifically, the saddle point at ca. (θ = 38, r = 3.2) suggests a transition state for the process of a certain water molecule exchanging an Hbond partner. The work of Laage and Hynes65,66 suggests that water reorientation occurs via a large amplitude angular jump. It is tempting to speculate that such a jump occurs between the two basins in the PMFHB, Figure 2B and S2B. The PMFHB at the saddle point suggests a free energy of ∼12 kJ/mol at the transition state relative to the Hbonded basin, which is not far from the estimate for the free energy barrier given by Laage and Hynes of 8−9 kJ/mol.66 Note, however, that to properly extract kinetic parameters from the PMFHB, appropriate corrections should be introduced carefully,67 in a procedure that goes beyond the scope of the current work. To evaluate the strength of the Hbond we now consider the thermally averaged free energy change associated with the (hypothetical) process of reorienting the donor−acceptor ensemble from the completely random distribution. Thus, we consider the entire configurational ensemble of donor−acceptor pairs in the first solvation shell, eq 2, to determine Hbond strength, Figure 3B. Table 1 lists the related changes in thermodynamic quantities for Hbonds in pure water and in aqueous trehalose solution: ΔG, ΔH, and TΔS. We compare values for water−water (WW) Hbonds, trehalose−water Hbonds (TW), and trehalose−trehalose Hbonds formed between hydroxyl groups on different molecules (TT). The general trend in free energy change, ΔG, between the different Hbond types considered is similar to the trend in ΔGPMF, Table 1. Although the enthalpy changes that we found using eq 2 and those derived from ΔGPMF for the WW Hbond in pure water were (perhaps fortuitously) close, the associated entropic change derived from the PMFHB is quantitatively and qualitatively inconsistent with the idea of entropy loss upon Hbond formation. Indeed this inconsistency naturally arises from the fact the PMFHB scheme takes into account only one configuration (the global free energy minimum), while the thermodynamically averaged ΔG incorporates the entire Hbonded configurational ensemble (Table 1). Accordingly, the corresponding entropic drop ΔS properly accounts for the expected orientational constriction associated with the directed Hbond configuration. The values we derive from simulations for WW Hbonds in pure water are in good agreement with several experimental studies. Of particular relevance are spectroscopic experiments aimed at dissecting the measurements into contributions from two distinct populations: Hbonded and non-Hbonded. Silverstein and Dill13 have compiled available experimental data sets of Hbond stability, with estimates in the range of [−6, −13] kJ/mol for ΔH and [−12, −3] J/mol K for ΔS. Later, Smith et al.16 used X-ray absorption spectroscopy to measure the energy required to locally rearrange (or distort) Hbonds and determined it to be −6.3 ± 2.1 kJ/mol. Our analysis of simulations of pure TIP4Pew water60 gives results close to those of the relevant experiments, ΔH = −5.92 ± 0.09 kJ/mol. The entropic component we find for WW Hbonds in pure water is negative, ΔS = −7.5 ± 0.3 J/mol K. This negative entropy agrees well with the range of estimates reported by Silverstein and Dill.13 We find that water−water Hbonds are somewhat stronger (by ca. 3%) in trehalose solution than they are in pure water,

Figure 3A. This distance range, common also in other analysis schemes,52 is analogous to integrating over the first solvation shell for acceptor−donor interaction, as evident by the first peak of the radial distribution function, g(r). By contrast, another population noticeable at larger r (see Figures 2A-B and Figure 3A) is already part of the second solvation shell and therefore does not describe a direct Hbond, see Figure S2 in the SI for further details. The populations at r < rmin are excluded from the integration because they represent states that are inaccessible due to the repulsive (steric) interaction between electronegative atoms. The elimination of additional regions due to steric interactions is possible, as long as the random distribution is normalized accordingly. In the examples studied here we did not find large differences in the accessible configuration space between different studied Hbonds and accordingly conserved the integration limits according to the criteria given above. The expression for ΔG in eq 2 is fundamental in our methodology. It describes the free energy associated with the formation of preferentially oriented hydrogen bonds out of an orientationally completely random distribution of donors and acceptors with the same volume density. We note that eq 2 is analogous to the way free energy changes are calculated in the information-theoretic approach by using the surprisal of a given distribution.53,54 The entropic and enthalpic contributions to this free energy change, ΔS and ΔH, are derived from linear fits to ΔG vs T. We have also tested more elaborate nonlinear relationships between ΔG and T but found that they give the same results within statistical error. 2.4. Simulations Details. To demonstrate how we extract thermodynamic properties from simulations of mixtures containing Hbonds, we performed molecular dynamics (MD) simulations of two systems: pure water solution and 0.8m trehalose solution, at several temperatures in the range [289, 307] K. All MD simulations were performed using the GROMACS package55−59 using the TIP4Pew water model60 (a variant of TIP4P61) and the CHARMM36 force field for trehalose.62,63 See Section 1 in the SI for more details of the computational methods.

3. RESULTS AND DISCUSSION The calculation of PMFHB is an important step in the procedure we use to evaluate ΔG. As described in section 2.2, we can ascribe a meaning to any transition along a path in the PMFHB(r,θ) landscape between any two configurations. Of specific interest is the global minimum, which represents the configuration with lowest free energy with respect to the value at large distances (bulk). In pure water, the transition between the bulk and the global minimum of the Hbond configurational free energy landscape indicates a free energy change of ΔGPMF = −10.26 ± 0.02 kJ/mol, see PMFHB in Figure 2B−C and Table 1. From the temperature dependence of ΔGPMF, we resolve the corresponding ΔHPMF = −5.6 ± 0.7 kJ/mol and ΔSPMF = 16 ± 2 J/mol K. Interestingly, this positive entropy change can be related to an apparent increase in the number of (degenerate) accessible states Ω = 6.6 through the Boltzmann entropy S = k ln Ω. This apparent number of states is close to the six indistinguishable configurations of two neighboring water molecules expected for a perfect tetrahedral arrangement, which can be derived in a way similar to that described in refs 1 and 64. In the context of the PMFHB, this can be understood as the number of degenerate realizations of a bonded donor−acceptor pair at a given angle. 2854

DOI: 10.1021/acs.jctc.7b00238 J. Chem. Theory Comput. 2017, 13, 2851−2857

Article

Journal of Chemical Theory and Computation Table 1. Hydrogen Bond Thermodynamic Parameters from MD Simulationsa system

Hbond type

ΔGPMF

ΔG

ΔH

TΔS

pure water 0.8 M trehalose

WW WW TW TT

−10.26 (0.02) −10.67 (0.03) −9.08 (0.07) −9.0 (0.2)

−3.7 (0.1) −3.8 (0.2) −3.5 (0.5) −3.4 (0.9)

−5.92 (0.09) −6.0 (0.2) −3.6 (0.4) −2.8 (0.6)

−2.24 (0.09) −2.2 (0.2) −0.2 (0.4) 0.6 (0.6)

All quantities are given in units of kJ/mol and at T = 298 K. Numbers in brackets represent standard deviations. Values of ΔGPMF were derived from eq 1. Other values are derived from eq 2 and the temperature dependence of ΔG. a

determination of energetic and entropic components of the free energy of hydrogen bond formation. Akin to the Born energy of solvation, the hydrogen bond energy evaluated here is the free energy associated with forming (or “charging up”) the strongly directed bond starting from a completely random orientational configuration. This is somewhat different than the “unbound” reference state used in many other studies of the hydrogen bond. The methodology should prove useful in future endeavors to study the role of Hbonds in complex solutions that contain many nonaqueous components, such as solutions of interacting biomacromolecules, as well as deep eutectic solvents (DES). When applied to MD simulations performed over a range of temperatures and concentrations, this scheme should help to elucidate the role of Hbonds in molecular mechanisms underlying different solvation phenomena.

Table 1. This observation is consistent with other computational45 as well as experimental68 studies, that found shorter and more linear water−water hydrogen bonds in sugar solutions as well as in the presence of other osmolytes.34−37 By contrast, the trehalose-associated Hbonds are weaker than WW Hbonds, with TT Hbonds found to be the weakest. Still, these solute−solute Hbonds can lead to appreciable self-association of trehalose,44−46 as in the case of other small osmolytes.37,49,69 Notably, the enthalpic component of trehalose Hbonds, that is significantly smaller than in water, is partially compensated by a corresponding smaller entropic loss, so that ΔG for TW and TT Hbonds is only ∼8−10% different from water−water hydrogen bonds, Table 1. Interestingly, in the case of TT Hbonds, the entropy change is close to zero (within statistical error) or even favorable, which is unintuitive. Although resolving the precise underlying reason would require further investigation, we can propose several possibilities: (i) a hydrophobic-like contribution, i.e. the Hbonding between trehalose molecules may affect water structure around the molecules; (ii) upon TT Hbond formation, water molecules previously bound to these molecules may be “released”, thereby effectively contributing to the Hbond entropy increase; or (iii) since the TT Hbond is energetically weaker than WW bonds, there may be little change in the sugars’ angular distribution and so little entropy to be lost from directing the bond.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge. The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.7b00238. Details of the molecular dynamics simulations, plot of the random distributions for the configurations of Hbondforming donor and acceptor pairs, and figure of the radial distribution function of oxygen atoms in pure water solution along with the respective free energy plots (PDF)

4. SUMMARY AND CONCLUSIONS We present a methodology to quantitatively analyze Hbond thermodynamics in atomistic simulations, including MD and other computational techniques. First, we evaluate the probabilities in the configurational ensemble of donor−acceptor pairs in the first solvation shell and compare it to the reference ensemble of completely random pairs. This yields the donor− acceptor configurational free energy landscape, PMFHB. Then, the free energy of formation of the Hbond is defined as the reversible work of forging the observed distribution from the completely random distribution, as defined in eq 2. This definition pleasingly accounts for the formation of Hbonds with preferred orientations from the entire configurational ensemble of acceptor− donor pairs. The methodology we present here affords several advantages. First, the method does not introduce any externally calibrated force fields to analyze the strength of hydrogen bonds but rather relies on distributions of molecular orientations as determined in the simulations themselves to evaluate the free energy. The methodology is moreover easily extended to any molecular moieties, since it does not assume any a priori constraints on the definition of the hydrogen bond beyond the fundamental IUPAC definition.47,48 Specifically, the method does not require selection of hydrogen bonded or unbonded-pair populations. In contrast to several other proposed methodologies, the only requirement is that the configurational donor−acceptor landscape will be well sampled. Finally, the method allows an easy

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Daniel Harries: 0000-0002-3057-9485 Present Address †

Faculty of Mechanical Engineering, TechnionIsrael Institute of Technology, Haifa 32000, Israel. Author Contributions

The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS

The financial support from the Israel Science Foundation (ISF grant No. 1538/13) is gratefully acknowledged. The Fritz Haber Research Center is supported by the Minerva Foundation, Munich, Germany. 2855

DOI: 10.1021/acs.jctc.7b00238 J. Chem. Theory Comput. 2017, 13, 2851−2857

Article

Journal of Chemical Theory and Computation

■

(25) Stumpe, M. C.; Grubmuller, H. Aqueous Urea Solutions: Structure, Energetics, and Urea Aggregation. J. Phys. Chem. B 2007, 111, 6220−6228. (26) Markovitch, O.; Agmon, N. Structure and Energetics of the Hydronium Hydration Shells. J. Phys. Chem. A 2007, 111, 2253−2256. (27) Silverstein, K. a T.; Haymet, a D. J.; Dill, K. a. Molecular Model of Hydrophobic Solvation. J. Chem. Phys. 1999, 111, 8000. (28) Sharp, K. A.; Vanderkooi, J. M. Water in the Half Shell: Structure of Water, Focusing on Angular Structure and Solvation. Acc. Chem. Res. 2010, 43, 231−239. (29) Sciortino, F.; Fornili, S. L. Hydrogen Bond Cooperativity in Simulated Water: Time Dependence Analysis of Pair Interactions. J. Chem. Phys. 1989, 90, 2786−2792. (30) Tropsha, a; Bowen, J. P.; Brown, F. K.; Kizer, J. S. Do Interhelical Side Chain-Backbone Hydrogen Bonds Participate in Formation of Leucine Zipper Coiled Coils? Proc. Natl. Acad. Sci. U. S. A. 1991, 88, 9488−9492. (31) Sharp, K. A.; Madan, B. Hydrophobic Effect, Water Structure, and Heat Capacity Changes. J. Phys. Chem. B 1997, 101, 4343−4348. (32) Vanzi, F.; Madan, B.; Sharp, K. Effect of the Protein Denaturants Urea and Guanidinium on Water Structure: A Structural and Thermodynamic Study. J. Am. Chem. Soc. 1998, 120, 10748−10753. (33) Hochachka, P. W.; Somero, G. N. Biochemical Adaptation; Oxford University Press: New York, NY, 2002; pp 217−289. (34) Sharp, K. A.; Madan, B.; Manas, E.; Vanderkooi, J. M. Water Structure Changes Induced by Hydrophobic and Polar Solutes Revealed by Simulations and Infrared Spectroscopy. J. Chem. Phys. 2001, 114, 1791−1796. (35) Zou, Q.; Bennion, B. J.; Daggett, V.; Murphy, K. P. The Molecular Mechanism of Stabilization of Proteins by TMAO and Its Ability to Counteract the Effects of Urea. J. Am. Chem. Soc. 2002, 124, 1192−1202. (36) Politi, R.; Sapir, L.; Harries, D. The Impact of Polyols on Water Structure in Solution: A Computational Study. J. Phys. Chem. A 2009, 113, 7548−7555. (37) Towey, J. J.; Dougan, L. Structural Examination of the Impact of Glycerol on Water Structure. J. Phys. Chem. B 2012, 116, 1633−1641. (38) Elbein, A. D.; Pan, Y. T.; Pastuszak, I.; Carroll, D. New Insights on Trehalose: A Multifunctional Molecule. Glycobiology 2003, 13, 17R− 27R. (39) Xie, G. F.; Timasheff, S. N. The Thermodynamic Mechanism of Protein Stabilization by Trehalose. Biophys. Chem. 1997, 64, 25−43. (40) Jain, N. K.; Roy, I. Effect of Trehalose on Protein Structure. Protein Sci. 2009, 18, 24−36. (41) Politi, R.; Harries, D. Enthalpically Driven Peptide Stabilization by Protective Osmolytes. Chem. Commun. 2010, 46, 6449−6451. (42) Crowe, L. M.; Reid, D. S.; Crowe, J. H. Is Trehalose Special for Preserving Dry Biomaterials? Biophys. J. 1996, 71, 2087−2093. (43) Clegg, J. S. Cryptobiosis  a Peculiar State of Biological Organization. Comp. Biochem. Physiol., Part B: Biochem. Mol. Biol. 2001, 128, 613−624. (44) Lerbret, A.; Bordat, P.; Affouard, F.; Descamps, M.; Migliardo, F. How Homogeneous Are the Trehalose, Maltose, and Sucrose Water Solutions? An Insight from Molecular Dynamics Simulations. J. Phys. Chem. B 2005, 109, 11046−11057. (45) Sapir, L.; Harries, D. Linking Trehalose Self-Association with Binary Aqueous Solution Equation of State. J. Phys. Chem. B 2011, 115, 624−634. (46) Winther, L. R.; Qvist, J.; Halle, B. Hydration and Mobility of Trehalose in Aqueous Solution. J. Phys. Chem. B 2012, 116, 9196−9207. (47) Arunan, E.; Desiraju, G. R.; Klein, R. a.; Sadlej, J.; Scheiner, S.; Alkorta, I.; Clary, D. C.; Crabtree, R. H.; Dannenberg, J. J.; Hobza, P.; Kjaergaard, H. G.; Legon, A. C.; Mennucci, B.; Nesbitt, D. J. Definition of the Hydrogen Bond (IUPAC Recommendations 2011). Pure Appl. Chem. 2011, 83, 1637−1641. (48) Arunan, E.; Desiraju, G. R.; Klein, R. A.; Sadlej, J.; Scheiner, S.; Alkorta, I.; Clary, D. C.; Crabtree, R. H.; Dannenberg, J. J.; Hobza, P.; Kjaergaard, H. G.; Legon, A. C.; Mennucci, B.; Nesbitt, D. J. Defining the Hydrogen Bond: An Account (IUPAC Technical Report). Pure Appl. Chem. 2011, 83, 1619−1636.

REFERENCES

(1) Pauling, L. The Nature of the Chemical Bond and the Structure of Molecules and Crystals: An Introduction to Modern Structural Chemistry, 3rd ed.; Cornell University Press: 1960. (2) Chandler, D. Interfaces and the Driving Force of Hydrophobic Assembly. Nature 2005, 437, 640−647. (3) Hindman, J. C. Proton Resonance Shift of Water in the Gas and Liquid States. J. Chem. Phys. 1966, 44, 4582. (4) Walrafen, G. E.; Fisher, M. R.; Hokmabadi, M. S.; Yang, W.-H. Temperature Dependence of the Low- and High-Frequency Raman Scattering from Liquid Water. J. Chem. Phys. 1986, 85, 6970−6982. (5) Hare, D. E.; Sorensen, C. M. Raman Spectroscopic Study of Dilute HOD in Liquid H2O in the Temperature Range − 31.5 to 160 °C. J. Chem. Phys. 1990, 93, 6954. (6) Espinosa, E.; Molins, E.; Lecomte, C. Hydrogen Bond Strengths Revealed by Topological Analyses of Experimentally Observed Electron Densities. Chem. Phys. Lett. 1998, 285, 170−173. (7) Carey, D. M. Measurement of the Raman Spectrum of Liquid Water. J. Chem. Phys. 1998, 108, 2669−2675. (8) Suresh, S. J.; Naik, V. M. Hydrogen Bond Thermodynamic Properties of Water from Dielectric Constant Data. J. Chem. Phys. 2000, 113, 9727−9732. (9) Modig, K.; Pfrommer, B.; Halle, B. Temperature-Dependent Hydrogen-Bond Geometry in Liquid Water. Phys. Rev. Lett. 2003, 90, 25−28. (10) Starr, F. W.; Nielsen, J. K.; Stanley, H. E. Hydrogen-Bond Dynamics for the Extended Simple Point-Charge Model of Water. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 2000, 62, 579− 587. (11) van der Spoel, D.; van Maaren, P. J.; Larsson, P.; Timneanu, N. Thermodynamics of Hydrogen Bonding in Hydrophilic and Hydrophobic Media Thermodynamics of Hydrogen Bonding in Hydrophilic and Hydrophobic Media. J. Phys. Chem. B 2006, 110, 4393−4398. (12) Hakem, I. F.; Boussaid, A.; Benchouk-Taleb, H.; Bockstaller, M. R. Temperature, Pressure, and Isotope Effects on the Structure and Properties of Liquid Water: A Lattice Approach. J. Chem. Phys. 2007, 127, 224106. (13) Silverstein, K. a T.; Haymet, a D. J.; Dill, K. a. The Strength of Hydrogen Bonds in Liquid Water and Around Nonpolar Solutes. J. Am. Chem. Soc. 2000, 122, 8037−8041. (14) Chaplin, M. Water’s Hydrogen Bond Strength. arXiv 2007, 1−20. (15) Muller, N. Search for A Realistic View of Hydrophobic Effects. Acc. Chem. Res. 1990, 23, 23−28. (16) Smith, J. D.; Cappa, C. D.; Wilson, K. R.; Messer, B. M.; Cohen, R. C.; Saykally, R. J. Energetics of Hydrogen Bond Network Rearrangements in Liquid Water. Science 2004, 306, 851−853. (17) Pople, J. A. Molecular Association in Liquids. II. A Theory of the Structure of Water. Proc. R. Soc. London, Ser. A 1951, 205, 163−178. (18) Sceats, M. G.; Stavola, M.; Rice, S. A. A Zeroth Order Random Network Model of Liquid Water. J. Chem. Phys. 1979, 70, 3927−3938. (19) Sceats, M. G.; Rice, S. A. The Enthalpy and Heat Capacity of Liquid Water and Ice Polymorphs from a Random Network Model. J. Chem. Phys. 1980, 72, 3248−3259. (20) Kumar, R.; Schmidt, J. R.; Skinner, J. L. Hydrogen Bonding Definitions and Dynamics in Liquid Water. J. Chem. Phys. 2007, 126, 204107. (21) McDonald, I. K.; Thornton, J. M. Satisfying Hydrogen-Bonding Potential in Proteins. J. Mol. Biol. 1994, 238, 777−793. (22) Luzar, A.; Chandler, D. Effect of Environment on Hydrogen Bond Dynamics in Liquid Water. Phys. Rev. Lett. 1996, 76, 928−931. (23) Gilman-Politi, R.; Harries, D. Unraveling the Molecular Mechanism of Enthalpy Driven Peptide Folding by Polyol Osmolytes. J. Chem. Theory Comput. 2011, 7, 3816−3828. (24) Lerbret, A.; Affouard, F.; Hédoux, A.; Krenzlin, S.; Siepmann, J.; Bellissent-Funel, M.-C.; Descamps, M. How Strongly Does Trehalose Interact with Lysozyme in the Solid State? Insights from Molecular Dynamics Simulation and Inelastic Neutron Scattering. J. Phys. Chem. B 2012, 116, 11103−11116. 2856

DOI: 10.1021/acs.jctc.7b00238 J. Chem. Theory Comput. 2017, 13, 2851−2857

Article

Journal of Chemical Theory and Computation (49) Dashnau, J. L.; Nucci, N. V.; Sharp, K. A.; Vanderkooi, J. M. Hydrogen Bonding and the Cryoprotective Properties of Glycerol/ water Mixtures. J. Phys. Chem. B 2006, 110, 13670−13677. (50) Dashnau, J. L.; Vanderkooi, J. M. Computational Approaches to Investigate How Biological Macromolecules Can Be Protected in Extreme Conditions. J. Food Sci. 2007, 72, R001−R010. (51) Gallagher, K. R.; Sharp, K. A. Analysis of Thermal Hysteresis Protein Hydration Using the Random Network Model. Biophys. Chem. 2003, 105, 195−209. (52) Gallagher, K. R.; Sharp, K. a A New Angle on Heat Capacity Changes in Hydrophobic Solvation. J. Am. Chem. Soc. 2003, 125, 9853− 9860. (53) Ben-Shaul, A.; Levine, R. D.; Bernstein, R. B. Entropy and Chemical Change. II. Analysis of Product Energy Distributions: Temperature and Entropy Deficiency. J. Chem. Phys. 1972, 57, 5427− 5447. (54) Procaccia, I.; Levine, R. D. Potential Work: A StatisticalMechanical Approach for Systems in Disequilibrium. J. Chem. Phys. 1976, 65, 3357. (55) Hess, B.; Kutzner, C.; Van Der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435−447. (56) Van Der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. GROMACS: Fast, Flexible, and Free. J. Comput. Chem. 2005, 26, 1701−1718. (57) Lindahl, E.; Hess, B.; van der Spoel, D. GROMACS 3. 0: A Package for Molecular Simulation and Trajectory Analysis. J. Mol. Model. 2001, 7, 306−317. (58) Berendsen, H. J. C.; van der Spoel, D.; van Drunen, R. ROMACS: A Message-Passing Parallel Molecular Dynamics Implementation. Comput. Phys. Commun. 1995, 91, 43−56. (59) Bjelkmar, P.; Larsson, P.; Cuendet, M. a.; Hess, B.; Lindahl, E. Implementation of the CHARMM Force Field in GROMACS: Analysis of Protein Stability Effects from Correction Maps, Virtual Interaction Sites, and Water Models. J. Chem. Theory Comput. 2010, 6, 459−466. (60) Horn, H. W.; Swope, W. C.; Pitera, J. W.; Madura, J. D.; Dick, T. J.; Hura, G. L.; Head-Gordon, T. Development of an Improved FourSite Water Model for Biomolecular Simulations: TIP4P-Ew. J. Chem. Phys. 2004, 120, 9665−9678. (61) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926. (62) Guvench, O.; Greene, S. N.; Kamath, G.; Brady, J. W.; Venable, R. M.; Pastor, R. W.; Mackerell, A. D. Additive Empirical Force Field for Hexopyranose Monosaccharides. J. Comput. Chem. 2008, 29, 2543− 2564. (63) Guvench, O.; Hatcher, E.; Venable, R. M.; Pastor, R. W.; MacKerell, A. D. CHARMM Additive All-Atom Force Field for Glycosidic Linkages between Hexopyranoses. J. Chem. Theory Comput. 2009, 5, 2353−2370. (64) DiMarzio, E. A.; Stillinger, F. H. Residual Entropy of Ice. J. Chem. Phys. 1964, 40, 1577. (65) Laage, D. A Molecular Jump Mechanism of Water Reorientation. Science (Washington, DC, U. S.) 2006, 311, 832−835. (66) Laage, D.; Hynes, J. T. On the Molecular Mechanism of Water Reorientation. J. Phys. Chem. B 2008, 112, 14230−14242. (67) Schenter, G. K.; Garrett, B. C.; Truhlar, D. G. Generalized Transition State Theory in Terms of the Potential of Mean Force. J. Chem. Phys. 2003, 119, 5828−5833. (68) Guo, F.; Friedman, J. M. Osmolyte-Induced Perturbations of Hydrogen Bonding between Hydration Layer Waters: Correlation with Protein Conformational Changes. J. Phys. Chem. B 2009, 113, 16632− 16642. (69) Towey, J.; Soper, A.; Dougan, L. Molecular Insight into the Hydrogen Bonding and Micro-Segregation of a Cryoprotectant Molecule. J. Phys. Chem. B 2012, 116, 13898−13904.

2857

DOI: 10.1021/acs.jctc.7b00238 J. Chem. Theory Comput. 2017, 13, 2851−2857