Protonation/deprotonation effects on the stability of the Trp-cage miniprotein

The Trp-cage miniprotein is a 20 amino acid peptide that exhibits many of the properties of globular proteins. In this protein, the hydrophobic core is formed by a buried Trp side chain. The folded state is stabilized by an ion pair between aspartic acid and an arginine side chain. The effect of protonating the aspartic acid on the Trp-cage miniprotein folding/unfolding equilibrium is studied by explicit solvent molecular dynamics simulations of the protein in the charged and protonated Asp9 states. Unbiased Replica Exchange Molecular Dynamics (REMD) simulations, spanning a wide temperature range, are carried out to the microsecond time scale, using the AMBER99SB forcefield in explicit TIP3P water. The protein structural ensembles are studied in terms of various order parameters that differentiate the folded and unfolded states. We observe that in the folded state the root mean square distance (rmsd) from the backbone of the NMR structure shows two highly populated basins close to the native state with peaks at 0.06 nm and 0.16 nm, which are consistent with previous simulations using the same forcefield. The fraction of folded replicas shows a drastic decrease because of the absence of the salt bridge. However, significant populations of conformations with the arginine side chain exposed to the solvent, but within the folded basin, are found. This shows the possibility to reach the folded state without formation of the ion pair. We also characterize changes in the unfolded state. The equilibrium populations of the folded and unfolded states are used to characterize the thermodynamics of the system. We find that the change in free energy difference due to the protonation of the Asp amino acid is 3 kJ mol(-1) at 297 K, favoring the charged state, and resulting in ΔpK(1) = 0.5 units for Asp9. We also study the differences in the unfolded state ensembles for the two charge states and find significant changes at low temperature, where the protonated Asp side chain makes multiple hydrogen bonds to the protein backbone.     

Structural properties of water: comparison of the SPC, SPCE, TIP4P, and TIP5P models of water

Molecular-dynamics simulations were carried out for the SPC, SPCE, TIP4P, and TIP5P models of water at 298 K. From these results we determine the following quantities: the absolute entropy using the two-particle approximation, the mean lifetime of the hydrogen bond, the mean number of hydrogen bonds per molecule, and the mean energy of the hydrogen bond. From the entropy calculations we find that nearly all contributions to the total entropy originates from the orientation effects. Moreover, we determine the contributions to the total entropy which originate from the first, second, and higher solvation shells. It is interesting that the limits between solvation shells are clearly visible. The first solvation shell (0.22 < r < 0.36 nm) contributes approximately 43 J mol K to the total entropy; the second solvation shell (0.36 < r < 0.60 nm) contributes approximately 12 J mol K, while contributions from the third and other solvation shells are very small, approximately 2 J mol K in summary. This indicates that water molecules are strongly ordered up to 0.55-0.6 nm around the central water molecule, and beyond this limit the ordering diminishes. The results of calculations (entropy and hydrogen bonds) are compared with the experimental data for the choosing of the best water model. We find that the SPC and TIP4P models reproduce the best experimental values, and we recommend these models for computer simulations of the aqueous solution of biomolecules.     

Characterization of the TIP4P-Ew water model: vapor pressure and boiling point

The liquid-vapor-phase equilibrium properties of the previously developed TIP4P-Ew water model have been studied using thermodynamic integration free-energy simulation techniques in the temperature range of 274-400 K. We stress that free-energy results from simulations need to be corrected in order to be compared to the experiment. This is due to the fact that the thermodynamic end states accessible through simulations correspond to fictitious substances (classical rigid liquids and classical rigid ideal gases) while experiments operate on real substances (liquids and real gases, with quantum effects). After applying analytical corrections the vapor pressure curve obtained from simulated free-energy changes is in excellent agreement with the experimental vapor pressure curve. The boiling point of TIP4P-Ew water under ambient pressure is found to be at 370.3+/-1.9 K, about 7 K higher than the boiling point of TIP4P water (363.7+/-5.1 K; from simulations that employ finite range treatment of electrostatic and Lennard-Jones interactions). This is in contrast to the approximately +15 K by which the temperature of the density maximum and the melting temperature of TIP4P-Ew are shifted relative to TIP4P, indicating that the temperature range over which the liquid phase of TIP4P-Ew is stable is narrower than that of TIP4P and resembles more that of real water. The quality of the vapor pressure results highlights the success of TIP4P-Ew in describing the energetic and entropic aspects of intermolecular interactions in liquid water.     

Folding and thermodynamic studies of Trp-cage based on polarized force field

Two replica exchange molecular dynamics (REMD) simulations were carried out to study the thermodynamics of a 20-residue Trp-cage folding based on a newly developed polarized protein-specific charge (PPC). Starting from a fully extended conformation, Trp-cage native conformation was successfully sampled using REMD based on a 3-step PPC update. Next, the obtained Trp-cage folded conformation was then used to calculate the PPC in which another REMD was performed to explore the thermodynamic stability of Trp-cage. The theoretical melting temperature T _m of ≈325 K was found to be in close agreement with experimental melting temperature, T _m of 315 K. This indicates that the PPC was correctly predicting the temperature dependence. The current study provides a direct proof of how electrostatic polarization affects protein folding.     

Correlations in liquid water for the TIP3P-Ewald, TIP4P-2005, TIP5P-Ewald, and SWM4-NDP models

Water is one of the simplest molecules in existence, but also one of the most important in biological and engineered systems. However, understanding the structure and dynamics of liquid water remains a major scientific challenge. Molecular dynamics simulations of liquid water were performed using the water models TIP3P-Ewald, TIP4P-2005, TIP5P-Ewald, and SWM4-NDP to calculate the radial distribution functions (RDFs), the relative angular distributions, and the excess enthalpies, entropies, and free energies. In addition, lower-order approximations to the entropy were considered, identifying the fourth-order approximation as an excellent estimate of the full entropy. The second-order and third-order approximations are ~20% larger and smaller than the true entropy, respectively. All four models perform very well in predicting the radial distribution functions, with the TIP5P-Ewald model providing the best match to the experimental data. The models also perform well in predicting the excess entropy, enthalpy, and free energy of liquid water. The TIP4P-2005 and SWM4-NDP models are more accurate than the TIP3P-Ewald and TIP5P-Ewald models in this respect. However, the relative angular distribution functions of the four water models reveal notable differences. The TIP5P-Ewald model demonstrates an increased preference for water molecules to act both as tetrahedral hydrogen bond donors and acceptors, whereas the SWM4-NDP model demonstrates an increased preference for water molecules to act as planar hydrogen bond acceptors. These differences are not uncovered by analysis of the RDFs or the commonly employed tetrahedral order parameter. However, they are expected to be very important when considering water molecules around solutes and are thus a key consideration in modelling solvent entropy.     

Hydration thermodynamic properties of amino acid analogues: a systematic comparison of biomolecular force fields and water models

We present an extensive study on hydration thermodynamic properties of analogues of 13 amino acid side chains at 298 K and 1 atm. The hydration free energies DeltaG, entropies DeltaS, enthalpies DeltaH, and heat capacities Deltac(P)() were determined for 10 combinations of force fields and water models. The statistical sampling was extended such that precisions of 0.3, 0.8, 0.8 kJ/mol and 25 J/(mol K) were reached for DeltaG, TDeltaS, DeltaH, and Deltac(P)(), respectively. The three force fields used in this study are AMBER99, GROMOS 53A6, and OPLS-AA; the five water models are SPC, SPC/E, TIP3P, TIP4P, and TIP4P-Ew. We found that the choice of water model strongly influences the accuracy of the calculated hydration entropies, enthalpies, and heat capacities, while differences in accuracy between the force fields are small. On the basis of an analysis of the hydrophobic analogues of the amino acid side chains, we discuss what properties of the water models are responsible for the observed discrepancies between computed and experimental values. The SPC/E water model performs best with all three biomolecular force fields.     

Capillary waves at the liquid-vapor interface and the surface tension of water

Capillary waves occurring at the liquid-vapor interface of water are studied using molecular dynamics simulations. In addition, the surface tension, determined thermodynamically from the difference in the normal and tangential pressure at the liquid-vapor interface, is compared for a number of standard three- and four-point water models. We study four three-point models (SPC/E, TIP3P, TIP3P-CHARMM, and TIP3P-Ew) and two four-point models (TIP4P and TIP4P-Ew). All of the models examined underestimate the surface tension; the TIP4P-Ew model comes closest to reproducing the experimental data. The surface tension can also be determined from the amplitude of capillary waves at the liquid-vapor interface by varying the surface area of the interface. The surface tensions determined from the amplitude of the logarithmic divergence of the capillary interfacial width and from the traditional thermodynamic method agree only if the density profile is fitted to an error function instead of a hyperbolic tangent function.     

Secondary structure bias in generalized Born solvent models: comparison of conformational ensembles and free energy of solvent polarization from explicit and implicit solvation

The effects of the use of three generalized Born (GB) implicit solvent models on the thermodynamics of a simple polyalanine peptide are studied via comparing several hundred nanoseconds of well-converged replica exchange molecular dynamics (REMD) simulations using explicit TIP3P solvent to REMD simulations with the GB solvent models. It is found that when compared to REMD simulations using TIP3P the GB REMD simulations contain significant differences in secondary structure populations, most notably an overabundance of alpha-helical secondary structure. This discrepancy is explored via comparison of the differences in the electrostatic component of the free energy of solvation (DeltaDeltaG(pol)) between TIP3P (via thermodynamic Integration calculations), the GB models, and an implicit solvent model based on the Poisson equation (PE). The electrostatic components of the solvation free energies are calculated using each solvent model for four representative conformations of Ala10. Since the PE model is found to have the best performance with respect to reproducing TIP3P DeltaDeltaG(pol) values, effective Born radii from the GB models are compared to effective Born radii calculated with PE (so-called perfect radii), and significant and numerous deviations in GB radii from perfect radii are found in all GB models. The effect of these deviations on the solvation free energy is discussed, and it is shown that even when perfect radii are used the agreement of GB with TIP3P DeltaDeltaG(pol) values does not improve. This suggests a limit to the optimization of the effective Born radius calculation and that future efforts to improve the accuracy of GB models must extend beyond such optimizations.     

Free energy of liquid water from a computer simulation via cell theory

A method to calculate the free energy of water from computer simulation is presented. Based on cell theory, it approximates the potential energy surface sampled in the simulation by an anisotropic six-dimensional harmonic potential to model the three hindered translations and three hindered rotations of a single rigid water molecule. The potential is parametrized from the magnitude of the forces and torques measured in the simulation. The entropy of these six harmonic oscillators is calculated and summed with a conformational term to give the total entropy. Combining this with the simulation enthalpy yields the free energy. The six water models examined are TIP3P, SPC, TIP4P, SPC/E, TIP5P, and TIP4P-Ew. The results reproduce experiment well: free energies for all models are within 1.6 kJ mol(-1) and entropies are within 3.6 J K(-1) mol(-1). Approximately two-thirds of the entropy comes from translation, a third from rotation, and 5% from conformation. Vibrational frequencies match those in the experimental infrared spectrum and assist in their assignment. Intermolecular quantum effects are found to be small, with free energies for the classical oscillator lying 0.5-0.7 kJ mol(-1) higher than in the quantum case. Molecular displacements and vibrational and zero point energies are also calculated. Altogether, these results validate the harmonic oscillator as a quantitative model for the liquid state.     

Solvation of Glucose, Trehalose, and Sucrose by the Soft Sticky Dipole-Quadrupole-Octupole Water Model

Water structure around sugars modeled by partial charges is compared for soft-sticky dipole-quadrupole-octupole (SSDQO), a fast single-site multipole model, and commonly used multi-site models in Monte Carlo simulations. Radial distribution functions and coordination numbers of all the models indicate similar hydration by hydrogen-bond donor and acceptor waters. However, the new optimized SSDQO1 parameters as well as TIP4P-Ew and TIP5P predict a "lone-pair" orientation for the water accepting the sugar hydroxyl hydrogen bond that is more consistent with the limited experimental data than the "dipole" orientation in SPC/E, which has important implications for studies of the cryoprotectant properties of sugars.     

Melting temperature of ice Ih calculated from coexisting solid-liquid phases

We carried out molecular-dynamics simulations by using the two-phase coexistence method with the constant pressure, particle number, and enthalpy ensemble to compute the melting temperature of proton-disordered hexagonal ice I(h) at 1-bar pressure. Four models of water were considered, including the widely used TIP4P [W. L. Jorgensen, J. Chandrasekha, J. D. Madura, R. W. Impey, and M. L. Klein, J. Chem. Phys.79, 926 (1983)] and TIP5P [M. W. Mahoney and W. L. Jorgensen J. Chem. Phys.112, 8910 (2000)] models, as well as recently improved TIP4P and TIP5P models for use with Ewald techniques-the TIP4P-Ew [W. Horn, W. C. Swope, J. W. Pitera, J. C. Madura, T. J. Dick, G. L. Hura, and T. Head-Gordon, J. Chem. Phys.120, 9665 (2004)] and TIP5P-Ew [S. W. Rick, J. Chem. Phys.120, 6085 (2004)] models. The calculated melting temperature at 1 bar is T(m) = 229 +/- 1 K for the TIP4P and T(m) = 272.0 +/- 0.6 K for the TIP5P ice I(h), both are consistent with previous simulations based on free-energy methods. For the TIP4P-Ew and TIP5P-Ew models, the calculated melting temperature is T(m) = 257.0 +/- 1.1 K and T(m) = 253.9 +/- 1.1 K, respectively.     

Solvation free energies of amino acid side chain analogs for common molecular mechanics water models

Quantitative free energy computation involves both using a model that is sufficiently faithful to the experimental system under study (accuracy) and establishing statistically meaningful measures of the uncertainties resulting from finite sampling (precision). In order to examine the accuracy of a range of common water models used for protein simulation for their solute/solvent properties, we calculate the free energy of hydration of 15 amino acid side chain analogs derived from the OPLS-AA parameter set with the TIP3P, TIP4P, SPC, SPC/E, TIP3P-MOD, and TIP4P-Ew water models. We achieve a high degree of statistical precision in our simulations, obtaining uncertainties for the free energy of hydration of 0.02-0.06 kcal/mol, equivalent to that obtained in experimental hydration free energy measurements of the same molecules. We find that TIP3P-MOD, a model designed to give improved free energy of hydration for methane, gives uniformly the closest match to experiment; we also find that the ability to accurately model pure water properties does not necessarily predict ability to predict solute/solvent behavior. We also evaluate the free energies of a number of novel modifications of TIP3P designed as a proof of concept that it is possible to obtain much better solute/solvent free energetic behavior without substantially negatively affecting pure water properties. We decrease the average error to zero while reducing the root mean square error below that of any of the published water models, with measured liquid water properties remaining almost constant with respect to our perturbations. This demonstrates there is still both room for improvement within current fixed-charge biomolecular force fields and significant parameter flexibility to make these improvements. Recent research in computational efficiency of free energy methods allows us to perform simulations on a local cluster that previously required large scale distributed computing, performing four times as much computational work in approximately a tenth of the computer time as a similar study a year ago.     

Predicting hydration free energies using all-atom molecular dynamics simulations and multiple starting conformations

Molecular dynamics simulations in explicit solvent were applied to predict the hydration free energies for 23 small organic molecules in blind SAMPL2 test. We found good agreement with experimental results, with an RMS error of 2.82 kcal/mol over the whole set and 1.86 kcal/mol over all the molecules except several hydroxyl-rich compounds where we find evidence for a systematic error in the force field. We tested two different solvent models, TIP3P and TIP4P-Ew, and obtained very similar hydration free energies for these two models; the RMS difference was 0.64 kcal/mol. We found that preferred conformation of the carboxylic acids in water differs from that in vacuum. Surprisingly, this conformational change is not adequately sampled on simulation timescales, so we apply an umbrella sampling technique to include free energies associated with the conformational change. Overall, the results of this test reveal that the force field parameters for some groups of molecules (such as hydroxyl-rich compounds) still need to be improved, but for most compounds, accuracy was consistent with that seen in our previous tests.     

The importance of polarizability in the modeling of solubility: quantifying the effect of solute polarizability on the solubility of small nonpolar solutes in popular models of water

In recent work by Paschek [J. Chem. Phys. 120, 6674 (2004)] and others [see H. Docherty et al., J. Chem. Phys. 125, 074510 (2006) for a review] it has been suggested that, when coupled to a simple Lennard-Jones model for various small nonpolar solute molecules, the most common models of water (e.g., SPC/E and TIP4P) fail to reproduce quantitatively the solubility of small nonpolar solute molecules in water due in part to failing to account for polarization of the solute molecule. Given the importance of such systems as test-case prototype models of the solubility of proteins and biomolecules, in this work, we investigate the impact of using a polarizable solute model with the SPC/E, TIP3P, TIP4P, TIP4P-Ew, and TIP4P/2005 rigid water models. Specifically we consider Ne, Ar, Kr, Xe, and methane as solutes. In all cases we observe that the use of a polarizable solute improves agreement between experiment and simulations, with the best agreement seen for the largest solutes, Kr, CH(4), and Xe and the modern reparametrizations of the TIP4P model, i.e., the TIP4P-Ew and TIP4P/2005 models.     

Computer simulation of the hydrogen bond lifetime and the mechanism of the structural rearrangement of water

The concept of a hydrogen bond lifetime is analyzed in a computer experiment. It is shown that different definitions of the lifetime of a hydrogen bond characterize definite stages in the microdynamics of a liquid. The lifetimes of hydrogen bonds are calculated for water using the TIP4P and TIP4B-HB model potentials. A number of features in the mechanism of the structural rearrangement of the nearest surrounding is determined by comparing them.     

Simulation of the pressure and temperature folding/unfolding equilibrium of a small RNA hairpin

We report molecular dynamics simulations of the equilibrium folding/unfolding thermodynamics of the RNA tetraloop in explicit solvent. A replica exchange molecular dynamics study of the r(CGUUGCCG) oligomer that forms a hairpin is performed for 226 ns per replica, using 52 replicas. We are able to show the unbiased folding of all replicas starting from extended conformations. The equilibrium pressure-temperature free energy of folding, DeltaG(P,T), is calculated from the averaged energy, pressure, and specific volume change upon folding of the oligomer as a function of T at constant volume. We find that this oligomer is destabilized by increasing hydrostatic pressure, similar to the behavior of globular proteins.     

Salt effects on poly(N-isopropylacrylamide) phase transition thermodynamics from NMR spectroscopy

NMR spectra were collected for cross-linked poly(N-isopropylacrylamide), poly(NIPAM), hydrogels in the presence of NaCl and CaCl2 aqueous solutions. Intensity variations in the 1H NMR signals of the polymer provide insight into the phase transition process. These data were used to observe a two-stage phase transition process. Thermodynamic quantities were obtained from a van't Hoff analysis of the temperature-dependent equilibrium constants, which were derived from the NMR data. The Delta H degrees and Delta S degrees values for the hydrogel in D2O are 3.4 kJ/mol and 11.2 J/mol.K for stage I, which is attributed to the formation of hydrophobic bonds between neighboring isopropyl groups. The formation of hydrogen bonds during stage II yielded Delta H degrees and Delta S degrees values of 14.8 kJ/mol and 48.4 J/mol.K in D2O. However, the corresponding Delta H degrees values in 150 mM NaCl and 150 mM CaCl2 are reduced to 1.5 and 1.8 kJ/mol for stage I of the dehydration process. This corresponds to the known effect of salts on hydrophobic bond energetics. The value of Delta S degrees also decreased to 4.9 and 5.9 J/mol.K in NaCl and CaCl2 solutions, respectively. However, the thermodynamic values during stage II were only slightly affected by the salts. The lower temperatures required to induce spontaneous precipitation implies that Delta G degrees of precipitation is reduced. With our measurement of equilibrium thermodynamics, we see that 150 mM NaCl and CaCl2 solutions have a greater effect on hydrophobic bond formation associated with the phase transition process. In this manner, these salts aid in solvent reorganization necessary to form the hydrophobic bond, and this suggests that the formation of hydrophobic bonds is a strong determining factor in the stability of poly(NIPAM) hydrogels in water.     

A Single‐Molecule Perspective on the Role of Solvent Hydrogen Bonds in Protein Folding and Chemical Reactions

We present an array of force spectroscopy experiments that aim to identify the role of solvent hydrogen bonds in protein folding and chemical reactions at the single‐molecule level. In our experiments we control the strength of hydrogen bonds in the solvent environment by substituting water (H₂O) with deuterium oxide (D₂O). Using a combination of force protocols, we demonstrate that protein unfolding, protein collapse, protein folding and a chemical reaction are affected in different ways by substituting H₂O with D₂O. We find that D₂O molecules form an integral part of the unfolding transition structure of the immunoglobulin module of human cardiac titin, I27. Strikingly, we find that D₂O is a worse solvent than H₂O for the protein I27, in direct contrast with the behaviour of simple hydrocarbons. We measure the effect of substituting H₂O with D₂O on the force dependent rate of reduction of a disulphide bond engineered within a single protein. Altogether, these experiments provide new information on the nature of the underlying interactions in protein folding and chemical reactions and demonstrate the power of single‐molecule techniques to identify the changes induced by a small change in hydrogen bond strength.     

Molecular dynamics study of the hydration of the hydroxyl radical at body temperature

Classical molecular dynamics (MD) simulation of ˙OH in liquid water at 37 °C has been performed using flexible models of the solute and solvent molecules. We derived the Morse function describing the bond stretching of the radical and the potential for ˙OH-H(2)O interactions, including short-range interactions of hydrogen atoms. Scans of the potential energy surface of the ˙OH-H(2)O complex have been performed using the DFT method with the B3LYP functional and the 6-311G(d,p) basis set. The DFT-derived partial charges, ±0.375e, and the equilibrium bond-length, 0.975 ?, of ˙OH resulted in the dipole moment of 1.76 D. The radical-water radial distribution functions revealed that ˙OH is not built into the solvent structure but it rather occupies distortions or cavities in the hydrogen-bonded network. The solvent structure at 37 °C has been found to be the same as that of pure water. The hydration cage of the radical comprises 13-14 water molecules. The estimated hydration enthalpy -42 ± 5 kJ mol(-1) is comparable with the experimental value -39 ± 6 kJ mol(-1) for 25 °C. Inspection of hydrogen bonds showed the importance of short-range interaction of hydrogen atoms and indicated that neglect of the angular condition greatly overestimates the number of the H-acceptor radical-water bonds. The mean number ?n = 0.85 of radical-water H-bonds has been calculated using geometric definition of H-bond and ?n = 0.62 has been obtained when the energetic condition, E(da)≤-8 kJ mol(-1), was additionally considered. The continuous lifetimes of 0.033 ps and 0.024 ps have been estimated for the radical H-donor and the H-acceptor bonds, respectively. Within statistical uncertainty the radical self-diffusion coefficient, (2.9 ± 0.6) × 10(-9) m(2) s(-1), is the same as (3.1 ± 0.5) × 10(-9) m(2) s(-1) calculated for water in solution and in pure solvent. To the best of our knowledge, this is the first study of the ˙OH(aq) properties at a biologically relevant body temperature.