Accelerated molecular dynamics simulations of protein folding

Folding of four fast‐folding proteins, including chignolin, Trp‐cage, villin headpiece and WW domain, was simulated via accelerated molecular dynamics (aMD). In comparison with hundred‐of‐microsecond timescale conventional molecular dynamics (cMD) simulations performed on the Anton supercomputer, aMD captured complete folding of the four proteins in significantly shorter simulation time. The folded protein conformations were found within 0.2–2.1 Å of the native NMR or X‐ray crystal structures. Free energy profiles calculated through improved reweighting of the aMD simulations using cumulant expansion to the second‐order are in good agreement with those obtained from cMD simulations. This allows us to identify distinct conformational states (e.g., unfolded and intermediate) other than the native structure and the protein folding energy barriers. Detailed analysis of protein secondary structures and local key residue interactions provided important insights into the protein folding pathways. Furthermore, the selections of force fields and aMD simulation parameters are discussed in detail. Our work shows usefulness and accuracy of aMD in studying protein folding, providing basic references in using aMD in future protein‐folding studies. © 2015 Wiley Periodicals, Inc.     

Thermal analysis of interpenetrating polymer networks through molecular dynamics simulations: a comparison with experiments

Journal of Thermal Analysis and Calorimetry - In this work, we verified the synthesis of a novel sequential interpenetrating polymer network, composed of poly(2-hexyl-ethylacrylate) and...     

Application of the accelerated molecular dynamics simulations to the folding of a small protein

In this paper, we further explore the applicability of the accelerated molecular dynamics simulation method using a bias potential. The method is applied to both simple model systems and real multidimensional systems. The method is also compared to replica exchange simulations in folding a small protein, Trp cage, using an all atom potential for the protein and an implicit model for the solvent. We show that the bias potential method allows quick searches of folding pathways. We also show that the choice of the bias potential has significant influence on the efficiency of the bias potential method.     

Folding processes of the B domain of protein A to the native state observed in all-atom ab initio folding simulations

Reaching the native states of small proteins, a necessary step towards a comprehensive understanding of the folding mechanisms, has remained a tremendous challenge to ab initio protein folding simulations despite the extensive effort. In this work, the folding process of the B domain of protein A (BdpA) has been simulated by both conventional and replica exchange molecular dynamics using AMBER FF03 all-atom force field. Started from an extended chain, a total of 40 conventional (each to 1.0 micros) and two sets of replica exchange (each to 200.0 ns per replica) molecular dynamics simulations were performed with different generalized-Born solvation models and temperature control schemes. The improvements in both the force field and solvent model allowed successful simulations of the folding process to the native state as demonstrated by the 0.80 A C(alpha) root mean square deviation (RMSD) of the best folded structure. The most populated conformation was the native folded structure with a high population. This was a significant improvement over the 2.8 A C(alpha) RMSD of the best nativelike structures from previous ab initio folding studies on BdpA. To the best of our knowledge, our results demonstrate, for the first time, that ab initio simulations can reach the native state of BdpA. Consistent with experimental observations, including Phi-value analyses, formation of helix II/III hairpin was a crucial step that provides a template upon which helix I could form and the folding process could complete. Early formation of helix III was observed which is consistent with the experimental results of higher residual helical content of isolated helix III among the three helices. The calculated temperature-dependent profile and the melting temperature were in close agreement with the experimental results. The simulations further revealed that phenylalanine 31 may play critical to achieve the correct packing of the three helices which is consistent with the experimental observation. In addition to the mechanistic studies, an ab initio structure prediction was also conducted based on both the physical energy and a statistical potential. Based on the lowest physical energy, the predicted structure was 2.0 A C(alpha) RMSD away from the experimentally determined structure.     

Experimental resolution of early steps in protein folding: testing molecular dynamics simulations

Time-resolved Tyr fluorescence spectroscopy coupled with a laser-induced temperature-jump (T-jump) was employed to follow the folding relaxation dynamics of the B-domain of Staphylococcal protein A. The single Tyr is located in helix 1 (H1) and is a sensitive probe of the structure of this helix and the overall helical bundle structure. The results from this study were compared to those from a complementary infrared T-jump study on this protein [Vu, D. M.; Myers, J. K.; Oas, T. G.; Dyer, R. B. Biochemistry 2004, 43, 3582]. Both methods detect a microsecond process that follows the cooperative relaxation of the helical bundle core. However, a fast process (10-7 s) that follows the relaxation of the individual helices was observed only with the infrared probe. Thus, fast formation of H1 is not observed, but rather H1 forms in the microsecond phase, concomitantly with the docking to (and stabilization by) the other two helices to form the helical bundle structure. This observation validates the results of several previous molecular dynamics simulations that predict H1 formation only in the final assembly of the helix bundle.     

The Trp cage: folding kinetics and unfolded state topology via molecular dynamics simulations

Using over 75 mus of molecular dynamics simulation, we have generated several thousand folding simulations of the 20-residue Trp cage at experimental temperature and solvent viscosity. A total of 116 independent folding simulations reach RMSDcalpha values below 3 A RMSDcalpha, some as close as 1.4 A RMSDcalpha. We estimate a folding time of 5.5+/-3.5 mus, a rate that is in reasonable agreement with experimental kinetics. Finally, we characterize both the folded and unfolded ensemble under native conditions and note that the average topology of the unfolded ensemble is very similar to the topology of the native state.     

The ensemble folding dynamics of EF-hand domain in parvalbumin from a Monte Carlo simulation

The helix-loop-helix (i.e., EF-hand) is the most common motif in the superfamily of Ca\(^{2+}\)-binding proteins. Parvalbumins, as the classical EF-hand proteins, are associated with the muscle relaxation/contration, the calcium buffering etc. The previous researches focus more of interest on the ion coupled/decoupled mechanism using molecular dynamics (MD) computations. We developed the novel approach instead of MD, in which the Landau free energy was introduced to describe the protein in term of the skeletal C\(_\alpha \) chain. The unfolding and folding processes were simulated by the Glauber algorithm under the repeated heating and cooling cycles. The high-quality crystal structure (2PVB) was as the trigger of non-equilibrium dynamics simulation for parvalbumin-\(\beta \) (Parv). The evolution of three dimensional structure during the folding was displayed for the residues 8–64 fragment in Parv. The phenomenon of helix nucleation was observed, that was, the 3\(_{10}\) helix at residues 35–37 and the front of \(\alpha \)-helix at residues 40–50 were firstly formed in the folding process. We also found two potential misfoldings which were caused by the distortion of the local structures at residues 35–37, 45–50.     

The influence of aqueous versus glassy solvents on protein dynamics: vibrational echo experiments and molecular dynamics simulations

Spectrally resolved infrared stimulated vibrational echo measurements are used to measure the vibrational dephasing of the CO stretching mode of carbonmonoxy-hemoglobin (HbCO), a myoglobin mutant (H64V), and a bacterial cytochrome c(552) mutant (Ht-M61A) in aqueous solution and trehalose glasses. The vibrational dephasing of the heme-bound CO is significantly slower for all three proteins embedded in trehalose glasses compared to that of aqueous protein solutions. All three proteins exhibit persistent but notably slower spectral diffusion when the protein surface is fixed by the glassy solvent. Frequency-frequency correlation functions (FFCFs) of the CO are extracted from the vibrational echo data to reveal that the structural dynamics, as sensed by the CO, of the three proteins in trehalose and aqueous solution are dominated by fast (tens of femtoseconds), motionally narrowed fluctuations. MD simulations of H64V in dynamic and "static" water are presented as models of the aqueous and glassy environments. FFCFs are calculated from the H64V simulations and qualitatively reproduce the important features of the experimentally extracted FFCFs. The suppression of long time scale (picoseconds to tens of picoseconds) frequency fluctuations (spectral diffusion) in the glassy solvent is the result of a damping of atomic displacements throughout the protein structure and is not limited to structural dynamics that occur only at the protein surface. The analysis provides evidence that some dynamics are coupled to the hydration shell of water, supporting the idea that the bioprotection offered by trehalose is due to its ability to immobilize the protein surface through a thin, constrained layer of water.     

A comparison of methods for melting point calculation using molecular dynamics simulations

Accurate and efficient prediction of melting points for complex molecules is still a challenging task for molecular simulation, although many methods have been developed. Four melting point computational methods, including one free energy-based method (the pseudo-supercritical path (PSCP) method) and three direct methods (two interface-based methods and the voids method) were applied to argon and a widely studied ionic liquid 1-n-butyl-3-methylimidazolium chloride ([BMIM][Cl]). The performance of each method was compared systematically. All the methods under study reproduce the argon experimental melting point with reasonable accuracy. For [BMIM][Cl], the melting point was computed to be 320 K using a revised PSCP procedure, which agrees with the experimental value 337-339 K very well. However, large errors were observed in the computed results using the direct methods, suggesting that these methods are inappropriate for large molecules with sluggish dynamics. The strengths and weaknesses of each method are discussed.     

The free energy of the metastable supersaturated vapor via restricted ensemble simulations. II. Effects of constraints and comparison with molecular dynamics simulations

Extensive restricted canonical ensemble Monte Carlo simulations [D. S. Corti and P. Debenedetti, Chem. Eng. Sci. 49, 2717 (1994)] were performed. Pressure, excess chemical potential, and excess free energy with respect to ideal gas data were obtained at different densities of the supersaturated Lennard-Jones (LJ) vapor at reduced temperatures from 0.7 to 1.0. Among different constraints imposed on the system studied, the one with the local minimum of the excess free energy was taken to be the approximated equilibrium state of the metastable LJ vapor. Also, a comparison of our results with molecular dynamic simulations [A. Linhart et al., J. Chem. Phys. 122, 144506 (2005)] was made.     

Self-assembly of a catechol-based macrocycle at the liquid-solid interface: experiments and molecular dynamics simulations

This combined experimental (STM, XPS) and molecular dynamics simulation study highlights the complex and subtle interplay of solvent effects and surface interactions on the 2-D self-assembly pattern of a Schiff-base macrocycle containing catechol moieties at the liquid-solid interface. STM imaging reveals a hexagonal ordering of the macrocycles at the n-tetradecane/Au(111) interface, compatible with a desorption of the lateral chains of the macrocycle. Interestingly, all the triangular-shaped macrocycles are oriented in the same direction, avoiding a close-packed structure. XPS experiments indicate the presence of a strong macrocycle-surface interaction. Also, MD simulations reveal substantial solvent effects. In particular, we find that co-adsorption of solvent molecules with the macrocycles induces desorption of lateral chains, and the solvent molecules act as spacers stabilizing the open self-assembly pattern.     

Probing the transition state ensemble of a protein folding reaction by pressure-dependent NMR relaxation dispersion

The F61A/A90G mutant of a redesigned form of apocytochrome b562 folds by an apparent two-state mechanism. We have used the pressure dependence of 15N NMR relaxation dispersion rate profiles to study the changes in volumetric parameters that accompany the folding reaction of this protein at 45 degrees C. The experiments were performed under conditions where the folding/unfolding equilibrium could be studied at each pressure without addition of denaturants. The exquisite sensitivity of the methodology to small changes in folding/unfolding rates facilitated the use of relatively low-pressure values (between 1 and 270 bar) so that pressure-induced changes to the unfolded state ensemble could be minimized. A volume change for unfolding of -81 mL/mol is measured (at 1 bar), a factor of 1.4 larger (in absolute value) than the volume difference between the transition state ensemble (TSE) and the unfolded state. Notably, the changes in the free energy difference between folded and unfolded states and in the activation free energy for folding were not linear with pressure. Thus, the difference in the isothermal compressibility upon unfolding (-0.11 mL mol(-1) bar(-1)) and, for the first time, the compressibility of the TSE relative to the unfolded state (0.15 mL mol(-1) bar(-1)) could be calculated. The results argue for a TSE that is collapsed but loosely packed relative to the folded state and significantly hydrated, suggesting that the release of water occurs after the rate-limiting step in protein folding. The notion of a collapsed and hydrated TSE is consistent with expectations based on earlier temperature-dependent folding studies, showing that the barrier to folding at 45 degrees C is entropic (Choy, W. Y.; Zhou, Z.; Bai, Y.; Kay, L. E. J. Am. Chem. Soc. 2005, 127, 5066-5072).     

A detailed reaction mechanism for hexamethyldisiloxane combustion via experiments and ReaxFF molecular dynamics simulations

Hexamethyldisiloxane (HMDSO) is one of the main impurities in the syngas produced from sewage and landfill plants. In order to utilize this syngas or control the characteristics of the generated silica particles, it is crucial to understand the chemical kinetics of HMDSO combustion. This study investigated the process of HMDSO combustion using synchrotron radiation mass spectrometry (SRMS), gas chromatography (GC), and ReaxFF molecular dynamics simulations. First, the force field used for ReaxFF simulation was validated by comparing the energies of different bond lengths, bond angles, and dihedral angles with the ones from DFT calculations. Good agreements were found. Then, ReaxFF simulations of HMDSO combustion with this force field were conducted under various conditions, which include different equivalence ratios (0.67, 1.0, and 1.5) and temperatures ranging from 2000 to 3500 K. The oxidation characteristics of HMDSO were analyzed, including the evolution of gas products and particle formation. Finally, based on the results from experiments and ReaxFF simulations, the reaction pathways, reaction lists, and reaction kinetics data during HMDSO combustion were obtained. A detailed reaction mechanism was proposed and validated by applying it in modeling the H₂/HMDSO/O₂ combustion systems. The temperature and part of the gas products such as CO and CO₂ as well as SiO could be well predicted.     

Adsorption of acetic acid on ice: experiments and molecular dynamics simulations

Adsorption study of acetic acid on ice surfaces was performed by combining experimental and theoretical approaches. The experiments were conducted between 193 and 223 K using a coated wall flow tube coupled to a mass spectrometric detection. Under our experimental conditions, acetic acid was mainly dimerized in the gas phase. The surface coverage increases with decreasing temperature and with increasing concentrations of acetic acid dimers. The obtained experimental surface coverages were fitted according to the BET theory in order to determine the enthalpy of adsorption deltaH(ads) and the mololayer capacity N(M(dimers)) of the acetic acid dimers on ice: deltaH(ads) = (-33.5 +/- 4.2) kJ mol(-1), N(M(dimers)) = (l1.27 +/- 0.25) x 10(14) dimers cm(-2). The adsorption characteristics of acetic acid on an ideal ice I(n)(0001) surface were also studied by means of classical molecular dynamics simulations in the same temperature range. The monolayer capacity, the configurations of the molecules in their adsorption sites, and the corresponding adsorption energies have been determined for both acetic acid monomers and dimers, and compared to the corresponding data obtained from the experiments. In addition, the theoretical results show that the interaction with the ice surface could be strong enough to break the acetic acid dimers that exist in the gas phase and leads to the stabilization of acetic acid monomers on ice.     

Solvation dynamics in protein environments: comparison of fluorescence upconversion measurements of coumarin 153 in monomeric hemeproteins with molecular dynamics simulations

The complexes of the fluorescence probe coumarin 153 with apomyoglobin and apoleghemoglobin are used as model systems to study solvation dynamics in proteins. Time-resolved Stokes shift experiments are compared with molecular dynamics simulations, and very good agreement is obtained. The solvation of the coumarin probe is very rapid with approximately 60% occurring within 300 fs and is attributed to interactions with water (or possibly to the protein itself). Differences in the solvation relaxation (or correlation) function C(t) for the two proteins are attributed to differences in their hemepockets.     

Ground state structures and excited state dynamics of pyrrole-water complexes: ab initio excited state molecular dynamics simulations

Structures of the ground state pyrrole-(H2O)n clusters are investigated using ab initio calculations. The charge-transfer driven femtosecond scale dynamics are studied with excited state ab initio molecular dynamics simulations employing the complete-active-space self-consistent-field method for pyrrole-(H2O)n clusters. Upon the excitation of these clusters, the charge density is located over the farthest water molecule which is repelled by the depleted pi-electron cloud of pyrrole ring, resulting in a highly polarized complex. For pyrrole-(H2O), the charge transfer is maximized (up to 0.34 a.u.) around approximately 100 fs and then oscillates. For pyrrole-(H2O)2, the initial charge transfer occurs through the space between the pyrrole and the pi H-bonded water molecule and then the charge transfer takes place from this water molecule to the sigma H-bonded water molecule. The total charge transfer from the pyrrole to the water molecules is maximized (up to 0.53 a.u.) around approximately 100 fs.     

MDAnalysis: A toolkit for the analysis of molecular dynamics simulations

MDAnalysis is an object‐oriented library for structural and temporal analysis of molecular dynamics (MD) simulation trajectories and individual protein structures. It is written in the Python language with some performance‐critical code in C. It uses the powerful NumPy package to expose trajectory data as fast and efficient NumPy arrays. It has been tested on systems of millions of particles. Many common file formats of simulation packages including CHARMM, Gromacs, Amber, and NAMD and the Protein Data Bank format can be read and written. Atoms can be selected with a syntax similar to CHARMM's powerful selection commands. MDAnalysis enables both novice and experienced programmers to rapidly write their own analytical tools and access data stored in trajectories in an easily accessible manner that facilitates interactive explorative analysis. MDAnalysis has been tested on and works for most Unix‐based platforms such as Linux and Mac OS X. It is freely available under the GNU General Public License from http://mdanalysis.googlecode.com . © 2011 Wiley Periodicals, Inc. J Comput Chem 2011     

Dual control cell reaction ensemble molecular dynamics: a method for simulations of reactions and adsorption in porous materials

We present a simulation tool to study fluid mixtures that are simultaneously chemically reacting and adsorbing in a porous material. The method is a combination of the reaction ensemble Monte Carlo method and the dual control volume grand canonical molecular dynamics technique. The method, termed the dual control cell reaction ensemble molecular dynamics method, allows for the calculation of both equilibrium and nonequilibrium transport properties in porous materials such as diffusion coefficients, permeability, and mass flux. Control cells, which are in direct physical contact with the porous solid, are used to maintain the desired reaction and flow conditions for the system. The simulation setup closely mimics an actual experimental system in which the thermodynamic and flow parameters are precisely controlled. We present an application of the method to the dry reforming of methane reaction within a nanoscale reactor model in the presence of a semipermeable membrane that was modeled as a porous material similar to silicalite. We studied the effects of the membrane structure and porosity on the reaction species permeability by considering three different membrane models. We also studied the effects of an imposed pressure gradient across the membrane on the mass flux of the reaction species. Conversion of syngas (H2/CO) increased significantly in all the nanoscale membrane reactor models considered. A brief discussion of further potential applications is also presented.     

A molecular dynamics study of the correlations between solvent-accessible surface, molecular volume, and folding state

We analyzed the correlations between molecular volume, solvent-accessible surface, and folding state (secondary structure content) for unfolded conformers of alpha (holo- and apomyoglobin) and beta (retinal-binding protein) proteins and a small water-soluble alanine-rich alpha-helical peptide. Conformers with different degrees of folding were obtained using molecular dynamics at constant temperature and pressure with implicit solvent (dielectric constant adjustment) for all four systems and with explicit solvent for the single helix peptide. Our results support the view that unfolded conformations are not necessary extended, that volume variation is not a good indication of folding state and that the simple model of water penetrating the interior of the protein does not explain the increase in volume upon unfolding.