A Comparison Between Elastic Network Interpolation and MD Simulation of 16S Ribosomal RNA

Abstract

In this paper a coarse-grained method called elastic network interpolation (ENI) is used to generate feasible transition pathways between two given conformations of the core central domain of 16S Ribosomal RNA (16S rRNA). The two given conformations are the extremes generated by a molecular dynamics (MD) simulation, which differ from each other by 10Å in root-mean-square deviation (RMSD). It takes only several hours to build an ENI pathway on a 1.5GHz Pentium with 512 MB memory, while the MD takes several weeks on high-performance multi-processor servers such as the SGI ORIGIN 2000/2100. It is shown that multiple ENI pathways capture the essential anharmonic motions of millions of timesteps in a particular MD simulation. A coarse-grained normal mode analysis (NMA) is performed on each intermediate ENI conformation, and the lowest 1% of the normal modes (representing about 40 degrees of freedom (DOF)) are used to parameterize fluctuations. This combined ENI/NMA method captures all intermediate conformations in the MD run with 1.5Å RMSD on average. In addition, if we restrict attention to the time interval of the MD run between the two extreme conformations, the RMSD between the closest ENI/NMA pathway and the MD results is about 1Å. These results may serve as a paradigm for reducedDOF dynamic simulations of large biological macromolecules as well as a method for the reduced-parameter interpretation of massive amounts of MD data.     

Rigid-cluster models of conformational transitions in macromolecular machines and assemblies

We present a rigid-body-based technique (called rigid-cluster elastic network interpolation) to generate feasible transition pathways between two distinct conformations of a macromolecular assembly. Many biological molecules and assemblies consist of domains which act more or less as rigid bodies during large conformational changes. These collective motions are thought to be strongly related with the functions of a system. This fact encourages us to simply model a macromolecule or assembly as a set of rigid bodies which are interconnected with distance constraints. In previous articles, we developed coarse-grained elastic network interpolation (ENI) in which, for example, only Calpha atoms are selected as representatives in each residue of a protein. We interpolate distance differences of two conformations in ENI by using a simple quadratic cost function, and the feasible conformations are generated without steric conflicts. Rigid-cluster interpolation is an extension of the ENI method with rigid-clusters replacing point masses. Now the intermediate conformations in an anharmonic pathway can be determined by the translational and rotational displacements of large clusters in such a way that distance constraints are observed. We present the derivation of the rigid-cluster model and apply it to a variety of macromolecular assemblies. Rigid-cluster ENI is then modified for a hybrid model represented by a mixture of rigid clusters and point masses. Simulation results show that both rigid-cluster and hybrid ENI methods generate sterically feasible pathways of large systems in a very short time. For example, the HK97 virus capsid is an icosahedral symmetric assembly composed of 60 identical asymmetric units. Its original Hessian matrix size for a Calpha coarse-grained model is >(300,000)(2). However, it reduces to (84)(2) when we apply the rigid-cluster model with icosahedral symmetry constraints. The computational cost of the interpolation no longer scales heavily with the size of structures; instead, it depends strongly on the minimal number of rigid clusters into which the system can be decomposed.     

NMMD: Efficient Cryo-EM Flexible Fitting Based on Simultaneous Normal Mode and Molecular Dynamics atomic displacements

Atomic models of cryo electron microscopy (cryo-EM) maps of biomolecular conformations are often obtained by flexible fitting of the maps with available atomic structures of other conformations (e.g., obtained by X-ray crystallography). This article presents a new flexible fitting method, NMMD, which combines normal mode analysis (NMA) and molecular dynamics simulation (MD). Given an atomic structure and a cryo-EM map to fit, NMMD simultaneously estimates global atomic displacements based on NMA and local displacements based on MD. NMMD was implemented by modifying EMfit, a flexible fitting method using MD only, in GENESIS 1.4. As EMfit, NMMD can be run with replica exchange umbrella sampling procedure. The new method was tested using a variety of EM maps (synthetic and experimental, with different noise levels and resolutions). The results of the tests show that adding normal modes to MD-based fitting makes the fitting faster (40% in average) and, in the majority of cases, more accurate.     

Study of the Villin headpiece folding dynamics by combining coarse-grained Monte Carlo evolution and all-atom molecular dynamics

The folding mechanism of the Villin headpiece (HP36) is studied by means of a novel approach which entails an initial coarse-grained Monte Carlo (MC) scheme followed by all-atom molecular dynamics (MD) simulations in explicit solvent. The MC evolution occurs in a simplified free-energy landscape and allows an efficient selection of marginally-compact structures which are taken as viable initial conformations for the MD. The coarse-grained MC structural representation is connected to the one with atomic resolution through a "fine-graining" reconstruction algorithm. This two-stage strategy is used to select and follow the dynamics of seven different unrelated conformations of HP36. In a notable case the MD trajectory rapidly evolves towards the folded state, yielding a typical root-mean-square deviation (RMSD) of the core region of only 2.4 A from the closest NMR model (the typical RMSD over the whole structure being 4.0 A). The analysis of the various MC-MD trajectories provides valuable insight into the details of the folding and mis-folding mechanisms and particularly about the delicate influence of local and nonlocal interactions in steering the folding process.     

Conformational pathway for the kissing complex-->extended dimer transition of the SL1 stem-loop from genomic HIV-1 RNA as monitored by targeted molecular dynamics techniques

HIV-1 retroviral genomic RNA dimerization is initiated by loop-loop interactions between the SL1 stem-loops of two identical RNA molecules. The SL1-SL1 unstable resulting kissing complex (KC) then refolds irreversibly into a more stable complex called extended dimer (ED). Although the structures of both types of complex have been determined, very little is known about the conformational pathway corresponding to the transition, owing to the difficulty of observing experimentally intermediate conformations. In this study, we applied targeted molecular dynamics simulation techniques (TMD) to the phosphorus atoms for monitoring this pathway for the backbone, and a two-step strategy was adopted. In a first step, called TMD(-1), the dimer structure was constrained to progressively move away from KC without indicating the direction, until the RMSD from KC reaches 36A. A total of 20 TMD(-1) simulations were performed under different initial conditions and different simulation parameters. For RMSD ranging between 0 and 22A, the whole set of TMD(-1) simulations follows a similar pathway, then divergences are observed. None of the simulations leads to the ED structure. At RMSD=22A, the dimers look like two parallel Us, still linked by the initial loop-loop interaction, but the strands of the stems (the arms of the Us) are positioned in such a manner that they can form intramolecular as well as intermolecular Watson-Crick base-pairs. This family of structure is called UU. In a second step (TMD simulations), 18 structures were picked up along the pathways generated with TMD(-1) and were constrained to move toward ED by decreasing progressively their RMSD from ED. We found that only structures from the UU family are able to easily reach ED-like conformations of the backbones without exhibiting a large constraint energy.     

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.     

A one-dimensional reaction coordinate for identification of transition states from explicit solvent P(fold)-like calculations

A properly identified transition state ensemble (TSE) in a molecular dynamics (MD) simulation can reveal a tremendous amount about how a protein folds and offer a point of comparison to experimentally derived Phi(F) values, which reflect the degree of structure in these transient states. In one such method of TSE identification, dubbed P(fold), MD simulations of individual protein structures taken from an unfolding trajectory are used to directly assess an input structure's probability of folding before unfolding, and P(fold) is, by definition, 0.5 for the TSE. Other, less computationally intensive methods, such as multidimensional scaling (MDS) of the pairwise root mean-squared deviation (RMSD) matrix of the conformations sampled in a thermal unfolding trajectory, have also been used to identify the TSE. Identification of the TSE is made from the original MD simulation without the need to run further simulations. Here we present a P(fold)-like study and describe methods for identification of the TSE through the derivation of a high fidelity, bounded, one-dimensional reaction coordinate for protein folding. These methods are applied to the engrailed homeodomain. The TSE identified by this approach is essentially identical to the TSE identified previously by MDS of the pairwise RMSD matrix. However, the cost of performing P(fold), or even our reduced P(fold)-like calculations, is at least 36,000 times greater than the MDS method.     

Opening and closing of a toroidal group II chaperonin revealed by a symmetry constrained elastic network model

Recently, the atomic structures of both the closed and open forms of Group 2 chaperonin protein Mm‐cpn were revealed through crystallography and cryo‐electron microscopy. This toroidal‐like chaperonin is composed of two eightfold rings that face back‐to‐back. To gain a computational advantage, we used a symmetry constrained elastic network model (SCENM), which requires only a repeated subunit structure and its symmetric connectivity to neighboring subunits to simulate the entire system. In the case of chaperonin, only six subunits (i.e., three from each ring) were used out of the eight subunits comprising each ring. A smooth and symmetric pathway between the open and closed conformations was generated by elastic network interpolation (ENI). To support this result, we also performed a symmetry‐constrained normal mode analysis (NMA), which revealed the intrinsic vibration features of the given structures. The NMA and ENI results for the representative single subunit were duplicated according to the symmetry pattern to reconstruct the entire assembly. To test the feasibility of the symmetry model, its results were also compared with those obtained from the full model. This study allowed the folding mechanism of chaperonin Mm‐cpn to be elucidated by SCENM in a timely manner.     

Kinetic studies of folding of the B-domain of staphylococcal protein A with molecular dynamics and a united-residue (UNRES) model of polypeptide chains

Langevin dynamics is used with our physics-based united-residue (UNRES) force field to study the folding pathways of the B-domain of staphylococcal protein A (1BDD (alpha; 46 residues)). With 400 trajectories of protein A started from the extended state (to gather meaningful statistics), and simulated for more than 35 ns each, 380 of them folded to the native structure. The simulations were carried out at the optimal folding temperature of protein A with this force field. To the best of our knowledge, this is the first simulation study of protein-folding kinetics with a physics-based force field in which reliable statistics can be gathered. In all the simulations, the C-terminal alpha-helix forms first. The ensemble of the native basin has an average RMSD value of 4 A from the native structure. There is a stable intermediate along the folding pathway, in which the N-terminal alpha-helix is unfolded; this intermediate appears on the way to the native structure in less than one-fourth of the folding pathways, while the remaining ones proceed directly to the native state. Non-native structures persist until the end of the simulations, but the native-like structures dominate. To express the kinetics of protein A folding quantitatively, two observables were used: (i) the average alpha-helix content (averaged over all trajectories within a given time window); and (ii) the fraction of conformations (averaged over all trajectories within a given time window) with Calpha RMSD values from the native structure less than 5 A (fraction of completely folded structures). The alpha-helix content grows quickly with time, and its variation fits well to a single-exponential term, suggesting fast two-state kinetics. On the other hand, the fraction of folded structures changes more slowly with time and fits to a sum of two exponentials, in agreement with the appearance of the intermediate, found when analyzing the folding pathways. This observation demonstrates that different qualitative and quantitative conclusions about folding kinetics can be drawn depending on which observable is monitored.     

Flexibility of the exportins Cse1p and Xpot depicted by elastic network model

Nucleocytoplasmic transport in eukaryotic cells involves many interactions between macromolecules, and has been an active area for many researchers. However, the precise mechanism still evades us and more efforts are needed to better understand it. In this study, the authors investigated exportins (Cse1p and Xpot) by elastic network interpolation (ENI) and elastic network based normal mode analysis (EN-NMA). Results of the study on Cse1p were in good agreement with the results obtained by molecular dynamics simulation in another study but with the benefit of time-efficiency. First, a formation of ring closure obtained by ENI was observed. Second, HEAT 1 to 3 and HEAT 14 to 17 had the largest values of root mean square deviation (RMSD) which indicated the flexibility of Cse1p during the transition. In the case of Xpot, a possible pathway from nuclear state to cytoplasmic state was shown, and the predicted pathway was also quantitatively analyzed in terms of RMSD. The results suggested two flexible regions of Xpot that might be important to the transporting mechanism. Moreover, the dominant mode of Xpot in the nuclear state obtained by EN-NMA not only showed the tendency to match the predicted pathway to the cytoplasmic state of Xpot, but also displayed the flexible regions of Xpot. A time-efficient computational approach was presented in this paper and the results indicated that the flexibility of tested exportins might be required to perform the biological function of transporting cargos.     

Coarse-grained modeling of the actin filament derived from atomistic-scale simulations

A coarse-grained (CG) procedure that incorporates the information obtained from all-atom molecular dynamics (MD) simulations is presented and applied to actin filaments (F-actin). This procedure matches the averaged values and fluctuations of the effective internal coordinates that are used to define a CG model to the values extracted from atomistic MD simulations. The fluctuations of effective internal coordinates in a CG model are computed via normal-mode analysis (NMA), and the computed fluctuations are matched with the atomistic MD results in a self-consistent manner. Each actin monomer (G-actin) is coarse-grained into four sites, and each site corresponds to one of the subdomains of G-actin. The potential energy of a CG G-actin contains three bonds, two angles, and one dihedral angle; effective harmonic bonds are used to describe the intermonomer interactions in a CG F-actin. The persistence length of a CG F-actin was found to be sensitive to the cut-off distance of assigning intermonomer bonds. Effective harmonic bonds for a monomer with its third nearest neighboring monomers are found to be necessary to reproduce the values of persistence length obtained from all-atom MD simulations. Compared to the elastic network model, incorporating the information of internal coordinate fluctuations enhances the accuracy and robustness for a CG model to describe the shapes of low-frequency vibrational modes. Combining the fluctuation-matching CG procedure and NMA, the achievable time- and length scales of modeling actin filaments can be greatly enhanced. In particular, a method is described to compute the force-extension curve using the CG model developed in this work and NMA. It was found that F-actin is easily buckled under compressive deformation, and a writhing mode is developed as a result. In addition to the bending and twisting modes, this novel writhing mode of F-actin could also play important roles in the interactions of F-actin with actin-binding proteins and in the force-generation process via polymerization.     

Mechanism of fiber assembly: treatment of Aβ peptide aggregation with a coarse-grained united-residue force field

The growth mechanism of β-amyloid (Aβ) peptide fibrils was studied by a physics-based coarse-grained united-residue model and molecular dynamics (MD) simulations. To identify the mechanism of monomer addition to an Aβ_1-40 fibril, we placed an unstructured monomer at a distance of 20 Å from a fibril template and allowed it to interact freely with the latter. The monomer was not biased towards fibril conformation by either the force field or the MD algorithm. With the use of a coarse-grained model with replica-exchange molecular dynamics, a longer timescale was accessible, making it possible to observe how the monomers probe different binding modes during their search for the fibril conformation. Although different assembly pathways were seen, they all follow a dock-lock mechanism with two distinct locking stages, consistent with experimental data on fibril elongation. Whereas these experiments have not been able to characterize the conformations populating the different stages, we have been able to describe these different stages explicitly by following free monomers as they dock onto a fibril template and to adopt the fibril conformation (i.e., we describe fibril elongation step by step at the molecular level). During the first stage of the assembly (“docking”), the monomer tries different conformations. After docking, the monomer is locked into the fibril through two different locking stages. In the first stage, the monomer forms hydrogen bonds with the fibril template along one of the strands in a two-stranded β-hairpin; in the second stage, hydrogen bonds are formed along the second strand, locking the monomer into the fibril structure. The data reveal a free-energy barrier separating the two locking stages. The importance of hydrophobic interactions and hydrogen bonds in the stability of the Aβ fibril structure was examined by carrying out additional canonical MD simulations of oligomers with different numbers of chains (4-16 chains), with the fibril structure as the initial conformation. The data confirm that the structures are stabilized largely by hydrophobic interactions and show that intermolecular hydrogen bonds are highly stable and contribute to the stability of the oligomers as well.     

Efficient RMSD measures for the comparison of two molecular ensembles. Root-mean-square deviation

Quantitative measures are presented for comparing the conformations of two molecular ensembles. The measures are based on Kabsch's formula for the root-mean-square deviation (RMSD) and the covariance matrix of atomic positions of isotropically distributed ensembles (IDE). By using a Taylor series expansion, it is shown that the RMSD can be expressed solely in terms of the IDE matrices. A fast approximate method is introduced for the pairwise RMSD determination whose computational cost scales linearly with the number of structures. A similarity measure for two structural ensembles that is based on the trace metric of the differences of powers of the IDE matrices is presented. The measures are illustrated for conformational ensembles generated by a molecular dynamics computer simulation of a partially folded A-state analog of ubiquitin.     

Effects of different force fields on the structural character of α synuclein β-hairpin peptide (35–56) in aqueous environment

The hallmark of Parkinson’s disease (PD) is the intracellular protein aggregation forming Lewy Bodies (LB) and Lewy neuritis which comprise mostly of a protein, alpha synuclein (α-syn). Molecular dynamics (MD) simulation methods can augment experimental techniques to understand misfolding and aggregation pathways with atomistic resolution. The quality of MD simulations for proteins and peptides depends greatly on the accuracy of empirical force fields. The aim of this work is to investigate the effects of different force fields on the structural character of β hairpin fragment of α-syn (residues 35–56) peptide in aqueous solution. Six independent MD simulations are done in explicit solvent using, AMBER03, AMBER99SB, GROMOS96 43A1, GROMOS96 53A6, OPLS-AA, and CHARMM27 force fields with CMAP corrections. The performance of each force field is assessed from several structural parameters such as root mean square deviation (RMSD), root mean square fluctuation (RMSF), radius of gyration (Rg), solvent accessible surface area (SASA), formation of β-turn, the stability of folded β-hairpin structure, and the favourable conformations obtained for different force fields. In this study, CMAP correction of CHARMM27 force field is found to overestimate the helical conformation, while GROMOS96 53A6 is found to most successfully capture the conformational dynamics of α-syn β-hairpin fragment as elicited from NMR.     

Molecular dynamics simulations of hydrophobic collapse of ubiquitin.

Nine nonnative conformations of ubiquitin, generated during two different thermal denaturation trajectories, were simulated under nearly native conditions (62 degrees C). The simulations included all protein and solvent atoms explicitly, and simulation times ranged from 1-2.4 ns. The starting structures had alpha-carbon root-mean-square deviations (RMSDs) from the crystal structure of 4-12 A and radii of gyration as high as 1.3 times that of the native state. In all but one case, the protein collapsed when the temperature was lowered and sampled conformations as compact as those reached in a control simulation beginning from the crystal structure. In contrast, the protein did not collapse when simulated in a 60% methanol:water mixture. The behavior of the protein depended on the starting structure: during simulation of the most native-like starting structures (<5 A RMSD to the crystal structure) the RMSD decreased, the number of native hydrogen bonds increased, and the secondary and tertiary structure increased. Intermediate starting structures (5-10 A RMSD) collapsed to the radius of gyration of the control simulation, hydrophobic residues were preferentially buried, and the protein acquired some native contacts. However, the protein did not refold. The least native starting structures (10-12 A RMSD) did not collapse as completely as the more native-like structures; instead, they experienced large fluctuations in radius of gyration and went through cycles of expansion and collapse, with improved burial of hydrophobic residues in successive collapsed states.     

Fast Side Chain Replacement in Proteins Using a Coarse-Grained Approach for Evaluating the Effects of Mutation During Evolution

For high-throughput structural genomic and evolutionary bioinformatics approaches, there is a clear need for fast methods to evaluate substitutions structurally. Coarse-grained methods are both powerful and fast, and a coarse-grained approach to position the substituted side chains is presented. Through the application of a coarse-grained method, a speed-up on the single- residue replacement, of at least sevenfold is achieved compared with modern all-atom approaches. At the same time, this approach maintains a small median RMSD from the leading all-atom approach (as measured in coarse-grained space), and predicts the conformation of point mutants with similar accuracy and generates biologically realistic side chain angles. This method is also substantially more predictable in its run time, making it useful for high-throughput studies of protein structural evolution. To demonstrate the utility of this method, it has been implemented in a forward simulation of sequences threaded through the SH2 domains, with selective pressures to fold and bind specifically. The relative substitution rates across the protein structure and at the binding interface are reflective of those observed in SH2 domain evolution. The algorithm has been implemented in C++, with the source code and binaries (currently supported for Linux systems) freely available as SARA at http://www.wyomingbioinformatics.org/LiberlesGroup/SARA .     

A new Hybrid Monte Carlo algorithm for protein potential function test and structure refinement

A new Hybrid Monte Carlo (HMC) algorithm has been developed to test protein potential functions and, ultimately, refine protein structures. The main principle of this algorithm is, in each cycle, a new trial conformation is generated by carrying out a short period of molecular dynamics (MD) iterations with a set of random parameters (including the MD time step, the number of MD steps, the MD temperature, and the seed for initial MD velocity assignment); then to accept or reject the new conformation on the basis of the Metropolis criterion. The novelty in this paper is that the potential in MD iterations is different from that in the MC step. In the former, it is a molecular mechanics potential, in the latter it is a knowledge-based potential (KBP). Directed by the KBP, the MD iteration is used to search conformational space for realistic conformations with low KBP energy. It circumvents the difficulty in using KBP functions directly in MD simulation, as KBP functions are typically incomplete, and do not always have continuous derivatives required for the calculation of the forces. The new algorithm has been tested in explorations of conformational space. In these test calculations the KBP energy was found to drop below the value for the native conformation, and the correlation between the root mean square deviation (RMSD) and the KBP energy was shown to be different from the test results in other references. At the present time, the algorithm is useful for testing new KBP functions. Furthermore, if a KBP function can be found for which the native conformation has the lowest energy and the energy/RMSD correlation is good, then this new algorithm also will be a tool for refinement of the theory-based structural models.