An all-atom structure-based potential for proteins: bridging minimal models with all-atom empirical forcefields

Simulating Protein Fold Switching by Replica Exchange with Tunneling

Recent experiments suggest that an amino acid sequence encodes not only the native fold of a protein but also other forms that are essential for its function or are important during folding or association. These various forms populate a multifunnel folding and association landscape where mutations, changes in environment, or interaction with other molecules switch between the encoded folds. We introduce replica exchange with tunneling as a way to efficiently simulate switching between distinct folds of proteins and protein aggregates. The correctness and efficiency of our approach are demonstrated in a series of simulations covering a wide range of proteins, from a small 11-residue large designed peptide to two 56-residue large mutants of the A and B domains of protein G.

Capturing Transition States for tRNA Hybrid-State Formation in the Ribosome

In order to quantitatively describe the energetics of biomolecular rearrangements, it is necessary to identify reaction coordinates that accurately capture the relevant transition events. Here, we perform simulations of A-site tRNA movement (∼20 Å) during hybrid-state formation in the ribosome and quantify the ability of interatomic distances to capture the transition state ensemble. Numerous coordinates are found to be accurate indicators of the transition state, allowing tRNA rearrangements to be described as diffusion across a one-dimensional free-energy surface. In addition to providing insights into the physical-chemical relationship between biomolecular structure and dynamics, these results can help enable single-molecule techniques to probe the free-energy landscape of the ribosome.

Protein Domain-Swapping Can Be a Consequence of Functional Residues

Monomer topology has been implicated in domain-swapping, a potential first step on the route to disease-causing protein aggregation. Despite having the same topology (β1-α1-β2-β3-β4-β5), the cysteine protease inhibitor stefin-B domain swaps more readily than a single-chain variant of the heterodimeric sweet protein monellin (scMn). Here, we computationally study the folding of stefin-B and scMn in order to understand the molecular basis for the difference in their domain-swapping propensities. In agreement with experiments, our structure-based simulations show that scMn folds cooperatively without the population of an intermediate while stefin-B populates an equilibrium intermediate state. Since the simulation intermediate has only one domain structured (β3-β4-β5), it can directly lead to domain-swapping. Using computational variants of stefin-B, we show that the population of this intermediate is caused by regions of stefin-B that have been implicated in protease inhibition. We also find that the protease-binding regions are located on two structural elements and localized in space. In contrast, the residues that contribute to the sweetness of monellin are not localized to a few structural elements but are distributed over the protein fold. We conclude that the distributed functional residues of monellin do not induce large local perturbations in the protein structure, eliminating the formation of folding intermediates and in turn domain-swapping. On the other hand, the localized protease-binding regions of stefin-B promote the formation of a folding intermediate which can lead to domain-swapping. Thus, domain-swapping can be a direct consequence of the constraints that function imposes on the protein structure.

NTL9 Folding at Constant pH: The Importance of Electrostatic Interaction and pH Dependence

The folding process of the N-terminal domain of ribosomal protein L9 (NTL9) was investigated at constant-pH computer simulations. Evaluation of the role of electrostatic interaction during folding was carried out by including a Debye-Hückel potential into a Cα structure-based model (SBM). In this study, the charges of the ionizable residues and the electrostatic potential are susceptible to the solution conditions, such as pH and ionic strength, as well as to the presence of charged groups. Simulations were performed under different pHs, and the results were validated by comparing them with experimental values of pKa and with denaturation experiment data. Also, the free energy profiles, Φ-values, and folding routes were calculated for each condition. It was shown how charges vary along the folding under different pH, which is subject to different scenarios. This study reveals how simplified models can capture essential physical features, reproducing experimental results, and presenting the role of electrostatic interactions before, during, and after the transition state.

Theoretical Insights into the Biophysics of Protein Bi-stability and Evolutionary Switches

Deciphering the effects of nonsynonymous mutations on protein structure is central to many areas of biomedical research and is of fundamental importance to the study of molecular evolution. Much of the investigation of protein evolution has focused on mutations that leave a protein’s folded structure essentially unchanged. However, to evolve novel folds of proteins, mutations that lead to large conformational modifications have to be involved. Unraveling the basic biophysics of such mutations is a challenge to theory, especially when only one or two amino acid substitutions cause a large-scale conformational switch. Among the few such mutational switches identified experimentally, the one between the G_A all-α and G_B α+β folds is extensively characterized; but all-atom simulations using fully transferrable potentials have not been able to account for this striking switching behavior. Here we introduce an explicit-chain model that combines structure-based native biases for multiple alternative structures with a general physical atomic force field, and apply this construct to twelve mutants spanning the sequence variation between G_A and G_B. In agreement with experiment, we observe conformational switching from G_A to G_B upon a single L45Y substitution in the GA98 mutant. In line with the latent evolutionary potential concept, our model shows a gradual sequence-dependent change in fold preference in the mutants before this switch. Our analysis also indicates that a sharp G_A/G_B switch may arise from the orientation dependence of aromatic π-interactions. These findings provide physical insights toward rationalizing, predicting and designing evolutionary conformational switches.

The biological functions of globular proteins are intimately related to their folded structures and their associated conformational fluctuations. Evolution of new structures is an important avenue to new functions. Although many mutations do not change the folded state, experiments indicate that a single amino acid substitution can lead to a drastic change in the folded structure. The physics of this switch-like behavior remains to be elucidated. Here we develop a computational model for the relevant physical forces, showing that mutations can lead to new folds by passing through intermediate sequences where the old and new folds occur with varying probabilities. Our approach helps provide a general physical account of conformational switching in evolution and mutational effects on conformational dynamics.

Capturing transition paths and transition states for conformational rearrangements in the ribosome

To reveal the molecular determinants of biological function, one seeks to characterize the interactions that are formed in conformational and chemical transition states. In other words, what interactions govern the molecule's energy landscape? To accomplish this, it is necessary to determine which degrees of freedom can unambiguously identify each transition state. Here, we perform simulations of large-scale aminoacyl-transfer RNA (aa-tRNA) rearrangements during accommodation on the ribosome and project the dynamics along experimentally accessible atomic distances. From this analysis, we obtain evidence for which coordinates capture the correct number of barrier-crossing events and accurately indicate when the aa-tRNA is on a transition path. Although a commonly used coordinate in single-molecule experiments performs poorly, this study implicates alternative coordinates along which rearrangements are accurately described as diffusive movements across a one-dimensional free-energy profile. From this, we provide the theoretical foundation required for single-molecule techniques to uncover the energy landscape governing aa-tRNA selection by the ribosome.

Connecting thermal and mechanical protein (un)folding landscapes

Molecular dynamics simulations supplement single-molecule pulling experiments by providing the possibility of examining the full free energy landscape using many coordinates. Here, we use an all-atom structure-based model to study the force and temperature dependence of the unfolding of the protein filamin by applying force at both termini. The unfolding time-force relation τ(F) indicates that the force-induced unfolding behavior of filamin can be characterized into three regimes: barrier-limited low- and intermediate-force regimes, and a barrierless high-force regime. Slope changes of τ(F) separate the three regimes. We show that the behavior of τ(F) can be understood from a two-dimensional free energy landscape projected onto the extension X and the fraction of native contacts Q. In the low-force regime, the unfolding rate is roughly force-independent due to the small (even negative) separation in X between the native ensemble and transition state ensemble (TSE). In the intermediate-force regime, force sufficiently separates the TSE from the native ensemble such that τ(F) roughly follows an exponential relation. This regime is typically explored by pulling experiments. While X may fail to resolve the TSE due to overlap with the unfolded ensemble just below the folding temperature, the overlap is minimal at lower temperatures where experiments are likely to be conducted. The TSE becomes increasingly structured with force, whereas the average order of structural events during unfolding remains roughly unchanged. The high-force regime is characterized by barrierless unfolding, and the unfolding time approaches a limit of ∼10 μs for the highest forces we studied. Finally, a combination of X and Q is shown to be a good reaction coordinate for almost the entire force range.

Ensemble Structure of the Highly Flexible Complex Formed between Vesicular Stomatitis Virus Unassembled Nucleoprotein and its Phosphoprotein Chaperone

Nucleocapsid assembly is an essential process in the replication of the non-segmented, negative-sense RNA viruses (NNVs). Unassembled nucleoprotein (N(0)) is maintained in an RNA-free and monomeric form by its viral chaperone, the phosphoprotein (P), forming the N(0)-P complex. Our earlier work solved the structure of vesicular stomatitis virus complex formed between an N-terminally truncated N (NΔ21) and a peptide of P (P60) encompassing the N(0)-binding site, but how the full-length P interacts with N(0) remained unknown. Here, we combine several experimental biophysical methods including size exclusion chromatography with detection by light scattering and refractometry, small-angle X-ray and neutron scattering and nuclear magnetic resonance spectroscopy with molecular dynamics simulation and computational modeling to characterize the NΔ21(0)-PFL complex formed with dimeric full-length P. We show that for multi-molecular complexes, simultaneous multiple-curve fitting using small-angle neutron scattering data collected at varying contrast levels provides additional information and can help refine structural ensembles. We demonstrate that (a) vesicular stomatitis virus PFL conserves its high flexibility within the NΔ21(0)-PFL complex and interacts with NΔ21(0) only through its N-terminal extremity; (b) each protomer of P can chaperone one N(0) client protein, leading to the formation of complexes with stoichiometries 1N:P2 and 2N:P2; and (c) phosphorylation of residues Ser60, Thr62 and Ser64 provides no additional interactions with N(0) but creates a metal binding site in PNTR. A comparison with the structures of Nipah virus and Ebola virus N(0)-P core complex suggests a mechanism for the control of nucleocapsid assembly that is common to all NNVs.

SMOG 2: A Versatile Software Package for Generating Structure-Based Models

Molecular dynamics simulations with coarse-grained or simplified Hamiltonians have proven to be an effective means of capturing the functionally important long-time and large-length scale motions of proteins and RNAs. Originally developed in the context of protein folding, structure-based models (SBMs) have since been extended to probe a diverse range of biomolecular processes, spanning from protein and RNA folding to functional transitions in molecular machines. The hallmark feature of a structure-based model is that part, or all, of the potential energy function is defined by a known structure. Within this general class of models, there exist many possible variations in resolution and energetic composition. SMOG 2 is a downloadable software package that reads user-designated structural information and user-defined energy definitions, in order to produce the files necessary to use SBMs with high performance molecular dynamics packages: GROMACS and NAMD. SMOG 2 is bundled with XML-formatted template files that define commonly used SBMs, and it can process template files that are altered according to the needs of each user. This computational infrastructure also allows for experimental or bioinformatics-derived restraints or novel structural features to be included, e.g. novel ligands, prosthetic groups and post-translational/transcriptional modifications. The code and user guide can be downloaded at http://smog-server.org/smog2.

Steric interactions lead to collective tilting motion in the ribosome during mRNA-tRNA translocation

Translocation of mRNA and tRNA through the ribosome is associated with large-scale rearrangements of the head domain in the 30S ribosomal subunit. To elucidate the relationship between 30S head dynamics and mRNA–tRNA displacement, we apply molecular dynamics simulations using an all-atom structure-based model. Here we provide a statistical analysis of 250 spontaneous transitions between the A/P–P/E and P/P–E/E ensembles. Consistent with structural studies, the ribosome samples a chimeric ap/P–pe/E intermediate, where the 30S head is rotated ∼18°. It then transiently populates a previously unreported intermediate ensemble, which is characterized by a ∼10° tilt of the head. To identify the origins of head tilting, we analyse 781 additional simulations in which specific steric features are perturbed. These calculations show that head tilting may be attributed to specific steric interactions between tRNA and the 30S subunit (PE loop and protein S13). Taken together, this study demonstrates how molecular structure can give rise to large-scale collective rearrangements.

During protein elongation, the translocation of mRNA and tRNA molecules across the 30S ribosomal subunit is associated with large-scale motions of the 30S head domain. Here the authors carry out MD simulations to probe the associated steric interactions and identify novel tilting motions during the late stages of translocation.

Information and redundancy in the burial folding code of globular proteins within a wide range of shapes and sizes

Recent ab initio folding simulations for a limited number of small proteins have corroborated a previous suggestion that atomic burial information obtainable from sequence could be sufficient for tertiary structure determination when combined to sequence-independent geometrical constraints. Here, we use simulations parameterized by native burials to investigate the required amount of information in a diverse set of globular proteins comprising different structural classes and a wide size range. Burial information is provided by a potential term pushing each atom towards one among a small number L of equiprobable concentric layers. An upper bound for the required information is provided by the minimal number of layers L(min) still compatible with correct folding behavior. We obtain L(min) between 3 and 5 for seven small to medium proteins with 50 ≤ Nr ≤ 110 residues while for a larger protein with Nr = 141 we find that L ≥ 6 is required to maintain native stability. We additionally estimate the usable redundancy for a given L ≥ L(min) from the burial entropy associated to the largest folding-compatible fraction of "superfluous" atoms, for which the burial term can be turned off or target layers can be chosen randomly. The estimated redundancy for small proteins with L = 4 is close to 0.8. Our results are consistent with the above-average quality of burial predictions used in previous simulations and indicate that the fraction of approachable proteins could increase significantly with even a mild, plausible, improvement on sequence-dependent burial prediction or on sequence-independent constraints that augment the detectable redundancy during simulations.

Using the folding landscapes of proteins to understand protein function

Proteins fold on a biologically-relevant timescale because of a funnel-shaped energy landscape. This landscape is sculpted through evolution by selecting amino-acid sequences that stabilize native interactions while suppressing stable non-native interactions that occur during folding. However, there is strong evolutionary selection for functional residues and these cannot be chosen to optimize folding. Their presence impacts the folding energy landscape in a variety of ways. Here, we survey the effects of functional residues on folding by providing several examples. We then review how such effects can be detected computationally and be used as assays for protein function. Overall, an understanding of how functional residues modulate folding should provide insights into the design of natural proteins and their homeostasis.

Constructing sequence-dependent protein models using coevolutionary information

Recent developments in global statistical methodologies have advanced the analysis of large collections of protein sequences for coevolutionary information. Coevolution between amino acids in a protein arises from compensatory mutations that are needed to maintain the stability or function of a protein over the course of evolution. This gives rise to quantifiable correlations between amino acid sites within the multiple sequence alignment of a protein family. Here, we use the maximum entropy-based approach called mean field Direct Coupling Analysis (mfDCA) to infer a Potts model Hamiltonian governing the correlated mutations in a protein family. We use the inferred pairwise statistical couplings to generate the sequence-dependent heterogeneous interaction energies of a structure-based model (SBM) where only native contacts are considered. Considering the ribosomal S6 protein and its circular permutants as well as the SH3 protein, we demonstrate that these models quantitatively agree with experimental data on folding mechanisms. This work serves as a new framework for generating coevolutionary data-enriched models that can potentially be used to engineer key functional motions and novel interactions in protein systems.

Functionally Relevant Specific Packing Can Determine Protein Folding Routes

Functional residues can modulate the folding mechanisms of proteins. In some proteins, mutations to such residues can radically change the primary folding route. Is it possible then to learn more about the functional regions of a protein by investigating just its choice of folding route? The folding and the function of the protein Escherichia coli ribonuclease H (ecoRNase-H) have been extensively studied and its folding route is known to near-residue resolution. Here, we computationally study the folding of ecoRNase-H using molecular dynamics simulations of structure-based models of increasing complexity. The differences between a model that correctly predicts the experimentally determined folding route and a simpler model that does not can be attributed to a set of six aromatic residues clustered together in a region of the protein called CORE. This clustering, which we term "specific" packing, drives CORE to fold early and determines the folding route. Both the residues involved in specific packing and their packing are largely conserved across E. coli-like RNase-Hs from diverse species. Residue conservation is usually implicated in function. Here, the identified residues either are known to bind substrate in ecoRNase-H or pack against the substrate in the homologous human RNase-H where a substrate-bound crystal structure exists. Thus, the folding mechanism of ecoRNase-H is a byproduct of functional demands upon its sequence. Using our observations on specific packing, we suggest mutations to an engineered HIV RNase-H to make its function better. Our results show that understanding folding route choice in proteins can provide unexpected insights into their function.

Geometrical Frustration in Interleukin-33 Decouples the Dynamics of the Functional Element from the Folding Transition State Ensemble

Interleukin-33 (IL-33) is currently the focus of multiple investigations into targeting pernicious inflammatory disorders. This mediator of inflammation plays a prevalent role in chronic disorders such as asthma, rheumatoid arthritis, and progressive heart disease. In order to better understand the possible link between the folding free energy landscape and functional regions in IL-33, a combined experimental and theoretical approach was applied. IL-33 is a pseudo- symmetrical protein composed of three distinct structural elements that complicate the folding mechanism due to competition for nucleation on the dominant folding route. Trefoil 1 constitutes the majority of the binding interface with the receptor whereas Trefoils 2 and 3 provide the stable scaffold to anchor Trefoil 1. We identified that IL-33 folds with a three-state mechanism, leading to a rollover in the refolding arm of its chevron plots in strongly native conditions. In addition, there is a second slower refolding phase that exhibits the same rollover suggesting similar limitations in folding along parallel routes. Characterization of the intermediate state and the rate limiting steps required for folding suggests that the rollover is attributable to a moving transition state, shifting from a post- to pre-intermediate transition state as you move from strongly native conditions to the midpoint of the transition. On a structural level, we found that initially, all independent Trefoil units fold equally well until a Q_CA of 0.35 when Trefoil 1 will backtrack in order to allow Trefoils 2 and 3 to fold in the intermediate state, creating a stable scaffold for Trefoil 1 to fold onto during the final folding transition. The formation of this intermediate state and subsequent moving transition state is a result of balancing the difficulty in folding the functionally important Trefoil 1 onto the remainder of the protein. Taken together our results indicate that the functional element of the protein is geometrically frustrated, requiring the more stable elements to fold first, acting as a scaffold for docking of the functional element to allow productive folding to the native state.

From residue coevolution to protein conformational ensembles and functional dynamics

The analysis of evolutionary amino acid correlations has recently attracted a surge of renewed interest, also due to their successful use in de novo protein native structure prediction. However, many aspects of protein function, such as substrate binding and product release in enzymatic activity, can be fully understood only in terms of an equilibrium ensemble of alternative structures, rather than a single static structure. In this paper we combine coevolutionary data and molecular dynamics simulations to study protein conformational heterogeneity. To that end, we adapt the Boltzmann-learning algorithm to the analysis of homologous protein sequences and develop a coarse-grained protein model specifically tailored to convert the resulting contact predictions to a protein structural ensemble. By means of exhaustive sampling simulations, we analyze the set of conformations that are consistent with the observed residue correlations for a set of representative protein domains, showing that (i) the most representative structure is consistent with the experimental fold and (ii) the various regions of the sequence display different stability, related to multiple biologically relevant conformations and to the cooperativity of the coevolving pairs. Moreover, we show that the proposed protocol is able to reproduce the essential features of a protein folding mechanism as well as to account for regions involved in conformational transitions through the correct sampling of the involved conformers.

Calcium-controlled conformational choreography in the N-terminal half of adseverin

Adseverin is a member of the calcium-regulated gelsolin superfamily of actin-binding proteins. Here we report the crystal structure of the calcium-free N-terminal half of adseverin (iA1-A3) and the Ca(2+)-bound structure of A3, which reveal structural similarities and differences with gelsolin. Solution small-angle X-ray scattering combined with ensemble optimization revealed a dynamic Ca(2+)-dependent equilibrium between inactive, intermediate and active conformations. Increasing calcium concentrations progressively shift this equilibrium from a main population of inactive conformation to the active form. Molecular dynamics simulations of iA1-A3 provided insights into Ca(2+)-induced destabilization, implicating a critical role for the A2 type II calcium-binding site and the A2A3 linker in the activation process. Finally, mutations that disrupt the A1/A3 interface increase Ca(2+)-independent F-actin severing by A1-A3, albeit at a lower efficiency than observed for gelsolin domains G1-G3. Together, these data address the calcium dependency of A1-A3 activity in relation to the calcium-independent activity of G1-G3.

Prediction of the optimal set of contacts to fold the smallest knotted protein

Knotted protein chains represent a new motif in protein folds. They have been linked to various diseases, and recent extensive analysis of the Protein Data Bank shows that they constitute 1.5% of all deposited protein structures. Despite thorough theoretical and experimental investigations, the role of knots in proteins still remains elusive. Nonetheless, it is believed that knots play an important role in mechanical and thermal stability of proteins. Here, we perform a comprehensive analysis of native, shadow-specific and non-native interactions which describe free energy landscape of the smallest knotted protein (PDB id 2efv). We show that the addition of shadow-specific contacts in the loop region greatly enhances folding kinetics, while the addition of shadow-specific contacts along the C-terminal region (H3 or H4) results in a new folding route with slower kinetics. By means of direct coupling analysis (DCA) we predict non-native contacts which also can accelerate kinetics. Next, we show that the length of the C-terminal knot tail is responsible for the shape of the free energy barrier, while the influence of the elongation of the N-terminus is not significant. Finally, we develop a concept of a minimal contact map sufficient for 2efv protein to fold and analyze properties of this protein using this map.

The ribosome's energy landscape: Recent insights from computation

The ever-increasing capacity of computing resources has extended ribosome calculations from the study of small-scale fluctuations to large-scale barrier-crossing processes. As the field of computational/theoretical biophysics shifts focus to large-scale conformational transitions, there is a growing need for a systematic framework to interpret and analyze ribosome dynamics. To this end, energy landscape principles, largely developed for the study of biomolecular folding, have proven to be invaluable. These tools not only provide a foundation for describing simulations but can be used to reconcile experimental results, as well. In this review, I will discuss recent efforts to employ computational methods to reveal the characteristics of the ribosome's landscape and how these studies can help guide a new generation of experiments that more closely probe the underlying energetics. As a result of these investigations, general principles about ribosome function are beginning to emerge, including that: (1) small-scale fluctuations are the result of structure, rather than detailed energetics, (2) molecular flexibility leads to entropically favored rearrangements, and (3) tRNA dynamics may be accurately described as diffusive movement across an energy landscape.

MD Simulations of tRNA and Aminoacyl-tRNA Synthetases: Dynamics, Folding, Binding, and Allostery

While tRNA and aminoacyl-tRNA synthetases are classes of biomolecules that have been extensively studied for decades, the finer details of how they carry out their fundamental biological functions in protein synthesis remain a challenge. Recent molecular dynamics (MD) simulations are verifying experimental observations and providing new insight that cannot be addressed from experiments alone. Throughout the review, we briefly discuss important historical events to provide a context for how far the field has progressed over the past few decades. We then review the background of tRNA molecules, aminoacyl-tRNA synthetases, and current state of the art MD simulation techniques for those who may be unfamiliar with any of those fields. Recent MD simulations of tRNA dynamics and folding and of aminoacyl-tRNA synthetase dynamics and mechanistic characterizations are discussed. We highlight the recent successes and discuss how important questions can be addressed using current MD simulations techniques. We also outline several natural next steps for computational studies of AARS:tRNA complexes.

Generalized Manning Condensation Model Captures the RNA Ion Atmosphere

RNA is highly sensitive to the ionic environment and typically requires Mg(2+) to form compact structures. There is a need for models capable of describing the ion atmosphere surrounding RNA with quantitative accuracy. We present a model of RNA electrostatics and apply it within coarse-grained molecular dynamics simulation. The model treats Mg(2+) ions explicitly to account for ion-ion correlations neglected by mean-field theories. Since mean-field theories capture KCl well, it is treated implicitly by a generalized Manning counterion condensation model. The model extends Manning condensation to deal with arbitrary RNA conformations, nonlimiting KCl concentrations, and the ion inaccessible volume of RNA. The model is tested against experimental measurements of the excess Mg(2+) associated with the RNA, Γ(2+), because Γ(2+) is directly related to the Mg(2+)-RNA interaction free energy. The excellent agreement with experiment demonstrates that the model captures the ionic dependence of the RNA free energy landscape.

Neither Two-State nor Three-State: Dimerization of Lambda Cro Repressor

Lambda Cro repressor is one of the best studied dimeric transcription factors. However, there has still been an unsettled debate for decades about whether it is a two-state dimer or three-state dimer. We provide a new mechanism model that can reconcile these seemingly conflicting (mutually exclusive) experimental results. From simulations with all-atom structure-based model, we observe that the dimerization process of Lambda Cro repressor starts from one folded monomer with one unfolded monomer. Intrasubunit folding and intersubunit binding are partially coupled, in a fly casting manner.

GENESIS: a hybrid-parallel and multi-scale molecular dynamics simulator with enhanced sampling algorithms for biomolecular and cellular simulations

GENESIS (Generalized-Ensemble Simulation System) is a new software package for molecular dynamics (MD) simulations of macromolecules. It has two MD simulators, called ATDYN and SPDYN. ATDYN is parallelized based on an atomic decomposition algorithm for the simulations of all-atom force-field models as well as coarse-grained Go-like models. SPDYN is highly parallelized based on a domain decomposition scheme, allowing large-scale MD simulations on supercomputers. Hybrid schemes combining OpenMP and MPI are used in both simulators to target modern multicore computer architectures. Key advantages of GENESIS are (1) the highly parallel performance of SPDYN for very large biological systems consisting of more than one million atoms and (2) the availability of various REMD algorithms (T-REMD, REUS, multi-dimensional REMD for both all-atom and Go-like models under the NVT, NPT, NPAT, and NPγT ensembles). The former is achieved by a combination of the midpoint cell method and the efficient three-dimensional Fast Fourier Transform algorithm, where the domain decomposition space is shared in real-space and reciprocal-space calculations. Other features in SPDYN, such as avoiding concurrent memory access, reducing communication times, and usage of parallel input/output files, also contribute to the performance. We show the REMD simulation results of a mixed (POPC/DMPC) lipid bilayer as a real application using GENESIS. GENESIS is released as free software under the GPLv2 licence and can be easily modified for the development of new algorithms and molecular models. WIREs Comput Mol Sci 2015, 5:310–323. doi: 10.1002/wcms.1220

Contribution to the prediction of the fold code: application to immunoglobulin and flavodoxin cases

Background

Folding nucleus of globular proteins formation starts by the mutual interaction of a group of hydrophobic amino acids whose close contacts allow subsequent formation and stability of the 3D structure. These early steps can be predicted by simulation of the folding process through a Monte Carlo (MC) coarse grain model in a discrete space. We previously defined MIRs (Most Interacting Residues), as the set of residues presenting a large number of non-covalent neighbour interactions during such simulation. MIRs are good candidates to define the minimal number of residues giving rise to a given fold instead of another one, although their proportion is rather high, typically [15-20]% of the sequences. Having in mind experiments with two sequences of very high levels of sequence identity (up to 90%) but different folds, we combined the MIR method, which takes sequence as single input, with the “fuzzy oil drop” (FOD) model that requires a 3D structure, in order to estimate the residues coding for the fold. FOD assumes that a globular protein follows an idealised 3D Gaussian distribution of hydrophobicity density, with the maximum in the centre and minima at the surface of the “drop”. If the actual local density of hydrophobicity around a given amino acid is as high as the ideal one, then this amino acid is assigned to the core of the globular protein, and it is assumed to follow the FOD model. Therefore one obtains a distribution of the amino acids of a protein according to their agreement or rejection with the FOD model.

Results

We compared and combined MIR and FOD methods to define the minimal nucleus, or keystone, of two populated folds: immunoglobulin-like (Ig) and flavodoxins (Flav). The combination of these two approaches defines some positions both predicted as a MIR and assigned as accordant with the FOD model. It is shown here that for these two folds, the intersection of the predicted sets of residues significantly differs from random selection. It reduces the number of selected residues by each individual method and allows a reasonable agreement with experimentally determined key residues coding for the particular fold. In addition, the intersection of the two methods significantly increases the specificity of the prediction, providing a robust set of residues that constitute the folding nucleus.

Exploring the balance between folding and functional dynamics in proteins and RNA

As our understanding of biological dynamics continues to be refined, it is becoming clear that biomolecules can undergo transitions between ordered and disordered states as they execute functional processes. From a computational perspective, studying disorder events poses a challenge, as they typically occur on long timescales, and the associated molecules are often large (i.e., hundreds of residues). These size and time requirements make it advantageous to use computationally inexpensive models to characterize large-scale dynamics, where more highly detailed models can provide information about individual sub-steps associated with function. To reduce computational demand, one often uses a coarse-grained representation of the molecule or a simplified description of the energetics. In order to use simpler models to identify transient disorder in RNA and proteins, it is imperative that these models can accurately capture structural fluctuations about folded configurations, as well as the overall stability of each molecule. Here, we explore a class of simplified model for which all non-hydrogen atoms are explicitly represented. We find that this model can provide a consistent description of protein folding and native-basin dynamics for several representative biomolecules. We additionally show that the native-basin fluctuations of tRNA and the ribosome are robust to variations in the model. Finally, the extended variable loop in tRNAIle is predicted to be very dynamic, which may facilitate biologically-relevant rearrangements. Together, this study provides a foundation that will aid in the application of simplified models to study disorder during function in ribonucleoprotein (RNP) assemblies.

Conformational changes in a hyperthermostable glycoside hydrolase: enzymatic activity is a consequence of the loop dynamics and protonation balance

Endo-β-1, 4-mannanase from Thermotoga petrophila (TpMan) is a modular hyperthermostable enzyme involved in the degradation of mannan-containing polysaccharides. The degradation of these polysaccharides represents a key step for several industrial applications. Here, as part of a continuing investigation of TpMan, the region corresponding to the GH5 domain (TpManGH5) was characterized as a function of pH and temperature. The results indicated that the enzymatic activity of the TpManGH5 is pH-dependent, with its optimum activity occurring at pH 6. At pH 8, the studies demonstrated that TpManGH5 is a molecule with a nearly spherical tightly packed core displaying negligible flexibility in solution, and with size and shape very similar to crystal structure. However, TpManGH5 experiences an increase in radius of gyration in acidic conditions suggesting expansion of the molecule. Furthermore, at acidic pH values, TpManGH5 showed a less globular shape, probably due to a loop region slightly more expanded and flexible in solution (residues Y88 to A105). In addition, molecular dynamics simulations indicated that conformational changes caused by pH variation did not change the core of the TpManGH5, which means that only the above mentioned loop region presents high degree of fluctuations. The results also suggested that conformational changes of the loop region may facilitate polysaccharide and enzyme interaction. Finally, at pH 6 the results indicated that TpManGH5 is slightly more flexible at 65°C when compared to the same enzyme at 20°C. The biophysical characterization presented here is well correlated with the enzymatic activity and provide new insight into the structural basis for the temperature and pH-dependent activity of the TpManGH5. Also, the data suggest a loop region that provides a starting point for a rational design of biotechnological desired features.

Integrating molecular dynamics simulations with chemical probing experiments using SHAPE-FIT

Integration and calibration of molecular dynamics simulations with experimental data remain a challenging endeavor. We have developed a novel method to integrate chemical probing experiments with molecular simulations of RNA molecules by using a native structure-based model. Selective 2'-hydroxyl acylation by primer extension (SHAPE) characterizes the mobility of each residue in the RNA. Our method, SHAPE-FIT, automatically optimizes the potential parameters of the force field according to measured reactivities from SHAPE. The optimized parameter set allows simulations of dynamics highly consistent with SHAPE probing experiments. Such atomistic simulations, thoroughly grounded in experiment, can open a new window on RNA structure-function relations.

What protein folding teaches us about biological function and molecular machines

Protein folding was the first area of molecular biology for which a systematic statistical-mechanical analysis of dynamics was developed. As a result, folding is described as a process by which a disordered protein chain diffuses across a high-dimensional energy landscape and finally reaches the folded ensemble. Folding studies have produced countless theoretical concepts that are generalizable to other biomolecular processes, such as the functional dynamics of molecular assemblies. Common themes in folding and function include the dominant role of excluded volume, that a balance between energetic roughness and geometrical effects guides dynamics, and that folding/functional landscapes are relatively smooth. Here, we discuss how insights into protein folding have been applied to investigate the functional dynamics of biomolecular assemblies.

Folding mechanisms of Trefoil Knot proteins studied by molecular dynamics simulations and Go-model

Most proteins need to avoid the complex topologies when folding into the native structures, but some proteins with nontrivial topologies have been found in nature. Here we used protein unfolding simulations under high temperature and all-atom Gō-model to investigate the folding mechanisms for two trefoil knot proteins. Results show that, the contacts in β-sheet are important to the formation of knot protein, and if these contacts disappeared, the knot protein would be easy to untie. In the Gō-model simulations, the folding processes of the two knot proteins are similar. The compact structures of the two knot proteins with the native contacts in β-sheet are formed in transition state, and the intermediate state has loose C-terminal. This model also reveals the detailed folding mechanisms for the two proteins.

Folding and self-assembly of a small heterotetramer

Designed miniproteins offer a possibility to study folding and association of protein complexes, both experimentally and in silico. Using replica exchange molecular dynamics and the coarse-grain UNRES force field, we have simulated the folding and self-assembly of the heterotetramer BBAThet1, comparing it with that of the homotetramer BBAT1 which has the same size and ββα-fold. For both proteins, association of the tetramer precedes and facilitates folding of the individual chains.

Native structure-based modeling and simulation of biomolecular systems per mouse click

BACKGROUND

Molecular dynamics (MD) simulations provide valuable insight into biomolecular systems at the atomic level. Notwithstanding the ever-increasing power of high performance computers current MD simulations face several challenges: the fastest atomic movements require time steps of a few femtoseconds which are small compared to biomolecular relevant timescales of milliseconds or even seconds for large conformational motions. At the same time, scalability to a large number of cores is limited mostly due to long-range interactions. An appealing alternative to atomic-level simulations is coarse-graining the resolution of the system or reducing the complexity of the Hamiltonian to improve sampling while decreasing computational costs. Native structure-based models, also called Gō-type models, are based on energy landscape theory and the principle of minimal frustration. They have been tremendously successful in explaining fundamental questions of, e.g., protein folding, RNA folding or protein function. At the same time, they are computationally sufficiently inexpensive to run complex simulations on smaller computing systems or even commodity hardware. Still, their setup and evaluation is quite complex even though sophisticated software packages support their realization.

RESULTS

Here, we establish an efficient infrastructure for native structure-based models to support the community and enable high-throughput simulations on remote computing resources via GridBeans and UNICORE middleware. This infrastructure organizes the setup of such simulations resulting in increased comparability of simulation results. At the same time, complete workflows for advanced simulation protocols can be established and managed on remote resources by a graphical interface which increases reusability of protocols and additionally lowers the entry barrier into such simulations for, e.g., experimental scientists who want to compare their results against simulations. We demonstrate the power of this approach by illustrating it for protein folding simulations for a range of proteins.

CONCLUSIONS

We present software enhancing the entire workflow for native structure-based simulations including exception-handling and evaluations. Extending the capability and improving the accessibility of existing simulation packages the software goes beyond the state of the art in the domain of biomolecular simulations. Thus we expect that it will stimulate more individuals from the community to employ more confidently modeling in their research.

Order and disorder control the functional rearrangement of influenza hemagglutinin

Influenza hemagglutinin (HA), a homotrimeric glycoprotein crucial for membrane fusion, undergoes a large-scale structural rearrangement during viral invasion. X-ray crystallography has shown that the pre- and postfusion configurations of HA2, the membrane-fusion subunit of HA, have disparate secondary, tertiary, and quaternary structures, where some regions are displaced by more than 100 Å. To explore structural dynamics during the conformational transition, we studied simulations of a minimally frustrated model based on energy landscape theory. The model combines structural information from both the pre- and postfusion crystallographic configurations of HA2. Rather than a downhill drive toward formation of the central coiled-coil, we discovered an order-disorder transition early in the conformational change as the mechanism for the release of the fusion peptides from their burial sites in the prefusion crystal structure. This disorder quickly leads to a metastable intermediate with a broken threefold symmetry. Finally, kinetic competition between the formation of the extended coiled-coil and C-terminal melting results in two routes from this intermediate to the postfusion structure. Our study reiterates the roles that cracking and disorder can play in functional molecular motions, in contrast to the downhill mechanical interpretations of the "spring-loaded" model proposed for the HA2 conformational transition.

Pierced Lasso Bundles are a new class of knot-like motifs

A four-helix bundle is a well-characterized motif often used as a target for designed pharmaceutical therapeutics and nutritional supplements. Recently, we discovered a new structural complexity within this motif created by a disulphide bridge in the long-chain helical bundle cytokine leptin. When oxidized, leptin contains a disulphide bridge creating a covalent-loop through which part of the polypeptide chain is threaded (as seen in knotted proteins). We explored whether other proteins contain a similar intriguing knot-like structure as in leptin and discovered 11 structurally homologous proteins in the PDB. We call this new helical family class the Pierced Lasso Bundle (PLB) and the knot-like threaded structural motif a Pierced Lasso (PL). In the current study, we use structure-based simulation to investigate the threading/folding mechanisms for all the PLBs along with three unthreaded homologs as the covalent loop (or lasso) in leptin is important in folding dynamics and activity. We find that the presence of a small covalent loop leads to a mechanism where structural elements slipknot to thread through the covalent loop. Larger loops use a piercing mechanism where the free terminal plugs through the covalent loop. Remarkably, the position of the loop as well as its size influences the native state dynamics, which can impact receptor binding and biological activity. This previously unrecognized complexity of knot-like proteins within the helical bundle family comprises a completely new class within the knot family, and the hidden complexity we unraveled in the PLBs is expected to be found in other protein structures outside the four-helix bundles. The insights gained here provide critical new elements for future investigation of this emerging class of proteins, where function and the energetic landscape can be controlled by hidden topology, and should be take into account in ab initio predictions of newly identified protein targets.

Structure-based simulations reveal concerted dynamics of GPCR activation

G protein-coupled receptors (GPCRs) are a vital class of proteins that transduce biological signals across the cell membrane. However, their allosteric activation mechanism is not fully understood; crystal structures of active and inactive receptors have been reported, but the functional pathway between these two states remains elusive. Here, we use structure-based (Gō-like) models to simulate activation of two GPCRs, rhodopsin and the β₂ adrenergic receptor (β₂AR). We used data-derived reaction coordinates that capture the activation mechanism for both proteins, showing that activation proceeds through quantitatively different paths in the two systems. Both reaction coordinates are determined from the dominant concerted motions in the simulations so the technique is broadly applicable. There were two surprising results. First, the main structural changes in the simulations were distributed throughout the transmembrane bundle, and not localized to the obvious areas of interest, such as the intracellular portion of Helix 6. Second, the activation (and deactivation) paths were distinctly nonmonotonic, populating states that were not simply interpolations between the inactive and active structures. These transitions also suggest a functional explanation for β₂AR's basal activity: it can proceed through a more broadly defined path during the observed transitions.

Drastic changes in conformational dynamics of the antiterminator M2-1 regulate transcription efficiency in Pneumovirinae

The M2-1 protein of human metapneumovirus (HMPV) is a zinc-binding transcription antiterminator which is highly conserved among pneumoviruses. We report the structure of tetrameric HMPV M2-1. Each protomer features a N-terminal zinc finger domain and an α-helical tetramerization motif forming a rigid unit, followed by a flexible linker and an α-helical core domain. The tetramer is asymmetric, three of the protomers exhibiting a closed conformation, and one an open conformation. Molecular dynamics simulations and SAXS demonstrate a dynamic equilibrium between open and closed conformations in solution. Structures of adenosine monophosphate- and DNA- bound M2-1 establish the role of the zinc finger domain in base-specific recognition of RNA. Binding to 'gene end' RNA sequences stabilized the closed conformation of M2-1 leading to a drastic shift in the conformational landscape of M2-1. We propose a model for recognition of gene end signals and discuss the implications of these findings for transcriptional regulation in pneumoviruses.DOI: http://dx.doi.org/10.7554/eLife.02674.001.

Simulating movement of tRNA through the ribosome during hybrid-state formation

Biomolecular simulations provide a means for exploring the relationship between flexibility, energetics, structure, and function. With the availability of atomic models from X-ray crystallography and cryoelectron microscopy (cryo-EM), and rapid increases in computing capacity, it is now possible to apply molecular dynamics (MD) simulations to large biomolecular machines, and systematically partition the factors that contribute to function. A large biomolecular complex for which atomic models are available is the ribosome. In the cell, the ribosome reads messenger RNA (mRNA) in order to synthesize proteins. During this essential process, the ribosome undergoes a wide range of conformational rearrangements. One of the most poorly understood transitions is translocation: the process by which transfer RNA (tRNA) molecules move between binding sites inside of the ribosome. The first step of translocation is the adoption of a "hybrid" configuration by the tRNAs, which is accompanied by large-scale rotations in the ribosomal subunits. To illuminate the relationship between these rearrangements, we apply MD simulations using a multi-basin structure-based (SMOG) model, together with targeted molecular dynamics protocols. From 120 simulated transitions, we demonstrate the viability of a particular route during P/E hybrid-state formation, where there is asynchronous movement along rotation and tRNA coordinates. These simulations not only suggest an ordering of events, but they highlight atomic interactions that may influence the kinetics of hybrid-state formation. From these simulations, we also identify steric features (H74 and surrounding residues) encountered during the hybrid transition, and observe that flexibility of the single-stranded 3'-CCA tail is essential for it to reach the endpoint. Together, these simulations provide a set of structural and energetic signatures that suggest strategies for modulating the physical-chemical properties of protein synthesis by the ribosome.

Reduced model captures Mg(2+)-RNA interaction free energy of riboswitches

The stability of RNA tertiary structures depends heavily on Mg(2+). The Mg(2+)-RNA interaction free energy that stabilizes an RNA structure can be computed experimentally through fluorescence-based assays that measure Γ2+, the number of excess Mg(2+) associated with an RNA molecule. Previous explicit-solvent simulations predict that the majority of excess Mg(2+) ions interact closely and strongly with the RNA, unlike monovalent ions such as K(+), suggesting that an explicit treatment of Mg(2+) is important for capturing RNA dynamics. Here we present a reduced model that accurately reproduces the thermodynamics of Mg(2+)-RNA interactions. This model is able to characterize long-timescale RNA dynamics coupled to Mg(2+) through the explicit representation of Mg(2+) ions. KCl is described by Debye-Hückel screening and a Manning condensation parameter, which represents condensed K(+) and models its competition with condensed Mg(2+). The model contains one fitted parameter, the number of condensed K(+) ions in the absence of Mg(2+). Values of Γ2+ computed from molecular dynamics simulations using the model show excellent agreement with both experimental data on the adenine riboswitch and previous explicit-solvent simulations of the SAM-I riboswitch. This agreement confirms the thermodynamic accuracy of the model via the direct relation of Γ2+ to the Mg(2+)-RNA interaction free energy, and provides further support for the predictions from explicit-solvent calculations. This reduced model will be useful for future studies of the interplay between Mg(2+) and RNA dynamics.

From structure to function: the convergence of structure based models and co-evolutionary information

Understanding protein folding and function is one of the most important problems in biological research. Energy landscape theory and the folding funnel concept have provided a framework to investigate the mechanisms associated to these processes. Since protein energy landscapes are in most cases minimally frustrated, structure based models (SMBs) have successfully determined the geometrical features associated with folding and functional transitions. However, structural information is limited, particularly with respect to different functional configurations. This is a major limitation for SBMs. Alternatively, statistical methods to study amino acid co-evolution provide information on residue-residue interactions useful for the study of structure and function. Here, we show how the combination of these two methods gives rise to a novel way to investigate the mechanisms associated with folding and function. We use this methodology to explore the mechanistic aspects of protein translocation in the integral membrane protease FtsH. Dual basin-SBM simulations using the open and closed state of this hexameric motor reveals a functionally important paddling motion in the catalytic cycle. We also find that Direct Coupling Analysis (DCA) predicts physical contacts between AAA and peptidase domains of the motor, which are crucial for the open to close transition. Our combined method, which uses structural information from the open state experimental structure and co-evolutionary couplings, suggests that this methodology can be used to explore the functional landscape of complex biological macromolecules previously inaccessible to methods dependent on experimental structural information. This efficient way to sample the conformational space of large systems creates a theoretical/computational framework capable of better characterizing the functional landscape in large biomolecular assemblies.

Toward rationally redesigning bacterial two-component signaling systems using coevolutionary information

A challenge in molecular biology is to distinguish the key subset of residues that allow two-component signaling (TCS) proteins to recognize their correct signaling partner such that they can transiently bind and transfer signal, i.e., phosphoryl group. Detailed knowledge of this information would allow one to search sequence space for mutations that can be used to systematically tune the signal transmission between TCS partners as well as potentially encode a TCS protein to preferentially transfer signals to a nonpartner. Motivated by the notion that this detailed information is found in sequence data, we explore the sequence coevolution between signaling partners to better understand how mutations can positively or negatively alter their ability to transfer signal. Using direct coupling analysis for determining evolutionarily conserved protein-protein interactions, we apply a metric called the direct information score to quantify mutational changes in the interaction between TCS proteins and demonstrate that it accurately correlates with experimental mutagenesis studies probing the mutational change in measured in vitro phosphotransfer. Furthermore, by subtracting from our metric an appropriate null model corresponding to generic, conserved features in TCS signaling pairs, we can isolate the determinants that give rise to interaction specificity and recognition, which are variable among different TCS partners. Our methodology forms a potential framework for the rational design of TCS systems by allowing one to quickly search sequence space for mutations or even entirely new sequences that can increase or decrease our metric, as a proxy for increasing or decreasing phosphotransfer ability between TCS proteins.