Bayesian Energy Landscape Tilting: Towards Concordant Models of Molecular Ensembles

Predicting biological structure has remained challenging for systems such as disordered proteins that take on myriad conformations. Hybrid simulation/experiment strategies have been undermined by difficulties in evaluating errors from computational model inaccuracies and data uncertainties. Building on recent proposals from maximum entropy theory and nonequilibrium thermodynamics, we address these issues through a Bayesian energy landscape tilting (BELT) scheme for computing Bayesian hyperensembles over conformational ensembles. BELT uses Markov chain Monte Carlo to directly sample maximum-entropy conformational ensembles consistent with a set of input experimental observables. To test this framework, we apply BELT to model trialanine, starting from disagreeing simulations with the force fields ff96, ff99, ff99sbnmr-ildn, CHARMM27, and OPLS-AA. BELT incorporation of limited chemical shift and ³J measurements gives convergent values of the peptide’s α, β, and PP_II conformational populations in all cases. As a test of predictive power, all five BELT hyperensembles recover set-aside measurements not used in the fitting and report accurate errors, even when starting from highly inaccurate simulations. BELT’s principled framework thus enables practical predictions for complex biomolecular systems from discordant simulations and sparse data.     

Predicting molecular properties of α-synuclein using force fields for intrinsically disordered proteins

Independent force field validation is an essential practice to keep track of developments and for performing meaningful Molecular Dynamics simulations. In this work, atomistic force fields for intrinsically disordered proteins (IDP) are tested by simulating the archetypical IDP α-synuclein in solution for 2.5 μs. Four combinations of protein and water force fields were tested: ff19SB/OPC , ff19SB/TIP4P-D , ff03CMAP/TIP4P-D , and a99SB-disp/TIP4P-disp , with four independent repeat simulations for each combination. We compare our simulations to the results of a 73 μs simulation using the a99SB-disp/TIP4P-disp combination, provided by D. E. Shaw Research. From the trajectories, we predict a range of experimental observations of α-synuclein and compare them to literature data. This includes protein radius of gyration and hydration, intramolecular distances, NMR chemical shifts, and ³J-couplings. Both ff19SB/TIP4P-D and a99SB-disp/TIP4P-disp produce extended conformational ensembles of α-synuclein that agree well with experimental radius of gyration and intramolecular distances while a99SB-disp/TIP4P-disp reproduces a balanced α-synuclein secondary structure content. It was found that ff19SB/OPC and ff03CMAP/TIP4P-D produce overly compact conformational ensembles and show discrepancies in the secondary structure content compared to the experimental data.     

Flexibility of the Bacterial Chaperone Trigger Factor in Microsecond-Timescale Molecular Dynamics Simulations

The bacterial chaperone trigger factor (TF) is the first chaperone to be encountered by a nascent protein chain as it emerges from the ribosome exit tunnel. Experimental results suggest that TF possesses considerable conformational flexibility, and in an attempt to provide an atomic-level view of this flexibility, we have performed independent 1.5-μs molecular dynamics simulations of TF in explicit solvent using two different simulation force fields (OPLS-AA/L and AMBER ff99SB-ILDN). Both simulations indicate that TF possesses tremendous flexibility, with huge excursions from the crystallographic conformation caused by reorientations of the protein’s constituent domains; both simulations also predict the formation of extensive contacts between TF’s PPIase domain and the Arm 1 domain that is involved in nascent-chain binding. In the OPLS simulation, however, TF rapidly settles into a very compact conformation that persists for at least 1 μs, whereas in the AMBER simulation, it remains highly dynamic; additional simulations in which the two force fields were swapped suggest that these differences are at least partly attributable to sampling issues. The simulation results provide potential rationalizations of a number of experimental observations regarding TF’s conformational behavior and have implications for using simulations to model TF’s function on translating ribosomes.     

Effects of force fields on the conformational and dynamic properties of amyloid β(1‐40) dimer explored by replica exchange molecular dynamics simulations

The conformational space and structural ensembles of amyloid beta (Aβ) peptides and their oligomers in solution are inherently disordered and proven to be challenging to study. Optimum force field selection for molecular dynamics (MD) simulations and the biophysical relevance of results are still unknown. We compared the conformational space of the Aβ(1‐40) dimers by 300 ns replica exchange MD simulations at physiological temperature (310 K) using: the AMBER‐ff99sb‐ILDN, AMBER‐ff99sb*‐ILDN, AMBER‐ff99sb‐NMR, and CHARMM22* force fields. Statistical comparisons of simulation results to experimental data and previously published simulations utilizing the CHARMM22* and CHARMM36 force fields were performed. All force fields yield sampled ensembles of conformations with collision cross sectional areas for the dimer that are statistically significantly larger than experimental results. All force fields, with the exception of AMBER‐ff99sb‐ILDN (8.8 ± 6.4%) and CHARMM36 (2.7 ± 4.2%), tend to overestimate the α‐helical content compared to experimental CD (5.3 ± 5.2%). Using the AMBER‐ff99sb‐NMR force field resulted in the greatest degree of variance (41.3 ± 12.9%). Except for the AMBER‐ff99sb‐NMR force field, the others tended to under estimate the expected amount of β‐sheet and over estimate the amount of turn/bend/random coil conformations. All force fields, with the exception AMBER‐ff99sb‐NMR, reproduce a theoretically expected β‐sheet‐turn‐β‐sheet conformational motif, however, only the CHARMM22* and CHARMM36 force fields yield results compatible with collapse of the central and C‐terminal hydrophobic cores from residues 17‐21 and 30‐36. Although analyses of essential subspace sampling showed only minor variations between force fields, secondary structures of lowest energy conformers are different.     

Evaluating Molecular Mechanical Potentials for Helical Peptides and Proteins

Multiple variants of the AMBER all-atom force field were quantitatively evaluated with respect to their ability to accurately characterize helix-coil equilibria in explicit solvent simulations. Using a global distributed computing network, absolute conformational convergence was achieved for large ensembles of the capped A₂₁ and F_s helical peptides. Further assessment of these AMBER variants was conducted via simulations of a flexible 164-residue five-helix-bundle protein, apolipophorin-III, on the 100 ns timescale. Of the contemporary potentials that had not been assessed previously, the AMBER-99SB force field showed significant helix-destabilizing tendencies, with beta bridge formation occurring in helical peptides, and unfolding of apolipophorin-III occurring on the tens of nanoseconds timescale. The AMBER-03 force field, while showing adequate helical propensities for both peptides and stabilizing apolipophorin-III, (i) predicts an unexpected decrease in helicity with ALA→ARG⁺ substitution, (ii) lacks experimentally observed 3₁₀ helical content, and (iii) deviates strongly from average apolipophorin-III NMR structural properties. As is observed for AMBER-99SB, AMBER-03 significantly overweighs the contribution of extended and polyproline backbone configurations to the conformational equilibrium. In contrast, the AMBER-99φ force field, which was previously shown to best reproduce experimental measurements of the helix-coil transition in model helical peptides, adequately stabilizes apolipophorin-III and yields both an average gyration radius and polar solvent exposed surface area that are in excellent agreement with the NMR ensemble.     

Residue-specific α-helix propensities from molecular simulation

Formation of α-helices is a fundamental process in protein folding and assembly. By studying helix formation in molecular simulations of a series of alanine-based peptides, we obtain the temperature-dependent α-helix propensities of all 20 naturally occurring residues with two recent additive force fields, Amber ff03w and Amber ff99SB¹. Encouragingly, we find that the overall helix propensity of many residues is captured well by both energy functions, with Amber ff99SB¹ being more accurate. Nonetheless, there are some residues that deviate considerably from experiment, which can be attributed to two aspects of the energy function: i), variations of the charge model used to determine the atomic partial charges, with residues whose backbone charges differ most from alanine tending to have the largest error; ii), side-chain torsion potentials, as illustrated by the effect of modifications to the torsion angles of I, L, D, N. We find that constrained refitting of residue charges for charged residues in Amber ff99SB¹ significantly improves their helix propensity. The resulting parameters should more faithfully reproduce helix propensities in simulations of protein folding and disordered proteins.     

Combining NMR ensembles and molecular dynamics simulations provides more realistic models of protein structures in solution and leads to better chemical shift prediction

While chemical shifts are invaluable for obtaining structural information from proteins, they also offer one of the rare ways to obtain information about protein dynamics. A necessary tool in transforming chemical shifts into structural and dynamic information is chemical shift prediction. In our previous work we developed a method for 4D prediction of protein ¹H chemical shifts in which molecular motions, the 4th dimension, were modeled using molecular dynamics (MD) simulations. Although the approach clearly improved the prediction, the X-ray structures and single NMR conformers used in the model cannot be considered fully realistic models of protein in solution. In this work, NMR ensembles (NMRE) were used to expand the conformational space of proteins (e.g. side chains, flexible loops, termini), followed by MD simulations for each conformer to map the local fluctuations. Compared with the non-dynamic model, the NMRE+MD model gave 6–17% lower root-mean-square (RMS) errors for different backbone nuclei. The improved prediction indicates that NMR ensembles with MD simulations can be used to obtain a more realistic picture of protein structures in solutions and moreover underlines the importance of short and long time-scale dynamics for the prediction. The RMS errors of the NMRE+MD model were 0.24, 0.43, 0.98, 1.03, 1.16 and 2.39 ppm for ¹Hα, ¹HN, ¹³Cα, ¹³Cβ, ¹³CO and backbone ¹⁵N chemical shifts, respectively. The model is implemented in the prediction program 4DSPOT, available at .     

Evaluating amber force fields using computed NMR chemical shifts

NMR chemical shifts can be computed from molecular dynamics (MD) simulations using a template matching approach and a library of conformers containing chemical shifts generated from ab initio quantum calculations. This approach has potential utility for evaluating the force fields that underlie these simulations. Imperfections in force fields generate flawed atomic coordinates. Chemical shifts obtained from flawed coordinates have errors that can be traced back to these imperfections. We use this approach to evaluate a series of AMBER force fields that have been refined over the course of two decades (ff94, ff96, ff99SB, ff14SB, ff14ipq, and ff15ipq). For each force field a series of MD simulations are carried out for eight model proteins. The calculated chemical shifts for the ¹H, ¹⁵N, and ¹³C^a atoms are compared with experimental values. Initial evaluations are based on root mean squared (RMS) errors at the protein level. These results are further refined based on secondary structure and the types of atoms involved in nonbonded interactions. The best chemical shift for identifying force field differences is the shift associated with peptide protons. Examination of the model proteins on a residue by residue basis reveals that force field performance is highly dependent on residue position. Examination of the time course of nonbonded interactions at these sites provides explanations for chemical shift differences at the atomic coordinate level. Results show that the newer ff14ipq and ff15ipq force fields developed with the implicitly polarized charge method perform better than the older force fields.     

A molecular dynamics simulation investigation of the relative stability of the cyclic peptide octreotide and its deprotonated and its (CF3)-Trp substituted analogs in different solvents

The cyclic octa-peptide octreotide and its derivatives are used as diagnostics and therapeutics in relation to particular types of cancers. This led to investigations of their conformational properties using spectroscopic, NMR and CD, methods. A CF₃-substituted derivative, that was designed to stabilize the dominant octreotide conformer responsible for receptor binding, turned out to have a lower affinity. The obtained spectroscopic data were interpreted as to show an increased flexibility of the CF₃ derivative compared to the unsubstituted octreotide, which could then explain the lower affinity.In this article, we use MD simulation without and with time-averaged NOE distance and time-averaged local-elevation ³J-coupling restraining representing experimental NMR data to determine the conformational properties of the different peptides in the different solvents for which experimental data are available, that are compatible with the NOE atom–atom distance bounds and the ³J_HNHα-couplings as derived from the NMR measurements. The conformational ensembles show that the CF₃ substitution in combination with the change of solvent from water to methanol leads to a decrease in flexibility and a shift in the populations of the dominant conformers that are compatible with the experimental data.     

Analysis of Sub-τc and Supra-τc Motions in Protein Gβ1 Using Molecular Dynamics Simulations

The functions of proteins depend on the dynamical behavior of their native states on a wide range of timescales. To investigate these dynamics in the case of the small protein Gβ1, we analyzed molecular dynamics simulations with the model-free approach of nuclear magnetic relaxation. We found amplitudes of fast timescale motions (sub-τ_c, where τ_c is the rotational correlation time) consistent with S² obtained from spin relaxation measurements as well as amplitudes of slow timescale motions (supra-τ_c) in quantitative agreement with S² order parameters derived from residual dipolar coupling measurements. The slow timescale motions are associated with the large variations of the ³J couplings that follow transitions between different conformational substates. These results provide further characterization of the large structural fluctuations in the native states of proteins that occur on timescales longer than the rotational correlation time.     

GROMOS96 43a1 performance in predicting oligosaccharide conformational ensembles within glycoproteins

In previous work [Pol-Fachin, L.; Fernandes, C. L.; Verli, H.; Carbohydr. Res.2009, 344, 491-500], we had demonstrated that GROMOS96 43a1 force field and Löwdin HF/6-31G∗∗-derived atomic charges, adequately represent a glycoprotein’s conformational ensemble in aqueous solutions, taking as the starting geometries NMR-determined structures. Based on such data, the present work intends to evaluate the use of the main solution conformations of isolated disaccharides, to build the carbohydrate moiety of glycoproteins, for which no previous experimental information is available. The observed results suggested that the entire glycoprotein scaffold appears unable to promote major modifications in the conformational behavior of glycosidic linkages. Additionally, when compared to energy contour plots, the results support the use of solution ensembles, to refine vacuum conformations of carbohydrate databases in the assembling of glycoproteins 3D structures. Finally, such approach is applied to build a full glycosylated model for COX-1 and COX-2 enzymes.     

Microsecond simulations of mdm2 and its complex with p53 yield insight into force field accuracy and conformational dynamics

The p53‐MDM2 complex is both a major target for cancer drug development and a valuable model system for computational predictions of protein‐ligand binding. To investigate the accuracy of molecular simulations of MDM2 and its complex with p53, we performed a number of long (200 ns to 1 µs) explicit‐solvent simulations using a range of force fields. We systematically compared nine popular force fields (AMBER ff03, ff12sb, ff14sb, ff99sb, ff99sb‐ildn, ff99sb‐ildn‐nmr, ff99sb‐ildn‐phi, CHARMM22*, and CHARMM36) against experimental chemical shift data, and found similarly accurate results, with microsecond simulations achieving better agreement compared to 200‐ns trajectories. Although the experimentally determined apo structure has a closed binding cleft, simulations in all force fields suggest the apo state of MDM2 is highly flexible, and able to sample holo‐like conformations, consistent with a conformational selection model. Initial structuring of the MDM2 lid region, known to competitively bind the binding cleft, is also observed in long simulations. Taken together, these results show molecular simulations can accurately sample conformations relevant for ligand binding. We expect this study to inform future computational work on folding and binding of MDM2 ligands. Proteins 2015; 83:1665–1676. © 2015 Wiley Periodicals, Inc.     

Structural Interpretation of Paramagnetic Relaxation Enhancement-Derived Distances for Disordered Protein States

Paramagnetic relaxation enhancement (PRE) is a powerful technique for studying transient tertiary organizations of unfolded and partially folded proteins. The heterogeneous and dynamic nature of disordered protein states, together with the r⁻⁶ dependence of PRE, presents significant challenges for reliable structural interpretation of PRE-derived distances. Without additional knowledge of accessible conformational substates, ensemble-simulation-based protocols have been used to calculate structure ensembles that appear to be consistent with the PRE distance restraints imposed on the ensemble level with the proper r⁻⁶ weighting. However, rigorous assessment of the reliability of such protocols has been difficult without intimate knowledge of the true nature of disordered protein states. Here we utilize sets of theoretical PRE distances derived from simulated structure ensembles that represent the folded, partially folded and unfolded states of a small protein to investigate the efficacy of ensemble-simulation-based structural interpretation of PRE distances. The results confirm a critical limitation that, due to r⁻⁶ weighting, only one or a few members need to satisfy the distance restraints and the rest of the ensemble are essentially unrestrained. Consequently, calculated structure ensembles will appear artificially heterogeneous no matter whether the PRE distances are derived from the folded, partially unfolded or unfolded state. Furthermore, the nature of the heterogeneous ensembles is largely determined by the protein model employed in structure calculation and reflects little on the true nature of the underlying disordered state. These findings suggest that PRE measurements on disordered protein states alone generally do not contain enough information for a reliable structural interpretation and that the latter will require additional knowledge of accessible conformational substates. Interestingly, when a very large number of PRE measurements is available, faithful structural interpretation might be possible with intermediate ensemble sizes under ideal conditions.     

Structure and Disorder in an Unfolded State under Nondenaturing Conditions from Ensemble Models Consistent with a Large Number of Experimental Restraints

Obtaining detailed structural models of disordered states of proteins under nondenaturing conditions is important for a better understanding of both functional intrinsically disordered proteins and unfolded states of folded proteins. Extensive experimental characterization of the drk N-terminal SH3 domain unfolded state has shown that, although it appears to be highly disordered, it possesses significant nonrandom secondary and tertiary structure. In our previous attempts to generate structural models of the unfolded state using the program ENSEMBLE, we were limited by insufficient experimental restraints and conformational sampling. In this study, we have vastly expanded our experimental restraint set to include ¹H-¹⁵N residual dipolar couplings, small-angle X-ray scattering measurements, nitroxide paramagnetic relaxation enhancements, O₂-induced ¹³C paramagnetic shifts, hydrogen-exchange protection factors, and ¹⁵N R₂ data, in addition to the previously used nuclear Overhauser effects, amino terminal Cu²⁺-Ni²⁺ binding paramagnetic relaxation enhancements, J-couplings, chemical shifts, hydrodynamic radius, and solvent accessibility restraints. We have also implemented a new ensemble calculation methodology that uses iterative conformational sampling and seeks to calculate the simplest possible ensemble models. As a result, we can now generate ensembles that are consistent with much larger experimental data sets than was previously possible. Although highly heterogeneous and having broad molecular size distributions, the calculated drk N-terminal SH3 domain unfolded-state ensembles have very different properties than expected for random or statistical coils and possess significant nonnative α-helical structure and both native-like and nonnative tertiary structure.     

Differences in β-strand Populations of Monomeric Aβ40 and Aβ42

Using homonuclear ¹H NOESY spectra, with chemical shifts, ³J_H^N_H^α scalar couplings, residual dipolar couplings, and ¹H-¹⁵N NOEs, we have optimized and validated the conformational ensembles of the amyloid-β 1–40 (Aβ40) and amyloid-β 1–42 (Aβ42) peptides generated by molecular dynamics simulations. We find that both peptides have a diverse set of secondary structure elements including turns, helices, and antiparallel and parallel β-strands. The most significant difference in the structural ensembles of the two peptides is the type of β-hairpins and β-strands they populate. We find that Aβ42 forms a major antiparallel β-hairpin involving the central hydrophobic cluster residues (16–21) with residues 29–36, compatible with known amyloid fibril forming regions, whereas Aβ40 forms an alternative but less populated antiparallel β-hairpin between the central hydrophobic cluster and residues 9–13, that sometimes forms a β-sheet by association with residues 35–37. Furthermore, we show that the two additional C-terminal residues of Aβ42, in particular Ile-41, directly control the differences in the β-strand content found between the Aβ40 and Aβ42 structural ensembles. Integrating the experimental and theoretical evidence accumulated over the last decade, it is now possible to present monomeric structural ensembles of Aβ40 and Aβ42 consistent with available information that produce a plausible molecular basis for why Aβ42 exhibits greater fibrillization rates than Aβ40.     

Probing the helical forms of Ca2+‐guluronan junction zones in alginate gels by molecular dynamics 1: Duplexes

A molecular dynamics investigation of the helical forms adopted by (1→4)‐α‐L ‐guluronan in explicit water environment was carried out. Single chains and duplexes were modeled at 300 K starting both from 2₁ or 3₂ helical conformations and in the presence of a neutralizing amount of Ca²⁺ ions. All systems were allowed full conformational freedom. The initial perfect helices with integral screw symmetries were lost at the very beginning of simulations and two distinct behaviors were observed: At equilibrium the 2₁ models mostly retained the 2₁ local helical conformations while exploring the 3₂ ones the rest of the time. In duplexes the two chains, which behaved similarly, were well extended and slightly twisted. By contrast, the chains in 3₂ duplex models were dissimilar and explored a much broader conformational space in which 2₁ and 3₂ local helical conformations were dominant and equally represented but the 3₁ and other conformations were also present. The wide variety of conformations revealed in this study is consistent with the general difficulty in obtaining crystals of Ca²⁺‐guluronate with suitable lateral dimensions for crystallographic studies. © 2013 Wiley Periodicals, Inc. Biopolymers 99: 562–571, 2013.     

Distance-Based Metrics for Comparing Conformational Ensembles of Intrinsically Disordered Proteins

Intrinsically disordered proteins are proteins whose native functional states represent ensembles of highly diverse conformations. Such ensembles are a challenge for quantitative structure comparisons because their conformational diversity precludes optimal superimposition of the atomic coordinates necessary for deriving common similarity measures such as the root mean-square deviation of these coordinates. Here, we introduce superimposition-free metrics that are based on computing matrices of the Cα-Cα distance distributions within ensembles and comparing these matrices between ensembles. Differences between two matrices yield information on the similarity between specific regions of the polypeptide, whereas the global structural similarity is captured by the root mean-square difference between the medians of the Cα-Cα distance distributions of two ensembles. Together, our metrics enable rigorous investigations of structure-function relationships in conformational ensembles of intrinsically disordered proteins derived using experimental restraints or by molecular simulations and for proteins containing both structured and disordered regions.     

Uncovering Specific Electrostatic Interactions in the Denatured States of Proteins

The stability and folding of proteins are modulated by energetically significant interactions in the denatured state that is in equilibrium with the native state. These interactions remain largely invisible to current experimental techniques, however, due to the sparse population and conformational heterogeneity of the denatured-state ensemble under folding conditions. Molecular dynamics simulations using physics-based force fields can in principle offer atomistic details of the denatured state. However, practical applications are plagued with the lack of rigorous means to validate microscopic information and deficiencies in force fields and solvent models. This study presents a method based on coupled titration and molecular dynamics sampling of the denatured state starting from the extended sequence under native conditions. The resulting denatured-state pK_as allow for the prediction of experimental observables such as pH- and mutation-induced stability changes. I show the capability and use of the method by investigating the electrostatic interactions in the denatured states of wild-type and K12M mutant of NTL9 protein. This study shows that the major errors in electrostatics can be identified by validating the titration properties of the fragment peptides derived from the sequence of the intact protein. Consistent with experimental evidence, our simulations show a significantly depressed pK_a for Asp⁸ in the denatured state of wild-type, which is due to a nonnative interaction between Asp⁸ and Lys¹². Interestingly, the simulation also shows a nonnative interaction between Asp⁸ and Glu⁴⁸ in the denatured state of the mutant. I believe the presented method is general and can be applied to extract and validate microscopic electrostatics of the entire folding energy landscape.