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.     

Coarse-grained protein molecular dynamics simulations

A limiting factor in biological science is the time-scale gap between experimental and computational trajectories. At this point, all-atom explicit solvent molecular dynamics (MD) are clearly too expensive to explore long-range protein motions and extract accurate thermodynamics of proteins in isolated or multimeric forms. To reach the appropriate time scale, we must then resort to coarse graining. Here we couple the coarse-grained OPEP model, which has already been used with activated methods, to MD simulations. Two test cases are studied: the stability of three proteins around their experimental structures and the aggregation mechanisms of the Alzheimer's Abeta16-22 peptides. We find that coarse-grained isolated proteins are stable at room temperature within 50 ns time scale. Based on two 220 ns trajectories starting from disordered chains, we find that four Abeta16-22 peptides can form a three-stranded beta sheet. We also demonstrate that the reptation move of one chain over the others, first observed using the activation-relaxation technique, is a kinetically important mechanism during aggregation. These results show that MD-OPEP is a particularly appropriate tool to study qualitatively the dynamics of long biological processes and the thermodynamics of molecular assemblies.     

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.     

Monte Carlo vs molecular dynamics for all-atom polypeptide folding simulations

An efficient Monte Carlo (MC) algorithm including concerted rotations is directly compared to molecular dynamics (MD) in all-atom statistical mechanics folding simulations of small polypeptides. The previously reported algorithm "concerted rotations with flexible bond angles" (CRA) has been shown to successfully locate the native state of small polypeptides. In this study, the folding of three small polypeptides (trpzip2/H1/Trp-cage) is investigated using MC and MD, for a combined sampling time of approximately 10(11) MC configurations and 8 micros, respectively. Both methods successfully locate the experimentally determined native states of the three systems, but they do so at different speed, with 2-2.5 times faster folding of the MC runs. The comparison reveals that thermodynamic and dynamic properties can reliably be obtained by both and that results from folding simulations do not depend on the algorithm used. Similar to previous comparisons of MC and MD, it is found that one MD integration step of 2 fs corresponds to one MC scan, revealing the good sampling of MC. The simplicity and efficiency of the MC method will enable its future use in folding studies involving larger systems and the combination with replica exchange algorithms.     

Ab initio folding of helix bundle proteins using molecular dynamics simulations

We have demonstrated that ab initio fast folding simulations at 400 K using a GB implicit solvent model with an all-atom based force field can describe the spontaneous formation of nativelike structures for the 36-residue villin headpiece and the 46-residue fragment B of Staphylococcal protein A. An implicit solvent model combined with high-temperature MD makes it possible to perform direct folding simulations of small- to medium-sized proteins by reducing the computational requirements tremendously. In the early stage of folding of the villin headpiece and protein A, initial hydrophobic collapse and rapid formation of helices were found to play important roles. For protein A, the third helix forms first in the early stage of folding and exhibits higher stability. The free energy profiles calculated from the folding simulations suggested that both of the helix-bundle proteins show a two-state thermodynamic behavior and protein A exhibits rather broad native basins.     

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.     

Modern protein force fields behave comparably in molecular dynamics simulations

Several molecular dynamics simulations were performed on three proteins--bovine apo-calbindin D9K, human interleukin-4 R88Q mutant, and domain IIA of bacillus subtilis glucose permease--with each of the AMBER94, CHARMM22, and OPLS-AA force fields as implemented in CHARMM. Structural and dynamic properties such as solvent-accessible surface area, radius of gyration, deviation from their respective experimental structures, secondary structure, and backbone order parameters are obtained from each of the 2-ns simulations for the purpose of comparing the protein portions of these force fields. For one of the proteins, the interleukin-4 mutant, two independent simulations were performed using the CHARMM22 force field to gauge the sensitivity of some of these properties to the specific trajectory. In general, the force fields tested performed remarkably similarly with differences on the order of those found for the two independent trajectories of interleukin-4 with CHARMM22. When all three proteins are considered together, no force field showed any consistent trend in variations for most of the properties monitored in the study.     

Solvent effect on the folding dynamics and structure of E6-associated protein characterized from ab initio protein folding simulations

Solvent effect on protein conformation and folding mechanism of E6-associated protein (E6ap) peptide are investigated using a recently developed charge update scheme termed as adaptive hydrogen bond-specific charge (AHBC). On the basis of the close agreement between the calculated helix contents from AHBC simulations and experimental results, we observed based on the presented simulations that the two ends of the peptide may simultaneously take part in the formation of the helical structure at the early stage of folding and finally merge to form a helix with lowest backbone RMSD of about 0.9 A? in 40% 2,2,2-trifluoroethanol solution. However, in pure water, the folding may start at the center of the peptide sequence instead of at the two opposite ends. The analysis of the free energy landscape indicates that the solvent may determine the folding clusters of E6ap, which subsequently leads to the different final folded structure. The current study demonstrates new insight to the role of solvent in the determination of protein structure and folding dynamics.     

Phonostat: thermostatting phonons in molecular dynamics simulations

Thermostat algorithms in a molecular dynamics simulation maintain an average temperature of a system by regulating the atomic velocities rather than the internal degrees of freedom. Herein, we present a "phonostat" algorithm that can regulate the total energy in a given internal degree of freedom. In this algorithm, the modal energies are computed at each time step using a mode-tracking scheme and then the system is driven by an external driving force of desired frequency and amplitude. The rate and amount of energy exchange between the phonostat and the system is controlled by two distinct damping parameters. Two different schemes for controlling the external driving force amplitude are also presented. In order to test our algorithm, the method is applied initially to a simple anharmonic oscillator for which the role of various phonostat parameters can be carefully tested. We then apply the phonostat to a more realistic (10,0) carbon nanotube system and show how such an approach can be used to regulate energy of highly anharmonic modes.     

Unbinding of nicotine from the acetylcholine binding protein: steered molecular dynamics simulations

The ligand binding/unbinding process is critical to our understanding of the pharmacology of both the nicotinic acetylcholine receptor (nAChR) and the acetylcholine binding protein (AChBP). Steered molecular dynamics simulations were performed to learn about the unbinding process of the full agonist nicotine. Three different pulling models were designed to investigate the possible binding/unbinding pathways: radial and tangent models, and also a mixed model. Of the three, the tangent pulling model finally failed to dissociate nicotine from the ligand binding pocket. The efficiency of the pulling force profiles was superior, and the opening of the C-loop was smaller in the mixed pulling model than that in the radial model. The most favorable pathway for the cholinergic agonist nicotine to enter or leave the binding pocket is through the principal binding side, following a curvilinear track. Noticeably, it has been seen that the unbinding of the nicotine is concomitant with a global rotation of the protein-ligand complex which could be caused by the interactions of the ligand with protein at the tangent direction.     

Cyan fluorescent protein: molecular dynamics, simulations, and electronic absorption spectrum

The dynamics and electronic absorption spectrum of enhanced cyan fluorescent protein (ECFP), a mutant of green fluorescent protein (GFP), have been studied by means of a 1 ns molecular dynamics (MD) simulation. The two X-ray conformations A' and B' of ECFP were considered. The chromophore was assumed to be neutral, and all titratable residues were taken in their standard protonation state at neutral pH. The protein was embedded in a box of water molecules (and counterions). The first result is that the two conformations A' and B' are found to be stable all along the simulation. Then, an analysis of the hydrogen-bond networks shows strong differences between the two conformations in the surroundings of the nitrogen atom of the indolic part of the chromophore. This is partly due to the imperfection in the beta barrel near the His148 residue, which allows the access of one solvent molecule inside the protein in conformation A'. Finally, quantum mechanical calculations of the electronic transition energies of the chromophore in the charge cloud of the protein and solvent water molecules were performed using the TDDFT method on 160 snapshots extracted every 5 ps of the MD trajectories. It is found that conformations A' and B' exhibit very similar spectra despite different H-bond networks involving the chromophore. This similarity is related to the weak charge transfer involved in the electronic transition and the weak electrostatic field created by ECFP near the chromophore, within the hypotheses made in the present simulation.     

Multibaric-multithermal ensemble molecular dynamics simulations

We present new generalized-ensemble molecular dynamics simulation algorithms, which we refer to as the multibaric-multithermal molecular dynamics. We describe three algorithms based on (1) the Nosé thermostat and the Andersen barostat, (2) the Nosé-Poincaré thermostat and the Andersen barostat, and (3) the Gaussian thermostat and the Andersen barostat. The multibaric-multithermal simulations perform random walks widely both in the potential-energy space and in the volume space. Therefore, one can calculate isobaric-isothermal ensemble averages in wide ranges of temperature and pressure from only one simulation run. We test the effectiveness of the multibaric-multithermal algorithm by applying it to a Lennard-Jones 12-6 potential system.     

Basin hopping simulations for all-atom protein folding

We investigate different protocols of the basin hopping technique for de novo protein folding. Using the protein free-energy force field PFF01 we report the reproducible all-atom folding of the 20-amino-acid tryptophan-cage protein [Protein Data Bank (PDB) code: 112y] and of the recently discovered 26-amino-acid potassium channel blocker (PDB code: 1wqc), which exhibits an unusual fold. We find that simulations with increasing cycle length and random starting temperatures perform best in comparison with other parametrizations. The basin hopping technique emerges as a simple but very efficient and robust workhorse for all-atom protein folding.     

Kinetic coefficient of steps at the Si(111) crystal-melt interface from molecular dynamics simulations

Nonequilibrium molecular dynamics simulations are applied to the investigation of step-flow kinetics at crystal-melt interfaces of silicon, modeled with the Stillinger-Weber potential [Phys. Rev. B 31, 5262 (1985)]. Step kinetic coefficients are calculated from crystallization rates of interfaces that are vicinals of the faceted (111) orientation. These vicinal interfaces contain periodic arrays of bilayer steps, and they are observed to crystallize in a step-flow growth mode at undercoolings lower than 40 K. Kinetic coefficients for both [110] and [121] oriented steps are determined for several values of the average step separation, in the range of 7.7-62.4 A. The values of the step kinetic coefficients are shown to be highly isotropic, and are found to increase with increasing step separation until they saturate at step separations larger than approximately 50 A. The largest step kinetic coefficients are found to be in the range of 0.7-0.8 m(sK), values that are more than five times larger than the kinetic coefficient for the rough (100) crystal-melt interface in the same system. The dependence of step mobility on step separation and the relatively large value of the step kinetic coefficient are discussed in terms of available theoretical models for crystal growth kinetics from the melt.     

Alamethicin channels in a membrane: molecular dynamics simulations

Alamethicin (Alm) is a 20 residue peptide which forms a kinked alpha-helix in membrane and membrane-mimetic environments. Ion channels formed by intramembraneous aggregates of Alm are thought to be formed by bundles of approximately parallel Alm helices surrounding a central bilayer pore. Different channel conductance levels correspond to different numbers of helices per bundle, ranging from N = 5 to N > 8. Calculation of the predicted pKA values of the ring of Glu18 sidechains at the C-terminal mouth of the pore suggests that at neutral pH most or all of these sidechains will remain protonated. Nanosecond molecular dynamics (MD) simulations of N = 5, 6, 7 and 8 bundles of Alm helices in a POPC bilayer have been run, corresponding to a total simulation time of 4 ns. These simulations explore the stability and conformational dynamics of these helix bundle channels when embedded in a full phospholipid bilayer in an aqueous environment. The structural and dynamic properties of water in these model channels are examined. As in earlier in vacuo simulations (J. Breed, R. Sankararamakrishnan, I. D. Kerr and M. S. P. Sansom, Biophys. J., 1996, 70, 1643) the dipole moments of water molecules within the pores are aligned antiparallel to the helix dipoles. This helps to contribute to the stability of the helix bundles.     

On simulations of complex interfaces: molecular dynamics simulations of stationary phases

Methodological considerations for molecular dynamics simulations of complex interfaces are presented in this article. A slab geometry is examined in the context of stationary phases where selectivity occurs predominantly in pores within silica beads. Specifically, we examine the Whelk-O1 interface with n-hexane/2-propanol, the TMA-(Pro)(2)-N(CH(3))-tether interface with n-hexane/2-propanol, and the C(18)H(37)Si interface with water/methanol. The following methodological issues are considered in detail: The assessment of solvent density within the confined region and excluded volume of the interface; the structural equilibration of surface-bound moieties; solvent equilibration for binary mixtures; surface size effects, and periodic boundary conditions; the treatment of electrostatic interactions; and the impact of pore size.     

Assessment of biomolecular force fields for molecular dynamics simulations in a protein crystal

Different biomolecular force fields (OPLS‐AA, AMBER03, and GROMOS96) in conjunction with SPC, SPC/E and TIP3P water models are assessed for molecular dynamics simulations in a tetragonal lysozyme crystal. The root mean square deviations for the C_a atoms of lysozymes are about 0.1 to 0.2 nm from OPLS‐AA and AMBER03, smaller than 0.4 nm from GROMOS96. All force fields exhibit similar pattern in B‐factors, whereas OPLS‐AA and AMBER03 accurately reproduce experimental measurements. Despite slight variations, the primary secondary structures are well conserved using different force fields. Water diffusion in the crystal is approximately ten‐fold slower than in bulk phase. The directional and average water diffusivities from OPLS‐AA and AMBER03 along with SPC/E model match fairly well with experimental data. Compared to GROMOS96, OPLS‐AA and AMBER03 predict larger hydrophilic solvent‐accessible surface area of lysozyme, more hydrogen bonds between lysozyme and water, and higher percentage of water in hydration shell. SPC, SPC/E and TIP3P water models have similar performance in most energetic and structural properties, but SPC/E outperforms in water diffusion. While all force fields overestimate the mobility and electrical conductivity of NaCl, a combination of OPLS‐AA for lysozyme and the Kirkwood‐Buff model for ions is superior to others. As attributed to the steric restraints and surface interactions, the mobility and conductivity in the crystal are reduced by one to two orders of magnitude from aqueous solution. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2010     

Estimation of protein folding probability from equilibrium simulations

The assumption that similar structures have similar folding probabilities (p(fold)) leads naturally to a procedure to evaluate p(fold) for every snapshot saved along an equilibrium folding-unfolding trajectory of a structured peptide or protein. The procedure utilizes a structurally homogeneous clustering and does not require any additional simulation. It can be used to detect multiple folding pathways as shown for a three-stranded antiparallel beta-sheet peptide investigated by implicit solvent molecular dynamics simulations.