Estimating the temperature dependence of peptide folding entropies and free enthalpies from total energies in molecular dynamics simulations

The temperature dependence of thermodynamic quantities, such as heat capacity, entropy and free enthalpy, may be obtained by using well-known equations that relate these quantities to the enthalpy of the molecular system of interest at a range of temperatures. In turn, the enthalpy of a molecular system can be estimated from molecular dynamics simulations of an appropriate model. To demonstrate this, we have investigated the temperature dependence of the enthalpy, heat capacity, entropy and free enthalpy of a system that consists of a beta-heptapeptide in methanol and have used the statistical mechanics relationships to describe the thermodynamics of the folding/unfolding equilibrium of the peptide. The results illustrate the power of current molecular simulation force fields and techniques in establishing the link between thermodynamic quantities and conformational distributions.     

温控分子动力学研究微管蛋白活性肽链的折叠机制

运用温控和常温分子动力学方法, 研究了微管蛋白活性部位Pep1-28肽链的折叠机制, 总模拟时间为380.0 ns. 对于温控分子动力学, 逐渐降温可以清晰显示Pep1-28肽链的折叠途径, 发生明显折叠的温度约为550 K, 其折叠和去折叠可逆机制为U(>1200 K)←→I1(1200-1000 K)←→I2(800 K)←→I3(600 K)←→I4(450 K)←→F1(400 K)←→F2(300 K), 其中U为去折叠态构象, I1、I2、I3和I4是折叠过程中的四个重要的中间态构象, F1和F2是两个结构相近的折叠态构象. 对于常温(300 K)分子动力学, 其构象转变和折叠过程相当迅速, 很难观察到有效、稳定的中间态构象. 尤其引人注意的是, 其折叠态结构陷入了能量局域极小点, 与温控(300 K)的有较大差异, 两者能量差高达297.53 kJ·mol-1. 可见, 温控分子动力学方法不仅清晰地显示多肽和蛋白质折叠过程的重要中间态构象, 为折叠和去折叠机制提供直接、可靠的依据, 而且还有助于跨越较高的构象能垒, 促使多肽和蛋白质折叠以形成全局能量最低的稳定结构.     

Precision of free energies calculated by molecular dynamics simulations of peptides in solution

Errors in free energies for molecular replacement and for conformation change of a small model peptide have been determined empirically by repeated simulations from different starting points. All calculations have been done using thermodynamic integration, in which the system's potential energy is coupled to a parameter λ, that is increased or decreased by a small amount at each step of the simulation. The effects of several factors that may alter the precision are evaluated. These factors include: the length of the simulation, the dependence of the potential energy on λ, the use of conformational restraints, and their magnitude and form. The methods used for restraint and conformational forcing are described in detail. The free energy change, calculated as the mean from several successive simulations with alternately increasing and decreasing λ, is found to be independent of the length of the simulations. As expected, longer simulations produce more precise results. The variation of the calculated free energies is found to consist of two parts, a random error and a systematic hysteresis, i.e., a dependence on the direction in which λ changes. The hysteresis varies as the inverse of the length of the simulation and the random error as the inverse square root The advantage of the use of a different (nonlinear) dependence of the attractive and repulsive parts of the nonbonded potential energy on the coupling parameter when “creating” particles in solution is found to be very large. This nonlinear coupling was found to be superior to the use of linear coupling and a nonlinear change of the coupling parameter with the simulation time. The hysteresis in conformational free energy calculations is found to increase markedly if too weak a forcing restraint is chosen. It is shown how to deconvolute the contribution of a torsional restraint from the dependence of the free energy on a torsion angle.     

Folding of the 25 residue Abeta(12-36) peptide in TFE/water: temperature-dependent transition from a funneled free-energy landscape to a rugged one

The free-energy landscape of the Alzheimer beta-amyloid peptide Abeta(12-36) in a 40% (v/v) 2,2,2-trifluoroethanol (TFE)/water solution was determined by using multicanonical molecular dynamics simulations. Simulations using this enhanced conformational sampling technique were initiated from a random unfolded polypeptide conformation. Our simulations reliably folded the peptide to the experimental NMR structure, which consists of two linked helices. The shape of the free energy landscape for folding was found to be strongly dependent on temperature: Above 325 K, the overall shape was funnel-like, with the bottom of the funnel coinciding exactly with the NMR structure. Below 325 K, on the other hand, the landscape became increasingly rugged, with the emergence of new conformational clusters connected by low free-energy pathways. Finally, our simulations reveal that water and TFE solvate the polypeptide in different ways: The hydrogen bond formation between TFE and Abeta was enhanced with decreasing temperature, while that between water and Abeta was depressed.     

Conformational evolution of ubiquitin ions in electrospray mass spectrometry: molecular dynamics simulations at gradually increasing temperatures

Evidence from cross section data indicates that ubiquitin +13 ions lose their secondary and tertiary structure in mass spectrometric experiments. These transitions from the folded state into the near linear final structure occur at the experimental temperatures on time scales that are far too long for conventional molecular dynamics simulations. In this study, an approach to mass spectrometric unfolding processes is developed and a detailed application to an ubiquitin +13 ion system is presented. The approach involves a sequence of molecular dynamics simulations at gradually increasing temperatures leading to identification of major intermediate states, and the unfolding pathway. The unfolding rate at any temperature can then be calculated by a Rice-Ramsperger-Kassel (RRK) approach. For ubiquitin +13, three interesting intermediate states were found and the final near linear geometry was computed. The several relevant energy barriers calculated for the process are in the range of 7 to 15 kcal mol(-1). The unfolding time scale at 300 K was computed to be 2 ms. Cross section calculations using a hard sphere scattering model were carried out for the final structure and found to be in good accord with the results of electrospray experiments supporting the theoretical model used. The approach employed here should be applicable to any other solvent-free protein system.     

Temperature dependence of dimerization and dewetting of large-scale hydrophobes: a molecular dynamics study

We studied by molecular dynamics simulations the temperature dependence of hydrophobic association and drying transition of large-scale solutes. Similar to the behavior of small solutes, we found the association process to be characterized by a large negative heat capacity change. The origin of this large change in heat capacity is the high fragility of hydrogen bonds between water molecules at the interface with hydrophobic solutes; an increase in temperature breaks more hydrogen bonds at the interface than in the bulk. With increasing temperature, both entropy and enthalpy changes for association strongly decrease, while the change in free energy weakly varies, exhibiting a small minimum at high temperatures. At around T=Ts=360 K, the change in entropy is zero, a behavior similar to the solvation of small nonpolar solutes. Unexpectedly, we find that at Ts, there is still a substantial orientational ordering of the interfacial water molecules relative to the bulk. Nevertheless, at this point, the change in entropy vanishes due to a compensating contribution of translational entropy. Thus, at Ts, there is rotational order and translational disorder of the interfacial water relative to bulk water. In addition, we studied the temperature dependence of the drying-wetting transition. By calculating the contact angle of water on the hydrophobic surface at different temperatures, we compared the critical distance observed in the simulations with the critical distance predicted by macroscopic theory. Although the deviations of the predicted from the observed values are very small (8-23%), there seems to be an increase in the deviations with an increase in temperature. We suggest that these deviations emerge due to increased fluctuations, characterizing finite systems, as the temperature increases.     

Inverse temperature transition of a biomimetic elastin model: reactive flux analysis of folding/unfolding and its coupling to solvent dielectric relaxation

The inverse temperature transition (ITT) of a biomimetic model for elastin, capped GVG(VPGVG) in liquid water, is investigated by a comprehensive classical molecular dynamics study. The temperature dependence of the solvation structure and dynamics of the octapeptide are compared using three common force fields, CHARMM, GROMOS, and OPLS. While these force fields differ in quantitative detail, they all predict this octapeptide to undergo a "folding transition" to closed conformations upon heating and a subsequent "unfolding transition" to open conformations at still higher temperatures, thus reproducing the ITT scenario. The peptide kinetics is analyzed within the reactive flux formalism applied to the largest-amplitude mode extracted from principal component analysis, and the solvent's dielectric fluctuations are obtained from the total water dipole autocorrelations. Most importantly, preliminary evidence for an intimate coupling of peptide folding/unfolding dynamics, and thus the ITT, and dielectric relaxation of bulk water is given, possibly being consistent with a "slave mode" picture.     

Temperature dependence of NMR order parameters and protein dynamics

The helical subdomain, HP36, of the F-actin-binding headpiece domain of chicken villin, is the smallest naturally occurring polypeptide that folds to a thermostable compact structure. Unconstrained molecular dynamics simulations and constrained molecular dynamics simulations using umbrella sampling are used to study the temperature dependence of internal motions of the backbone amide moieties of HP36. The potential of mean force (PMF) for the N-H bond vector, determined from the constrained simulations, is found to be temperature dependent. A simple analytical expression is derived that describes the temperature dependence of the PMF. The parameters of this model are obtained from the PMF, from the unconstrained molecular dynamics simulations, or from experimental values of the generalized order parameter. The results provide a linkage between experimental and theoretical measures of the temperature dependence of protein motions.     

Molecular dynamics study on hydrophobic effects in aqueous urea solutions.

Hydrophobic effects in aqueous urea were analyzed by molecular dynamics simulations. The contribution of solvents to the potential of mean force between two methane molecules was calculated by using molecular dynamics simulations and was compared with the solubility data of hydrocarbons in aqueous urea. Both the simulation results and the solubility data indicated that urea stabilizes methane-methane association. The stabilization was due to increasing the solvation free energies of small hydrocarbons such as methane by addition of urea. The solvation free energies of larger hydrocarbons, on the other hand, are decreased by addition of urea. This effect of the solute size on hydrophobic free energies in aqueous urea was also analyzed by using molecular dynamics simulations by means of division of the solvation process into two parts: the cavity formation and the introduction of the solute-solvent attractive interactions. In the cavity formation, urea increased hydrophobic free energies, and in the introduction of the solute-solvent attractive interactions, urea decreased hydrophobic free energies. The influence of urea on hydrophobic free energies was determined by the balance of effects of the two parts of the solvation process.     

Simulations of lipid crystals: Characterization of potential energy functions and parameters for lecithin molecules

We evaluate an empirical potential energy function and associated parameters for classical molecular dynamics simulations of lecithins, a common class of lipid. The physical accuracy of the force field was tested through its application to molecular dynamics simulations of the known crystal structures of lipid molecules. Average atomic positions and molecular conformation are well maintained during the simulations despite considerable thermal motion. Calculated isotropic temperature factors correlate highly with those from experiment.     

Ionic melts with waterlike anomalies: thermodynamic properties of liquid BeF2

Thermodynamic properties of liquid beryllium difluoride (BeF(2)) are studied using canonical ensemble molecular dynamics simulations of the transferable rigid ion model potential. The negative slope of the locus of points of maximum density in the temperature-pressure plane is mapped out. The excess entropy, computed within the pair correlation approximation, is found to show an anomalous increase with isothermal compression at low temperatures which will lead to diffusional as well as structural anomalies resembling those in water. The anomalous behavior of the entropy is largely connected with the behavior of the Be-F pair correlation function. The internal energy shows a T(35) temperature dependence. The pair correlation entropy shows a T(-25) temperature dependence only at high densities and temperatures. The correlation plots between internal energy and the pair correlation entropy for isothermal compression show the characteristic features expected of network-forming liquids with waterlike anomalies. The tagged particle potential energy distributions are shown to have a multimodal form at low temperatures and densities similar to those seen in other liquids with three-dimensional tetrahedral networks, such as water and silica.     

Simple continuous and discrete models for simulating replica exchange simulations of protein folding

The efficiency of temperature replica exchange (RE) simulations hinge on their ability to enhance conformational sampling at physiological temperatures by taking advantage of more rapid conformational interconversions at higher temperatures. While temperature RE is a parallel simulation technique that is relatively straightforward to implement, kinetics in the RE ensemble is complicated, and there is much to learn about how best to employ RE simulations in computational biophysics. Protein folding rates often slow down above a certain temperature due to entropic bottlenecks. This "anti-Arrhenius" behavior represents a challenge for RE. However, it is far from straightforward to systematically explore the impact of this on RE by brute force molecular simulations, since RE simulations of protein folding are very difficult to converge. To understand some of the basic mechanisms that determine the efficiency of RE, it is useful to study simplified low dimensionality systems that share some of the key characteristics of molecular systems. Results are presented concerning the efficiency of temperature RE on a continuous two-dimensional potential that contains an entropic bottleneck. Optimal efficiency was obtained when the temperatures of the replicas did not exceed the temperature at which the harmonic mean of the folding and unfolding rates is maximized. This confirms a result we previously obtained using a discrete network model of RE. Comparison of the efficiencies obtained using the continuous and discrete models makes it possible to identify non-Markovian effects, which slow down equilibration of the RE ensemble on the more complex continuous potential. In particular, the rate of temperature diffusion and also the efficiency of RE is limited by the time scale of conformational rearrangements within free energy basins.     

Comparison of kinetic Monte Carlo and molecular dynamics simulations of diffusion in a model glass former

We present results from kinetic Monte Carlo (KMC) simulations of diffusion in a model glass former. We find that the diffusion constants obtained from KMC simulations have Arrhenius temperature dependence, while the correct behavior, obtained from molecular dynamics simulations, can be super-Arrhenius. We conclude that the discrepancy is due to undersampling of higher-lying local minima in the KMC runs. We suggest that the relevant connectivity of minima on the potential energy surface is proportional to the energy density of the local minima, which determines the "inherent structure entropy." The changing connectivity with potential energy may produce a correlation between dynamics and thermodynamics.     

Towards an understanding of the molecular basis of hydrophobicity

Summary The modern view is stressed that the structuring of water around nonpolar solutes, a process called hydrophobic hydration, actually favors the solubility of nonpolar solutes in water, its associated positive free energy of transfer arising from the enthalpic input required to create a cavity in water to accommodate the solute. The results of a series of molecular dynamics simulations of methane in SPC/E water at different temperatures are reported. These results show the existence of a larger fraction of broken hydrogen bonds in the hydration-shell water of the nonpolar solutes with respect to the bulk water, the difference increasing with a rise in temperature. This supports Muller's modified hydration-shell hydrogen-bond model predictions, where hydration-shell water molecules have lower free energies of hydrogen-bond breaking than those in the bulk.This paper is based on a presentation given at the 14th Molecular Graphics and Modelling Society Conference, held in Cairns, Australia, August 27 September 1, 1995.     

Effects of temperature and pH on the helicity of a peptide adsorbed to colloidal silica

The conformation of a cationic -helical peptide (DDDDAAAARRRRR) adsorbed to anionic colloidal silica has been investigated by circular dichroism (CD) spectroscopy as a function of temperature and pH in order to examine how the structure of an adsorbed molecule responds to two simultaneous perturbations. Increased temperature destabilizes the helicity of the peptide in solution, while pH changes alter the substrate surface charge and the corresponding strength of the interaction with the peptide. Near neutral pH, the helicity of the adsorbed peptide, which is determined from the intensity of the CD signal at 222 nm, decreases with increasing temperature, similarly to the temperature-dependent behavior observed for the peptide in aqueous solution. By contrast, at basic pH and a strongly negative surface charge, the helicity of the adsorbed peptide increases with temperature. In order to elucidate the origin of the reversal of the temperature dependence of helicity, a statistical model for the conformation of the adsorbed peptide has been formulated based on the Lifson–Roig model for the helix–coil transition of the peptide in solution. The model provides insight into the trends in fractional helicity and reveals that the temperature dependence of the helicity of the adsorbed peptide results from a competition between the intramolecular interactions that promote helicity and the intermolecular interactions with the surface. The statistical model also enables estimation of the free energy contributions from specific aspects of the adsorption process. Through identification of a connection between the conformation of adsorbed peptide and the interactions of the peptide with the surface, this work suggests a route for the control of adsorbate conformation through peptide and surface engineering.     

多肽抑制剂抑制淀粉质多肽42构象转换的分子动力学模拟和结合自由能计算

应用分子动力学模拟和结合自由能计算方法研究了多肽抑制剂KLVFF、VVIA和LPFFD抑制淀粉质多肽42 (Aβ42)构象转换的分子机理. 结果表明, 三种多肽抑制剂均能够有效抑制Aβ42的二级结构由α-螺旋向β-折叠的构象转换. 另外, 多肽抑制剂降低了Aβ42分子内的疏水相互作用, 减少了多肽分子内远距离的接触, 有效抑制了Aβ42的疏水塌缩, 从而起到稳定其初始构象的作用. 这些抑制剂与Aβ42之间的疏水和静电相互作用(包括氢键)均有利于它们抑制Aβ42的构象转换. 此外, 抑制剂中的带电氨基酸残基可以增强其和Aβ42之间的静电相互作用(包括氢键), 并降低抑制剂之间的聚集, 从而大大增强对Aβ42构象转换的抑制能力. 但脯氨酸的引入会破坏多肽的线性结构, 从而大大降低其与Aβ42 之间的作用力. 上述分子模拟的结果揭示了多肽抑制剂KLVFF、VVIA和LPFFD抑制Aβ42构象转换的分子机理, 对于进一步合理设计Aβ的高效短肽抑制剂具有非常重要的理论指导意义.     

Calibration of force-field dependency in free energy landscapes of peptide conformations by quantum chemical calculations

The free energy landscapes of peptide conformations were calibrated by ab initio quantum chemical calculations, after the enhanced conformational diversity search using the multicanonical molecular dynamics simulations. Three different potentials of mean force for an isolated dipeptide were individually obtained by the multicanonical molecular dynamics simulations using the conventional force fields, AMBER parm94, AMBER parm96, and CHARMm22. Each potential of mean force was then calibrated based upon the umbrella sampling algorithm from the adiabatic energy map that was calculated separately by the ab initio molecular orbital method, and all of the calibrated potentials of mean force coincided well. The calibration method was also applied to the simulations of a peptide dimer in explicit water models, and it was shown that the calibrated free energy landscapes did not depend on the force field used in the classical simulations, as far as the conformational space was sampled well. The current calibration method fuses the classical free energy calculation with the quantum chemical calculation, and it should generally make simulations for biomolecular systems much more reliable when combining with enhanced conformational sampling.     

On the Debye-Waller factor of hexagonal ice: a computer simulation study

We investigate by molecular dynamics (MD) simulations the temperature dependence of the Debye-Waller (DW) factor of hexagonal ice with 25 different proton-disordered configurations. Each initial configuration is composed of 288 water molecules with no net dipole moment. The intermolecular interaction of water is described by TIP4P potential. Each production run of the simulation is 15 ns or longer. We observe a change in slope of the DW factor around 200 K, which cannot be explained within the framework of either classical or quantum harmonic approximation. Configurations generated by MD simulations are subjected to the steepest descent energy minimization. Analysis of the local energy minimum structures reveals that water molecules above 200 K jump to other lattice sites via some local energy minimum structures which contain some water molecules sitting on the locations other than the lattice sites. As time evolves, these defect molecules move back and forth to the lattice sites yielding defect-free structures. Those motions are responsible for the unusual increase in the DW factor at high temperatures. In making a transition from an energy-minimum structure to another one, a small number of water molecules are involved in a highly cooperative fashion. The larger DW factor at higher temperature arises from jump-like motions of water molecules among these locally stable configurations which may or may not be a family of the proton-disordered ice forms satisfying the "ice rule".