Principal component and normal mode analysis of proteins; a quantitative comparison using the GroEL subunit

Principal component analysis (PCA) and normal mode analysis (NMA) have emerged as two invaluable tools for studying conformational changes in proteins. To compare these approaches for studying protein dynamics, we have used a subunit of the GroEL chaperone, whose dynamics is well characterized. We first show that both PCA on trajectories from molecular dynamics (MD) simulations and NMA reveal a general dynamical behavior in agreement with what has previously been described for GroEL. We thus compare the reproducibility of PCA on independent MD runs and subsequently investigate the influence of the length of the MD simulations. We show that there is a relatively poor one-to-one correspondence between eigenvectors obtained from two independent runs and conclude that caution should be taken when analyzing principal components individually. We also observe that increasing the simulation length does not improve the agreement with the experimental structural difference. In fact, relatively short MD simulations are sufficient for this purpose. We observe a rapid convergence of the eigenvectors (after ca. 6 ns). Although there is not always a clear one-to-one correspondence, there is a qualitatively good agreement between the movements described by the first five modes obtained with the three different approaches; PCA, all-atoms NMA, and coarse-grained NMA. It is particularly interesting to relate this to the computational cost of the three methods. The results we obtain on the GroEL subunit contribute to the generalization of robust and reproducible strategies for the study of protein dynamics, using either NMA or PCA of trajectories from MD simulations.     

Molecular dynamics simulation of proteins under high pressure: Structure,function and thermodynamics

BackgroundMolecular dynamics (MD) simulation is well-recognized as a powerful tool to investigate protein structure, function, and thermodynamics. MD simulation is also used to investigate high pressure effects on proteins. For conducting better MD simulation under high pressure, the main issues to be addressed are: (i) protein force fields and water models were originally developed to reproduce experimental properties obtained at ambient pressure; and (ii) the timescale to observe the pressure effect is often much longer than that of conventional MD simulations.Scope of reviewFirst, we describe recent developments in MD simulation methodologies for studying the high-pressure structure and dynamics of protein molecules. These developments include force fields for proteins and water molecules, and enhanced simulation techniques. Then, we summarize recent studies of MD simulations of proteins in water under high pressure.Major conclusionsRecent MD simulations of proteins in solution under pressure have reproduced various phenomena identified by experiments using high pressure, such as hydration, water penetration, conformational change, helix stabilization, and molecular stiffening.General significanceMD simulations demonstrate differences in the properties of proteins and water molecules between ambient and high-pressure conditions. Comparing the results obtained by MD calculations with those obtained experimentally could reveal the mechanism by which biological molecular machines work well in collaboration with water molecules.     

Ensemble MD simulations restrained via crystallographic data: Accurate structure leads to accurate dynamics

Currently, the best existing molecular dynamics (MD) force fields cannot accurately reproduce the global free‐energy minimum which realizes the experimental protein structure. As a result, long MD trajectories tend to drift away from the starting coordinates (e.g., crystallographic structures). To address this problem, we have devised a new simulation strategy aimed at protein crystals. An MD simulation of protein crystal is essentially an ensemble simulation involving multiple protein molecules in a crystal unit cell (or a block of unit cells). To ensure that average protein coordinates remain correct during the simulation, we introduced crystallography‐based restraints into the MD protocol. Because these restraints are aimed at the ensemble‐average structure, they have only minimal impact on conformational dynamics of the individual protein molecules. So long as the average structure remains reasonable, the proteins move in a native‐like fashion as dictated by the original force field. To validate this approach, we have used the data from solid‐state NMR spectroscopy, which is the orthogonal experimental technique uniquely sensitive to protein local dynamics. The new method has been tested on the well‐established model protein, ubiquitin. The ensemble‐restrained MD simulations produced lower crystallographic R factors than conventional simulations; they also led to more accurate predictions for crystallographic temperature factors, solid‐state chemical shifts, and backbone order parameters. The predictions for ¹⁵N relaxation rates are at least as accurate as those obtained from conventional simulations. Taken together, these results suggest that the presented trajectories may be among the most realistic protein MD simulations ever reported. In this context, the ensemble restraints based on high‐resolution crystallographic data can be viewed as protein‐specific empirical corrections to the standard force fields.     

Mass-weighted molecular dynamics simulation of cyclic polypeptides.

A modified molecular dynamics (MD) method in which atomic masses are weighted was developed previously for studying the conformational flexibility of neuroregulating tetrapeptide Phe-Met-Arg-Phe-amide (FMRF-amide). The method has now been applied to longer and constrained molecules, namely a disulfide-linked cyclic hexapeptide, c[CYFQNC], and its linear and "pseudo-cyclic" analogues. The sampling of dehedral conformational space of teh linear hexapeptide in mass-weighted MD simulations was found to be improved significantly over conventional MD simulations, as in the case of the shorter FMRF-amide molecule studied previously. In the cyclic hexapeptide, the internal constraint of the molecule due to the intramolecular disulfide bond (hence the absence of free terminals in the molecule) does not adversely affect the significant improvement of conformational sampling in mass-weighted MD simulations over normal MD simulations. The pseudo-cyclic polypeptide is identical to the linear CYFQNC molecule in amino acid sequence (i.e., side chains of the cysteine residues are reduced), but the positions of its two terminal heavy atoms were held fixed in space such that the molecule has a nearly cyclic conformation. For this molecule, the mass-weighted MD simulation generated a wide range of polypeptide backbone conformations covering the internal dihedral degrees of freedom; moreover, the physical space of the pseudo-cyclic structure was also sampled in a complete revolution of the entire molecular fragment about the two fixed termini during the simulation. These characteristics suggest that mass-weighted MD can also be an extremely useful method for conformational analyses of constrained molecules and, in particular, for modeling loops on protein surfaces.     

Exploring the energy landscape of protein folding using replica-exchange and conventional molecular dynamics simulations

Two independent replica-exchange molecular dynamics (REMD) simulations with an explicit water model were performed of the Trp-cage mini-protein. In the first REMD simulation, the replicas started from the native conformation, while in the second they started from a nonnative conformation. Initially, the first simulation yielded results qualitatively similar to those of two previously published REMD simulations: the protein appeared to be over-stabilized, with the predicted melting temperature 50-150K higher than the experimental value of 315K. However, as the first REMD simulation progressed, the protein unfolded at all temperatures. In our second REMD simulation, which starts from a nonnative conformation, there was no evidence of significant folding. Transitions from the unfolded to the folded state did not occur on the timescale of these simulations, despite the expected improvement in sampling of REMD over conventional molecular dynamics (MD) simulations. The combined 1.42 micros of simulation time was insufficient for REMD simulations with different starting structures to converge. Conventional MD simulations at a range of temperatures were also performed. In contrast to REMD, the conventional MD simulations provide an estimate of Tm in good agreement with experiment. Furthermore, the conventional MD is a fraction of the cost of REMD and continuous, realistic pathways of the unfolding process at atomic resolution are obtained.     

Refinement of protein structure homology models via long, all-atom molecular dynamics simulations

Accurate computational prediction of protein structure represents a longstanding challenge in molecular biology and structure-based drug design. Although homology modeling techniques are widely used to produce low-resolution models, refining these models to high resolution has proven difficult. With long enough simulations and sufficiently accurate force fields, molecular dynamics (MD) simulations should in principle allow such refinement, but efforts to refine homology models using MD have for the most part yielded disappointing results. It has thus far been unclear whether MD-based refinement is limited primarily by accessible simulation timescales, force field accuracy, or both. Here, we examine MD as a technique for homology model refinement using all-atom simulations, each at least 100 μs long-more than 100 times longer than previous refinement simulations-and a physics-based force field that was recently shown to successfully fold a structurally diverse set of fast-folding proteins. In MD simulations of 24 proteins chosen from the refinement category of recent Critical Assessment of Structure Prediction (CASP) experiments, we find that in most cases, simulations initiated from homology models drift away from the native structure. Comparison with simulations initiated from the native structure suggests that force field accuracy is the primary factor limiting MD-based refinement. This problem can be mitigated to some extent by restricting sampling to the neighborhood of the initial model, leading to structural improvement that, while limited, is roughly comparable to the leading alternative methods.     

Molecular dynamics simulations highlight mobile regions in proteins: A novel suggestion for converting a murine V(H) domain into a more tractable species

The V(H) region of the murine antibody 1F7 has been identified as a single-domain chorismate mutase, but a tendency to denature and aggregate has hampered its biochemical characterization. Standard mutagenesis approaches targeting antibody chain dimerization areas have been exhausted. We describe a new approach to the problem, where we use molecular dynamics (MD) simulations to find the differences between the untractable protein and the known soluble V(H) domain from a llama antibody. MD simulations of proteins yield information on the relative stability and fluctuations of parts of the proteins. By comparing simulation results of two related proteins their differences in stability and fluctuations can be analyzed and may suggest mutations aimed at (de)stabilization of one of the two proteins. For the mouse versus llama simulations, this approach highlights an untried area in the protein which shows increased fluctuations. The replacement of this eight-residue segment with the corresponding llama sequence gave a chimeric mutant which shows significantly decreased fluctuations. We see this as a general scheme to generate suggestions for mutagenesis experiments, not only obviously generalizable to other immunoglobulin domains, but to other protein systems as well.     

Hybrid approaches to molecular simulation

Molecular dynamics (MD) simulation is an established method for studying the conformational changes that are important for protein function. Recent advances in hardware and software have allowed MD simulations over the same timescales as experiment, improving the agreement between theory and experiment to a large extent. However, running such simulations are costly, in terms of resources, storage, and trajectory analysis. There is still a place for techniques that involve short MD simulations. In order to overcome the sampling paucity of short time-scales, hybrid methods that include some form of MD simulation can exploit certain features of the system of interest, often combining experimental information in surprising ways. Here, we review some recent hybrid approaches to the simulation of proteins.     

Structure, dynamics and hydration of the nogalamycin-d(ATGCAT)2Complex determined by NMR and molecular dynamics simulations in solution.

The structure of the 1:1 nogalamycin:d(ATGCAT)2 complex has been determined in solution from high-resolution NMR data and restrained molecular dynamics (rMD) simulations using an explicit solvation model. The antibiotic intercalates at the 5'-TpG step with the nogalose lying along the minor groove towards the centre of the duplex. Many drug-DNA nuclear Overhauser enhancements (NOEs) in the minor groove are indicative of hydrophobic interactions over the TGCA sequence. Steric occlusion prevents a second nogalamycin molecule from binding at the symmetry-related 5'-CpA site, leading to the conclusion that the observed binding orientation in this complex is the preferred orientation free of the complication of end-effects (drug molecules occupy terminal intercalation sites in all X-ray structures) or steric interactions between drug molecules (other NMR structures have two drug molecules bound in close proximity), as previously suggested. Fluctuations in key structural parameters such as rise, helical twist, slide, shift, buckle and sugar pucker have been examined from an analysis of the final 500 ps of a 1 ns rMD simulation, and reveal that many sequence-dependent structural features previously identified by comparison of different X-ray structures lie within the range of dynamic fluctuations observed in the MD simulations. Water density calculations on MD simulation data reveal a time-averaged pattern of hydration in both the major and minor groove, in good agreement with the extensive hydration observed in two related X-ray structures in which nogalamycin is bound at terminal 5'-TpG sites. However, the pattern of hydration determined from the sign and magnitude of NOE and ROE cross-peaks to water identified in 2D NOESY and ROESY experiments identifies only a few "bound" water molecules with long residence times. These solvate the charged bicycloaminoglucose sugar ring, suggesting an important role for water molecules in mediating drug-DNA electrostatic interactions within the major groove. The high density of water molecules found in the minor groove in X-ray structures and MD simulations is found to be associated with only weakly bound solvent in solution.     

Molecular dynamics simulations of peptides and proteins with a continuum electrostatic model based on screened Coulomb potentials

A continuum electrostatics approach for molecular dynamics (MD) simulations of macromolecules is presented and analyzed for its performance on a peptide and a globular protein. The approach incorporates the screened Coulomb potential (SCP) continuum model of electrostatics, which was reported earlier. The model was validated in a broad set of tests some of which were based on Monte Carlo simulations that included single amino acids, peptides, and proteins. The implementation for large-scale MD simulations presented in this article is based on a pairwise potential that makes the electrostatic model suitable for fast analytical calculation of forces. To assess the suitability of the approach, a preliminary validation is conducted, which consists of (i) a 3-ns MD simulation of the immunoglobulin-binding domain of streptococcal protein G, a 56-residue globular protein and (ii) a 3-ns simulation of Dynorphin, a biological peptide of 17 amino acids. In both cases, the results are compared with those obtained from MD simulations using explicit water (EW) molecules in an all-atom representation. The initial structure of Dynorphin was assumed to be an alpha-helix between residues 1 and 9 as suggested from NMR measurements in micelles. The results obtained in the MD simulations show that the helical structure collapses early in the simulation, a behavior observed in the EW simulation and consistent with spectroscopic data that suggest that the peptide may adopt mainly an extended conformation in water. The dynamics of protein G calculated with the SCP implicit solvent model (SCP-ISM) reveals a stable structure that conserves all the elements of secondary structure throughout the entire simulation time. The average structures calculated from the trajectories with the implicit and explicit solvent models had a cRMSD of 1.1 A, whereas each average structure had a cRMSD of about 0.8A with respect to the X-ray structure. The main conformational differences of the average structures with respect to the crystal structure occur in the loop involving residues 8-14. Despite the overall similarity of the simulated dynamics with EW and SCP models, fluctuations of side-chains are larger when the implicit solvent is used, especially in solvent exposed side-chains. The MD simulation of Dynorphin was extended to 40 ns to study its behavior in an aqueous environment. This long simulation showed that the peptide has a tendency to form an alpha-helical structure in water, but the stabilization free energy is too weak, resulting in frequent interconversions between random and helical conformations during the simulation time. The results reported here suggest that the SCP implicit solvent model is adequate to describe electrostatic effects in MD simulation of both peptides and proteins using the same set of parameters. It is suggested that the present approach could form the basis for the development of a reliable and general continuum approach for use in molecular biology, and directions are outlined for attaining this long-term goal.     

A structure-based simulation approach for electron paramagnetic resonance spectra using molecular and stochastic dynamics simulations

Electron paramagnetic resonance (EPR) spectroscopy using site-directed spin-labeling is an appropriate technique to analyze the structure and dynamics of flexible protein regions as well as protein-protein interactions under native conditions. The analysis of a set of protein mutants with consecutive spin-label positions leads to the identification of secondary and tertiary structure elements. In the first place, continuous-wave EPR spectra reflect the motional freedom of the spin-label specifically linked to a desired site within the protein. EPR spectra calculations based on molecular dynamics (MD) and stochastic dynamics simulations facilitate verification or refinement of predicted computer-aided models of local protein conformations. The presented spectra simulation algorithm implies a specialized in vacuo MD simulation at 600 K with additional restrictions to sample the entire accessible space of the bound spin-label without large temporal effort. It is shown that the distribution of spin-label orientations obtained from such MD simulations at 600 K agrees well with the extrapolated motion behavior during a long timescale MD at 300 K with explicit water. The following potential-dependent stochastic dynamics simulation combines the MD data about the site-specific orientation probabilities of the spin-label with a realistic rotational diffusion coefficient yielding a set of trajectories, each more than 700 ns long, essential to calculate the EPR spectrum. Analyses of a structural model of the loop between helices E and F of bacteriorhodopsin are illustrated to demonstrate the applicability and potentials of the reported simulation approach. Furthermore, effects on the motional freedom of bound spin-labels induced by solubilization of bacteriorhodopsin with Triton X-100 are examined.     

Molecular dynamics simulations and membrane protein structure quality

Despite a growing repertoire of membrane protein structures (currently ∼120 unique structures), considerations of low resolution and crystallization in the absence of a lipid bilayer require the development of techniques to assess the global quality of membrane protein folds. This is also the case for assessment of, e.g. homology models of human membrane proteins based on structures of (distant) bacterial homologues. Molecular dynamics (MD) simulations may be used to help evaluate the quality of a membrane protein structure or model. We have used a structure of the bacterial ABC transporter MsbA which has the correct transmembrane helices but an incorrect handedness and topology of their packing to test simulation methods of quality assessment. An MD simulation of the MsbA model in a lipid bilayer is compared to a simulation of another bacterial ABC transporter, BtuCD. The latter structure has demonstrated good conformational stability in the same bilayer environment and over the same timescale (20 ns) as for the MsbA model simulation. A number of comparative analyses of the two simulations were performed to assess changes in the structural integrity of each protein. The results show a significant difference between the two simulations, chiefly due to the dramatic structural deformations of MsbA. We therefore propose that MD could become a useful quality control tool for membrane protein structural biology. In particular, it provides a way in which to explore the global conformational stability of a model membrane protein fold.     

Defining Coarse-Grained Representations of Large Biomolecules and Biomolecular Complexes from Elastic Network Models

Coarse-grained (CG) models of large biomolecular complexes enable simulations of these systems over long timescales that are not accessible for atomistic molecular dynamics (MD) simulations. A systematic methodology, called essential dynamics coarse-graining (ED-CG), has been developed for defining coarse-grained sites in a large biomolecule. The method variationally determines the CG sites so that key dynamic domains in the protein are preserved in the CG representation. The original ED-CG method relies on a principal component analysis (PCA) of a MD trajectory. However, for many large proteins and multi-protein complexes such an analysis may not converge or even be possible. This work develops a new ED-CG scheme using an elastic network model (ENM) of the protein structure. In this procedure, the low-frequency normal modes obtained by ENM are used to define dynamic domains and to define the CG representation accordingly. The method is then applied to several proteins, such as the HIV-1 CA protein dimer, ATP-bound G-actin, and the Arp2/3 complex. Numerical results show that ED-CG with ENM (ENM-ED-CG) is much faster than ED-CG with PCA because no MD is necessary. The ENM-ED-CG models also capture functional essential dynamics of the proteins almost as well as those using full MD with PCA. Therefore, the ENM-ED-CG method may be better suited to coarse-grain a very large biomolecule or biomolecular complex that is too computationally expensive to be simulated by conventional MD, or when a high resolution atomic structure is not even available.     

Event detection and sub‐state discovery from biomolecular simulations using higher‐order statistics: Application to enzyme adenylate kinase

Biomolecular simulations at millisecond and longer time‐scales can provide vital insights into functional mechanisms. Because post‐simulation analyses of such large trajectory datasets can be a limiting factor in obtaining biological insights, there is an emerging need to identify key dynamical events and relating these events to the biological function online, that is, as simulations are progressing. Recently, we have introduced a novel computational technique, quasi‐anharmonic analysis (QAA) (Ramanathan et al., PLoS One 2011;6:e15827), for partitioning the conformational landscape into a hierarchy of functionally relevant sub‐states. The unique capabilities of QAA are enabled by exploiting anharmonicity in the form of fourth‐order statistics for characterizing atomic fluctuations. In this article, we extend QAA for analyzing long time‐scale simulations online. In particular, we present HOST4MD—a higher‐order statistical toolbox for molecular dynamics simulations, which (1) identifies key dynamical events as simulations are in progress, (2) explores potential sub‐states, and (3) identifies conformational transitions that enable the protein to access those sub‐states. We demonstrate HOST4MD on microsecond timescale simulations of the enzyme adenylate kinase in its apo state. HOST4MD identifies several conformational events in these simulations, revealing how the intrinsic coupling between the three subdomains (LID, CORE, and NMP) changes during the simulations. Further, it also identifies an inherent asymmetry in the opening/closing of the two binding sites. We anticipate that HOST4MD will provide a powerful and extensible framework for detecting biophysically relevant conformational coordinates from long time‐scale simulations. Proteins 2012. © 2012 Wiley Periodicals, Inc.     

Characterization of Protein Flexibility Using Small-Angle X-Ray Scattering and Amplified Collective Motion Simulations

Large-scale flexibility within a multidomain protein often plays an important role in its biological function. Despite its inherent low resolution, small-angle x-ray scattering (SAXS) is well suited to investigate protein flexibility and determine, with the help of computational modeling, what kinds of protein conformations would coexist in solution. In this article, we develop a tool that combines SAXS data with a previously developed sampling technique called amplified collective motions (ACM) to elucidate structures of highly dynamic multidomain proteins in solution. We demonstrate the use of this tool in two proteins, bacteriophage T4 lysozyme and tandem WW domains of the formin-binding protein 21. The ACM simulations can sample the conformational space of proteins much more extensively than standard molecular dynamics (MD) simulations. Therefore, conformations generated by ACM are significantly better at reproducing the SAXS data than are those from MD simulations.     

Characterization of Protein Flexibility Using Small-Angle X-Ray Scattering and Amplified Collective Motion Simulations

Large-scale flexibility within a multidomain protein often plays an important role in its biological function. Despite its inherent low resolution, small-angle x-ray scattering (SAXS) is well suited to investigate protein flexibility and determine, with the help of computational modeling, what kinds of protein conformations would coexist in solution. In this article, we develop a tool that combines SAXS data with a previously developed sampling technique called amplified collective motions (ACM) to elucidate structures of highly dynamic multidomain proteins in solution. We demonstrate the use of this tool in two proteins, bacteriophage T4 lysozyme and tandem WW domains of the formin-binding protein 21. The ACM simulations can sample the conformational space of proteins much more extensively than standard molecular dynamics (MD) simulations. Therefore, conformations generated by ACM are significantly better at reproducing the SAXS data than are those from MD simulations.     

Low-mass molecular dynamics simulation: A simple and generic technique to enhance configurational sampling

CLN025 is one of the smallest fast-folding proteins. Until now it has not been reported that CLN025 can autonomously fold to its native conformation in a classical, all-atom, and isothermal–isobaric molecular dynamics (MD) simulation. This article reports the autonomous and repeated folding of CLN025 from a fully extended backbone conformation to its native conformation in explicit solvent in multiple 500-ns MD simulations at 277 K and 1 atm with the first folding event occurring as early as 66.1 ns. These simulations were accomplished by using AMBER forcefield derivatives with atomic masses reduced by 10-fold on Apple Mac Pros. By contrast, no folding event was observed when the simulations were repeated using the original AMBER forcefields of FF12SB and FF14SB. The results demonstrate that low-mass MD simulation is a simple and generic technique to enhance configurational sampling. This technique may propel autonomous folding of a wide range of miniature proteins in classical, all-atom, and isothermal–isobaric MD simulations performed on commodity computers—an important step forward in quantitative biology.     

Unfolding stabilities of two structurally similar proteins as probed by temperature-induced and force-induced molecular dynamics simulations

Unfolding stabilities of two homologous proteins, cardiotoxin III and short-neurotoxin (SNTX) belonging to three-finger toxin (TFT) superfamily, have been probed by means of molecular dynamics (MD) simulations. Combined analysis of data obtained from steered MD and all-atom MD simulations at various temperatures in near physiological conditions on the proteins suggested that overall structural stabilities of the two proteins were different from each other and the MD results are consistent with experimental data of the proteins reported in the literature. Rationalization for the differential structural stabilities of the structurally similar proteins has been chiefly attributed to the differences in the structural contacts between C- and N-termini regions in their three-dimensional structures, and the findings endorse the ‘CN network’ hypothesis proposed to qualitatively analyse the thermodynamic stabilities of proteins belonging to TFT superfamily of snake venoms. Moreover, the ‘CN network’ hypothesis has been revisited and the present study suggested that ‘CN network’ should be accounted in terms of ‘structural contacts’ and ‘structural strengths’ in order to precisely describe order of structural stabilities of TFTs.     

2H-NMR study and molecular dynamics simulation of the location, alignment, and mobility of pyrene in POPC bilayers

The alignment of pyrene in a 1-palmitoyl-2-oleoyl-phosphatidylcholine bilayer was investigated using two different approaches, namely solid-state (2)H-NMR spectroscopy and molecular dynamics (MD) simulations. Quadrupolar splittings from (2)H-NMR spectra of deuterated pyrene-d(10) in an oriented lipid bilayer give information about the orientation of C-D bonds with respect to the membrane normal. From MD simulations, geometric information is accessible via trajectories. By defining molecular and bond order parameters, the data from MD trajectories and NMR spectra can be compared straightforwardly. To ensure that the results from both methods are comparable, parameters of the experimental and the simulation setup were chosen to be as similar as possible. From simulations, we saw that pyrene prefers a position inside the lipid membrane near the headgroups and has no tendency to diffuse from one monolayer of the membrane to the other. The results from simulation and NMR show that the normal of the molecular plane is aligned nearly perpendicular to the bilayer normal. The long axis of pyrene lies preferentially parallel to the bilayer normal within a range of +/-30 degrees . The results from the two different methods are remarkably consistent. The good agreement can be explained by the fact that the different kind of motions of a pyrene molecule are already averaged within a few nanoseconds, which is the timescale covered by the MD simulation.