Convergence of replica exchange molecular dynamics

Structural characterization of an intrinsically disordered protein complex using integrated small-angle neutron scattering and computing

Characterizing structural ensembles of intrinsically disordered proteins (IDPs) and intrinsically disordered regions (IDRs) of proteins is essential for studying structure-function relationships. Due to the different neutron scattering lengths of hydrogen and deuterium, selective labeling and contrast matching in small-angle neutron scattering (SANS) becomes an effective tool to study dynamic structures of disordered systems. However, experimental timescales typically capture measurements averaged over multiple conformations, leaving complex SANS data for disentanglement. We hereby demonstrate an integrated method to elucidate the structural ensemble of a complex formed by two IDRs. We use data from both full contrast and contrast matching with residue-specific deuterium labeling SANS experiments, microsecond all-atom molecular dynamics (MD) simulations with four molecular mechanics force fields, and an autoencoder-based deep learning (DL) algorithm. From our combined approach, we show that selective deuteration provides additional information that helps characterize structural ensembles. We find that among the four force fields, a99SB-disp and CHARMM36m show the strongest agreement with SANS and NMR experiments. In addition, our DL algorithm not only complements conventional structural analysis methods but also successfully differentiates NMR and MD structures which are indistinguishable on the free energy surface. Lastly, we present an ensemble that describes experimental SANS and NMR data better than MD ensembles generated by one single force field and reveal three clusters of distinct conformations. Our results demonstrate a new integrated approach for characterizing structural ensembles of IDPs.

Atomic visualization of flipped-back conformations of high mannose glycans interacting with cargo lectins: An MD simulation perspective

Protein-carbohydrate interactions play a crucial role in mediating several biomolecular recognition events. We attempt to unravel its intricacies by understanding how carbohydrate-binding proteins interpret the glycan code. We aim to decipher lectin-mediated recognition in the endoplasmic reticulum (ER), which plays a crucial role in ER-mediated quality control (ER-QC). The ER-QC functions in three phases-protein folding, transport, and degradation. Altered protein QC leads to ER-related storage disorders. Cargo transport proteins-Ergic53 and Vip36-necessary for maintaining cellular homeostasis-are our primary focus. They recognize monoglucosylated/high mannose N-glycans on the folded glycoproteins. This article reports on the first dynamic investigation of the ER cargo lectins in complex with the high mannose glycans using an advanced sampling technique-replica exchange molecular dynamics to decipher the inherent conformational heterogeneity and the binding mechanism. The study involves simulations for the proteins complexed with three high mannose glycans-Man8B, Man9, and mono-glucosylated glycan. The recognition process is captured using MD simulations to achieve mechanistic insights and characterize the dynamics of glycans in their native and bound states via dihedral angle analysis. Results indicate that the flipped conformation of the glycans was crucial in differentiating their interaction with the proteins. Similar conformers of the glycans are preferred for Ergic53 and Vip36 in their glycan recognition events. Ergic53 preferred Man8B while it was Man9 for Vip36, in coherence with the previous experimental reports. These simulations provide a computational microscopic purview of the mechanism at both spatial and temporal scales. The results correlate with the published experimental data on the specificities of these lectins.

The role of conformational change and key glutamic acid residues in the ClC-ec1 antiporter

The triple glutamine (Q) mutant (QQQ) structure of a Cl-/H+ antiporter from Escherichia coli (ClC-ec1) displaying a novel backbone arrangement has been used to challenge the long-held notion that Cl-/H+ antiporters do not operate through large conformational motions. The QQQ mutant substitutes the glutamine residue for an external glutamate E148, an internal glutamate E203, and a third glutamate E113 that hydrogen-bonds with E203. However, it is unknown if QQQ represents a physiologically relevant state, as well as how the protonation of the wild-type glutamates relates to the global dynamics. We herein apply continuous constant-pH molecular dynamics to investigate the H+-coupled dynamics of ClC-ec1. Although any large-scale conformational rearrangement upon acidification would be due to the accumulation of excess charge within the protein, protonation of the glutamates significantly impacts mainly the local structure and dynamics. Despite the fact that the extracellular pore enlarges at acidic pHs, an occluded ClC-ec1 within the active pH range of 3.5-7.5 requires a protonated E148 to facilitate extracellular Cl- release. E203 is also involved in the intracellular H+ transfer as an H+ acceptor. The water wire connection of E148 with the intracellular solution is regulated by the charge states of the E113/E203 dyad with coupled proton titration. However, the dynamics extracted from our simulations are not QQQ-like, indicating that the QQQ mutant does not represent the behavior of the wild-type ClC-ec1. These findings reinforce the necessity of having a protonatable residue at the E203 position in ClC-ec1 and suggest that a higher level of complexity exists for the intracellular H+ transfer in Cl-/H+ antiporters.

Myofilament-associated proteins with intrinsic disorder (MAPIDs) and their resolution by computational modeling

The cardiac sarcomere is a cellular structure in the heart that enables muscle cells to contract. Dozens of proteins belong to the cardiac sarcomere, which work in tandem to generate force and adapt to demands on cardiac output. Intriguingly, the majority of these proteins have significant intrinsic disorder that contributes to their functions, yet the biophysics of these intrinsically disordered regions (IDRs) have been characterized in limited detail. In this review, we first enumerate these myofilament-associated proteins with intrinsic disorder (MAPIDs) and recent biophysical studies to characterize their IDRs. We secondly summarize the biophysics governing IDR properties and the state-of-the-art in computational tools toward MAPID identification and characterization of their conformation ensembles. We conclude with an overview of future computational approaches toward broadening the understanding of intrinsic disorder in the cardiac sarcomere.

Gaussian-Accelerated Molecular Dynamics with the Weighted Ensemble Method: A Hybrid Method Improves Thermodynamic and Kinetic Sampling

Gaussian-accelerated molecular dynamics (GaMD) is a well-established enhanced sampling method for molecular dynamics simulations that effectively samples the potential energy landscape of the system by adding a boost potential, which smoothens the surface and lowers the energy barriers between states. GaMD is unable to give time-dependent properties such as kinetics directly. On the other hand, the weighted ensemble (WE) method can efficiently sample transitions between states with its many weighted trajectories, which directly yield rates and pathways. However, convergence to equilibrium conditions remains a challenge for the WE method. Hence, we have developed a hybrid method that combines the two methods, wherein GaMD is first used to sample the potential energy landscape of the system and WE is subsequently used to further sample the potential energy landscape and kinetic properties of interest. We show that the hybrid method can sample both thermodynamic and kinetic properties more accurately and quickly compared to using either method alone.

Accelerating atomistic simulations of proteins using multiscale enhanced sampling with independent tempering

Efficient sampling of the conformational space is essential for quantitative simulations of proteins. The multiscale enhanced sampling (MSES) method accelerates atomistic sampling by coupling it to a coarse-grained (CG) simulation. Bias from coupling to the CG model is removed using Hamiltonian replica exchange, such that one could benefit simultaneously from the high accuracy of atomistic models and fast dynamics of CG ones. Here, we extend MSES to allow independent control of the effective temperatures of atomistic and CG simulations, by directly scaling the atomistic and CG Hamiltonians. The new algorithm, named MSES with independent tempering (MSES-IT), supports more sophisticated Hamiltonian and temperature replica exchange protocols to further improve the sampling efficiency. Using a small but nontrivial β-hairpin, we show that setting the effective temperature of CG model in all conditions to its melting temperature maximizes structural transition rates at the CG level and promotes more efficient replica exchange and diffusion in the condition space. As the result, MSES-IT drive faster reversible transitions at the atomic level and leads to significant improvement in generating converged conformational ensembles compared to the original MSES scheme.

Statistical Mechanical Approximations to More Efficiently Determine Polymorph Free Energy Differences for Small Organic Molecules

Methods to efficiently determine the relative stability of polymorphs of organic crystals are highly desired in crystal structure predictions (CSPs). Current methodologies include calculating the free energy of static lattice phonons, quasi-harmonic approximations (QHA), and computing the full thermodynamic cycle using replica exchange molecular dynamics (REMD). We found that 13 out of the 29 systems minimized from experimental crystal structures restructured to a lower energy minimum when heated and annealed using REMD, a phenomenon that QHA alone cannot capture. Here, we present a series of methods that are intermediate in accuracy and expense between QHA and computing the full thermodynamic cycle, which can save 42-80% of the computational cost and introduces, on this benchmark, a relatively small (0.16 ± 0.04 kcal/mol) error relative to the full thermodynamic cycle. In particular, a method that Boltzmann weights harmonic free energies from along the trajectory of REMD replica exchange appears to be an appropriate intermediate between QHA and the full thermodynamic cycle using MD when screening crystal polymorph stability.

Advances in Molecular Dynamics Simulations and Enhanced Sampling Methods for the Study of Protein Systems

Molecular dynamics (MD) simulation is a rigorous theoretical tool that when used efficiently could provide reliable answers to questions pertaining to the structure-function relationship of proteins. Data collated from protein dynamics can be translated into useful statistics that can be exploited to sieve thermodynamics and kinetics crucial for the elucidation of mechanisms responsible for the modulation of biological processes such as protein-ligand binding and protein-protein association. Continuous modernization of simulation tools enables accurate prediction and characterization of the aforementioned mechanisms and these qualities are highly beneficial for the expedition of drug development when effectively applied to structure-based drug design (SBDD). In this review, current all-atom MD simulation methods, with focus on enhanced sampling techniques, utilized to examine protein structure, dynamics, and functions are discussed. This review will pivot around computer calculations of protein-ligand and protein-protein systems with applications to SBDD. In addition, we will also be highlighting limitations faced by current simulation tools as well as the improvements that have been made to ameliorate their efficiency.

Calmodulin-Calcineurin Interaction beyond the Calmodulin-Binding Region Contributes to Calcineurin Activation

Calcineurin (CaN) is a calcium-dependent phosphatase involved in numerous signaling pathways. Its activation is in part driven by the binding of calmodulin (CaM) to a CaM recognition region (CaMBR) within CaN's regulatory domain (RD). However, secondary interactions between CaM and the CaN RD may be necessary to fully activate CaN. Specifically, it is established that the CaN RD folds upon CaM binding and a region C-terminal to CaMBR, the "distal helix", assumes an α-helix fold and contributes to activation [Dunlap, T. B., et al. (2013) Biochemistry 52, 8643-8651]. We hypothesized in that previous study that this distal helix can bind CaM in a region distinct from the canonical CaMBR. To test this hypothesis, we utilized molecular simulations, including replica-exchange molecular dynamics, protein-protein docking, and computational mutagenesis, to determine potential distal helix-binding sites on CaM's surface. We isolated a potential binding site on CaM (site D) that facilitates moderate-affinity interprotein interactions and predicted that mutation of site D residues K30 and G40 on CaM would weaken CaN distal helix binding. We experimentally confirmed that two variants (K30E and G40D) indicate weaker binding of a phosphate substrate p-nitrophenyl phosphate to the CaN catalytic site by a phosphatase assay. This weakened substrate affinity is consistent with competitive binding of the CaN autoinhibition domain to the catalytic site, which we suggest is due to the weakened distal helix-CaM interactions. This study therefore suggests a novel mechanism for CaM regulation of CaN that may extend to other CaM targets.

Massive-Scale Binding Free Energy Simulations of HIV Integrase Complexes Using Asynchronous Replica Exchange Framework Implemented on the IBM WCG Distributed Network

To perform massive-scale replica exchange molecular dynamics (REMD) simulations for calculating binding free energies of protein-ligand complexes, we implemented the asynchronous replica exchange (AsyncRE) framework of the binding energy distribution analysis method (BEDAM) in implicit solvent on the IBM World Community Grid (WCG) and optimized the simulation parameters to reduce the overhead and improve the prediction power of the WCG AsyncRE simulations. We also performed the first massive-scale binding free energy calculations using the WCG distributed computing grid and 301 ligands from the SAMPL4 challenge for large-scale binding free energy predictions of HIV-1 integrase complexes. In total there are ∼10000 simulated complexes, ∼1 million replicas, and ∼2000 μs of aggregated MD simulations. Running AsyncRE MD simulations on the WCG requires accepting a trade-off between the number of replicas that can be run (breadth) and the number of full RE cycles that can be completed per replica (depth). As compared with synchronous Replica Exchange (SyncRE) running on tightly coupled clusters like XSEDE, on the WCG many more replicas can be launched simultaneously on heterogeneous distributed hardware, but each full RE cycle requires more overhead. We compared the WCG results with that from AutoDock and more advanced RE simulations including the use of flattening potentials to accelerate sampling of selected degrees of freedom of ligands and/or receptors related to slow dynamics due to high energy barriers. We propose a suitable strategy of RE simulations to refine high throughput docking results which can be matched to corresponding computing resources: from HPC clusters, to small or medium-size distributed campus grids, and finally to massive-scale computing networks including millions of CPUs like the resources available on the WCG.

Molecular Dynamics Simulations of Conformational Conversions in Transformer Proteins

A relatively recently discovered class of proteins known as transformer proteins undergo large-scale conformational conversions that change their supersecondary structure. These structural transformations lead to different configurations that perform different functions. We describe computational methods using molecular dynamics simulations that allow the determination of the specific amino acids that facilitate the conformational transformations. These investigations provide guidance on the location and type of amino acid mutations that can either enhance or inhibit the structural transitions that allow transformer proteins to perform multiple functions.

Combining multi-scale modelling methods to decipher molecular motions of a branching sucrase from glycoside-hydrolase family 70

Among α-transglucosylases from Glycoside-Hydrolase family 70, the ΔN₁₂₃-GB-CD2 enzyme derived from the bifunctional DSR-E from L. citreum NRRL B-1299 is particularly interesting as it was the first described engineered Branching Sucrase, not able to elongate glucan polymers from sucrose substrate. The previously reported overall structural organization of this multi-domain enzyme is an intricate U-shape fold conserved among GH70 enzymes which showed a certain conformational variability of the so-called domain V, assumed to play a role in the control of product structures, in available X-ray structures. Understanding the role of functional dynamics on enzyme reaction and substrate recognition is of utmost interest although it remains a challenge for biophysical methods. By combining long molecular dynamics simulation (1μs) and multiple analyses (NMA, PCA, Morelet Continuous Wavelet Transform and Cross Correlations Dynamics), we investigated here the dynamics of ΔN₁₂₃-GB-CD2 alone and in interaction with sucrose substrate. Overall, our results provide the detailed picture at atomic level of the hierarchy of motions occurring along different timescales and how they are correlated, in agreement with experimental structural data. In particular, detailed analysis of the different structural domains revealed cooperative dynamic behaviors such as twisting, bending and wobbling through anti- and correlated motions, and also two structural hinge regions, of which one was unreported. Several highly flexible loops surrounding the catalytic pocket were also highlighted, suggesting a potential role in the acceptor promiscuity of ΔN123-GBD-CD2. Normal modes and essential dynamics underlined an interesting two-fold dynamic of the catalytic domain A, pivoting about an axis splitting the catalytic gorge in two parts. The comparison of the conformational free energy landscapes using principal component analysis of the enzyme in absence or in presence of sucrose, also revealed a more harmonic basin when sucrose is bound with a shift population of the bending mode, consistent with the substrate binding event.

Mycolyltransferase from Mycobacterium tuberculosis in covalent complex with tetrahydrolipstatin provides insights into antigen 85 catalysis

Mycobacterium tuberculosis antigen 85 (Ag85) enzymes catalyze the transfer of mycolic acid (MA) from trehalose monomycolate to produce the mycolyl arabinogalactan (mAG) or trehalose dimycolate (TDM). These lipids define the protective mycomembrane of mycobacteria. The current model of substrate binding within the active sites of Ag85s for the production of TDM is not sterically and geometrically feasible; additionally, this model does not account for the production of mAG. Furthermore, this model does not address how Ag85s limit the hydrolysis of the acyl-enzyme intermediate while catalyzing acyl transfer. To inform an updated model, we obtained an Ag85 acyl-enzyme intermediate structure that resembles the mycolated form. Here, we present a 1.45-Å X-ray crystal structure of M. tuberculosis Ag85C covalently modified by tetrahydrolipstatin (THL), an esterase inhibitor that suppresses M. tuberculosis growth and mimics structural attributes of MAs. The mode of covalent inhibition differs from that observed in the reversible inhibition of the human fatty-acid synthase by THL. Similarities between the Ag85-THL structure and previously determined Ag85C structures suggest that the enzyme undergoes structural changes upon acylation, and positioning of the peptidyl arm of THL limits hydrolysis of the acyl-enzyme adduct. Molecular dynamics simulations of the modeled mycolated-enzyme form corroborate the structural analysis. From these findings, we propose an alternative arrangement of substrates that rectifies issues with the previous model and suggest a direct role for the β-hydroxy of MA in the second half-reaction of Ag85 catalysis. This information affords the visualization of a complete mycolyltransferase catalytic cycle.

Insights into the Molecular Mechanisms of Alzheimer's and Parkinson's Diseases with Molecular Simulations: Understanding the Roles of Artificial and Pathological Missense Mutations in Intrinsically Disordered Proteins Related to Pathology

Amyloid-β and α-synuclein are intrinsically disordered proteins (IDPs), which are at the center of Alzheimer's and Parkinson's disease pathologies, respectively. These IDPs are extremely flexible and do not adopt stable structures. Furthermore, both amyloid-β and α-synuclein can form toxic oligomers, amyloid fibrils and other type of aggregates in Alzheimer's and Parkinson's diseases. Experimentalists face challenges in investigating the structures and thermodynamic properties of these IDPs in their monomeric and oligomeric forms due to the rapid conformational changes, fast aggregation processes and strong solvent effects. Classical molecular dynamics simulations complement experiments and provide structural information at the atomic level with dynamics without facing the same experimental limitations. Artificial missense mutations are employed experimentally and computationally for providing insights into the structure-function relationships of amyloid-β and α-synuclein in relation to the pathologies of Alzheimer's and Parkinson's diseases. Furthermore, there are several natural genetic variations that play a role in the pathogenesis of familial cases of Alzheimer's and Parkinson's diseases, which are related to specific genetic defects inherited in dominant or recessive patterns. The present review summarizes the current understanding of monomeric and oligomeric forms of amyloid-β and α-synuclein, as well as the impacts of artificial and pathological missense mutations on the structural ensembles of these IDPs using molecular dynamics simulations. We also emphasize the recent investigations on residual secondary structure formation in dynamic conformational ensembles of amyloid-β and α-synuclein, such as β-structure linked to the oligomerization and fibrillation mechanisms related to the pathologies of Alzheimer's and Parkinson's diseases. This information represents an important foundation for the successful and efficient drug design studies.

Inference of Calmodulin's Ca2+-Dependent Free Energy Landscapes via Gaussian Mixture Model Validation

A free energy landscape estimation method based on the well-known Gaussian mixture model (GMM) is used to compare the efficiencies of thermally enhanced sampling methods with respect to regular molecular dynamics. The simulations are carried out on two binding states of calmodulin, and the free energy estimation method is compared with other estimators using a toy model. We show that GMM with cross-validation provides a robust estimate that is not subject to overfitting. The continuous nature of Gaussians provides better estimates on sparse data than canonical histogramming. We find that diffusion properties determine the sampling method effectiveness, such that diffusion-dominated apo calmodulin is most efficiently sampled by regular molecular dynamics, while holo calmodulin, with its rugged free energy landscape, is better sampled by enhanced sampling methods.

Statistical Measures to Quantify Similarity between Molecular Dynamics Simulation Trajectories

Molecular dynamics simulation is commonly employed to explore protein dynamics. Despite the disparate timescales between functional mechanisms and molecular dynamics (MD) trajectories, functional differences are often inferred from differences in conformational ensembles between two proteins in structure-function studies that investigate the effect of mutations. A common measure to quantify differences in dynamics is the root mean square fluctuation (RMSF) about the average position of residues defined by Cα-atoms. Using six MD trajectories describing three native/mutant pairs of beta-lactamase, we make comparisons with additional measures that include Jensen-Shannon, modifications of Kullback-Leibler divergence, and local p-values from 1-sample Kolmogorov-Smirnov tests. These additional measures require knowing a probability density function, which we estimate by using a nonparametric maximum entropy method that quantifies rare events well. The same measures are applied to distance fluctuations between Cα-atom pairs. Results from several implementations for quantitative comparison of a pair of MD trajectories are made based on fluctuations for on-residue and residue-residue local dynamics. We conclude that there is almost always a statistically significant difference between pairs of 100 ns all-atom simulations on moderate-sized proteins as evident from extraordinarily low p-values.

Consistent Principal Component Modes from Molecular Dynamics Simulations of Proteins

Principal component analysis is a technique widely used for studying the movements of proteins using data collected from molecular dynamics simulations. In spite of its extensive use, the technique has a serious drawback: equivalent simulations do not afford the same PC-modes. In this article, we show that concatenating equivalent trajectories and calculating the PC-modes from the concatenated one significantly enhances the reproducibility of the results. Moreover, the consistency of the modes can be systematically improved by adding more individual trajectories to the concatenated one.

Aggregation of Chameleon Peptides: Implications of α-Helicity in Fibril Formation

We investigate the relationship between the inherent secondary structure and aggregation propensity of peptides containing chameleon sequences (i.e., sequences that can adopt either α or β structure depending on context) using a combination of replica exchange molecular dynamics simulations, ion-mobility mass spectrometry, circular dichroism, and transmission electron microscopy. We focus on an eight-residue long chameleon sequence that can adopt an α-helical structure in the context of the iron-binding protein from Bacillus anthracis (PDB id 1JIG ) and a β-strand in the context of the baculovirus P35 protein (PDB id 1P35 ). We show that the isolated chameleon sequence is intrinsically disordered, interconverting between α-helical and β-rich conformations. The inherent conformational plasticity of the sequence can be constrained by addition of flanking residues with a given secondary structure propensity. Intriguingly, we show that the chameleon sequence with helical flanking residues aggregates rapidly into fibrils, whereas the chameleon sequence with flanking residues that favor β-conformations has weak aggregation propensity. This work sheds new insights into the possible role of α-helical intermediates in fibril formation.

Structural Study of Cell Attachment Peptide Derived from Laminin by Molecular Dynamics Simulation

Peptides with cell attachment activity are beneficial component of biomaterials for tissue engineering. Conformational structure is one of the important factors for the biological activities. The EF1 peptide (DYATLQLQEGRLHFMFDLG) derived from laminin promotes cell spreading and cell attachment activity mediated by α2β1 integrin. Although the sequence of the EF2 peptide (DFATVQLRNGFPYFSYDLG) is homologous sequence to that of EF1, EF2 does not promote cell attachment activity. To determine whether there are structural differences between EF1 and EF2, we performed replica exchange molecular dynamics (REMD) simulations and conventional molecular dynamics (MD) simulations. We found that EF1 and EF2 had β-sheet structure as a secondary structure around the global minimum. However, EF2 had variety of structures around the global minimum compared with EF1 and has easily escaped from the bottom of free energy. The structural fluctuation of the EF1 is smaller than that of the EF2. The structural variation of EF2 is related to these differences in the structural fluctuation and the number of the hydrogen bonds (H-bonds). From the analysis of H-bonds in the β-sheet, the number of H-bonds in EF1 is larger than that in EF2 in the time scale of the conventional MD simulation, suggesting that the formation of H-bonds is related to the differences in the structural fluctuation between EF1 and EF2. From the analysis of other non-covalent interactions in the amino acid sequences of EF1 and EF2, EF1 has three pairs of residues with hydrophobic interaction, and EF2 has two pairs. These results indicate that several non-covalent interactions are important for structural stabilization. Consequently, the structure of EF1 is stabilized by H-bonds and pairs of hydrophobic amino acids in the terminals. Hence, we propose that non-covalent interactions around N-terminal and C-terminal of the peptides are crucial for maintaining the β-sheet structure of the peptides.

Folding Simulations of an α-Helical Hairpin Motif αtα with Residue-Specific Force Fields

α-Helical hairpin (two-helix bundle) is a structure motif composed of two interacting helices connected by a turn or a short loop. It is an important model for protein folding studies, filling the gap between isolated α-helix and larger all-α domains. Here, we present, for the first time, successful folding simulations of an α-helical hairpin. Our RSFF1 and RSFF2 force fields give very similar predicted structures of this αtα peptide, which is in good agreement with its NMR structure. Our simulations also give site-specific stability of α-helix formation in good agreement with amide hydrogen exchange experiments. Combining the folding free energy landscapes and analyses of structures sampled in five different ranges of the fraction of native contacts (Q), a folding mechanism of αtα is proposed. The most stable sites of Q9-E15 in helix-1 and E24-A30 in helix-2 close to the loop region act as the folding initiation sites. The formation of interhelix side-chain contacts also initiates near the loop region, but some residues in the central parts of the two helices also form contacts quite early. The two termini fold at a final stage, and the loop region remains flexible during the whole folding process. This mechanism is similar to the "zipping out" pathway of β-hairpin folding.

Improving Convergence of Replica-Exchange Simulations through Coupling to a High-Temperature Structure Reservoir

Parallel tempering or replica-exchange molecular dynamics (REMD) significantly increases the efficiency of conformational sampling for complex molecular systems. However, obtaining converged data with REMD remains challenging, especially for large systems with complex topologies. We propose a new variant to REMD where the replicas are also permitted to exchange with an ensemble of structures that have been generated in advance using high-temperature MD simulations, similar in spirit to J-walking methods. We tested this approach on two non-trivial model systems, a β-hairpin and a 3-stranded β-sheet and compared the results to those obtained from very long (>100 ns) standard REMD simulations. The resulting ensembles were indistinguishable, including relative populations of different conformations on the unfolded state. The use of the reservoir is shown to significantly reduce the time required for convergence.

Simulated Solute Tempering

For the enhanced conformational sampling in molecular dynamics (MD) simulations, we present "simulated solute tempering" (SST) which is an easy to implement variant of simulated tempering. SST extends conventional simulated tempering (CST) by key concepts of "replica exchange with solute tempering" (REST, Liu et al. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 13749). We have applied SST, CST, and REST to molecular dynamics (MD) simulations of an alanine octapeptide in explicit water. The weight parameters required for CST and SST are determined by two different formulas whose performance is compared. For SST only one of them yields a uniform sampling of the temperature space. Compared to CST and REST, SST provides the highest exchange probabilities between neighboring rungs in the temperature ladder. Concomitantly, SST leads to the fastest diffusion of the simulation system through the temperature space, in particular, if the "even-odd" exchange scheme is employed in SST. As a result, SST exhibits the highest sampling speed of the investigated tempering methods.

Resolution Exchange Simulation with Incremental Coarsening

We previously developed an algorithm, called "resolution exchange", which improves canonical sampling of atomic resolution models by swapping conformations between high- and low-resolution simulations. Here, we demonstrate a generally applicable incremental coarsening procedure and apply the algorithm to a larger peptide, met-enkephalin. In addition, we demonstrate a combination of resolution and temperature exchange, in which the coarser simulations are also at elevated temperatures. Both simulations are implemented in a "top-down" mode, to allow efficient allocation of CPU time among the different replicas.

Large-scale asynchronous and distributed multidimensional replica exchange molecular simulations and efficiency analysis

We describe methods to perform replica exchange molecular dynamics (REMD) simulations asynchronously (ASyncRE). The methods are designed to facilitate large scale REMD simulations on grid computing networks consisting of heterogeneous and distributed computing environments as well as on homogeneous high-performance clusters. We have implemented these methods on NSF (National Science Foundation) XSEDE (Extreme Science and Engineering Discovery Environment) clusters and BOINC (Berkeley Open Infrastructure for Network Computing) distributed computing networks at Temple University and Brooklyn College at CUNY (the City University of New York). They are also being implemented on the IBM World Community Grid. To illustrate the methods, we have performed extensive (more than 60 ms in aggregate) simulations for the beta-cyclodextrin-heptanoate host-guest system in the context of one- and two-dimensional ASyncRE, and we used the results to estimate absolute binding free energies using the binding energy distribution analysis method. We propose ways to improve the efficiency of REMD simulations: these include increasing the number of exchanges attempted after a specified molecular dynamics (MD) period up to the fast exchange limit and/or adjusting the MD period to allow sufficient internal relaxation within each thermodynamic state. Although ASyncRE simulations generally require long MD periods (>picoseconds) per replica exchange cycle to minimize the overhead imposed by heterogeneous computing networks, we found that it is possible to reach an efficiency similar to conventional synchronous REMD, by optimizing the combination of the MD period and the number of exchanges attempted per cycle.

Improving the replica-exchange molecular-dynamics method for efficient sampling in the temperature space

Replica-exchange molecular dynamics (REMD) is a popular sampling method in the molecular simulation. By frequently exchanging the replicas at different temperatures, the molecule can jump out of the minima and sample efficiently in the conformational space. Although REMD has been shown to be practical in a lot of applications, it does have a critical limitation. All the replicas at all the temperatures must be simulated for a period between the replica-exchange steps. This may be problematic for the reaction with high free energy barriers. In that case, too many replicas are required in the simulation. To reduce the calculation quantity and improve its performance, in this paper we propose a modified REMD method. During the simulation, each replica at each temperature can stay in either the active or inactive state and only switch between the states at the exchange step. In the active state, the replica moves freely in the canonical ensemble by the normal molecular dynamics, and in the inactive state, the replica is frozen temporarily until the next exchange step. The number of the replicas in the active states (active replicas) depends on the number of CPUs in the computer. Using the additional inactive replicas, one can perform an REMD simulation in a wider temperature space. The practical applications show that the modified REMD method is reliable. With the same number of active replicas, this REMD method can produce a more reasonable free energy surface around the free energy minima than the standard REMD method.

Replica exchange with guided annealing for accelerated sampling of disordered protein conformations

We critically examine a recently proposed convective replica exchange (cRE) method for enhanced sampling of protein conformation based on theoretical and numerical analysis. The results demonstrate that cRE and related replica exchange with guided annealing (RE-GA) schemes lead to unbalanced exchange attempt probabilities and break detailed balance whenever the system undergoes slow conformational transitions (relative to the temperature diffusion timescale). Nonetheless, numerical simulations suggest that approximate canonical ensembles can be generated for systems with small conformational transition barriers. This suggests that RE-GA maybe suitable for simulating intrinsically disordered proteins, an important class of newly recognized functional proteins. The efficacy of RE-GA is demonstrated by calculating the conformational ensembles of intrinsically disordered kinase inducible domain protein. The results show that RE-GA helps the protein to escape nonspecific compact states more efficiently and provides several fold speedups in generating converged and largely correct ensembles compared to the standard temperature RE.

Surface Isotope Segregation as a Probe of Temperature in Water Nanoclusters

Using ring polymer molecular dynamics simulations, we examine equilibrium and dynamical characteristics of solid-like, aqueous clusters that combine isotopic mixtures of HDO dilute in H2O, at temperatures intermediate between 50 and 175 K. In particular, we focus attention on the relative thermodynamic stabilities of the two isotopes at dangling hydrogen bond sites. The water octamer is analyzed as a reference system. For this aggregate, decreasing temperature yields a gradual stabilization of the light isotope at dangling sites in molecules acting as single-donor-double-acceptors of hydrogen bonds. At T ∼ 50 K, the imbalance between the corresponding quantum kinetic energies leads to a free energy difference between dangling and hydrogen bonded sites of the order of ∼2kBT. Similar free energy differences were found at dangling sites in Nw = 50 water clusters. The extent of the H/D segregation can be adequately monitored by modifications in the peak intensity of the high frequency shoulder of the stretching band of the infrared spectrum. These signals, in turn, represent a potential experimental signature of the elusive temperature of clusters in molecular beams.

Recognition of the HIV capsid by the TRIM5α restriction factor is mediated by a subset of pre-existing conformations of the TRIM5α SPRY domain

The binding of the TRIM5α restriction factor to the HIV capsid is mediated by the C-terminal SPRY domain of TRIM5α. Atomic-level details of this host-pathogen interaction, which involves mobile variable loops of the SPRY domain, remain unclear. Some of the key determinants of restriction are encompassed by the long and disordered v1 loop of the SPRY domain. We applied molecular modeling to elucidate the conformational repertoire of the v1 loop and its role in the interaction with the capsid. All-atom replica exchange molecular dynamics revealed multiple transient, interconverting states of the v1 loop consistent with the intrinsic disorder observed experimentally. The docking of the SPRY conformations representing 10 most populated states onto the high-resolution model of the assembled HIV-1 capsid revealed that a subset of v1 conformations produced plausible binding poses, in which the SPRY domain binds close to the pseudo-2-fold symmetry axis and the v1 loop spans the interhexamer gap. Such binding mode is well supported by the NMR binding data and known escape mutants. We speculate that the binding mode that involves interaction of the capsid with a subset of preexisting SPRY conformations arising from the intrinsic disorder of the v1 loop may explain the remarkable ability of TRIM5α to resist viral evasion by mutagenesis and to restrict divergent retroviruses.

Multidimensional Replica Exchange Molecular Dynamics Yields a Converged Ensemble of an RNA Tetranucleotide

A necessary step to properly assess and validate the performance of force fields for biomolecules is to exhaustively sample the accessible conformational space, which is challenging for large RNA structures. Given questions regarding the reliability of modeling RNA structure and dynamics with current methods, we have begun to use RNA tetranucleotides to evaluate force fields. These systems, though small, display considerable conformational variability and complete sampling with standard simulation methods remains challenging. Here we compare and discuss the performance of known variations of replica exchange molecular dynamics (REMD) methods, specifically temperature REMD (T-REMD), Hamiltonian REMD (H-REMD), and multidimensional REMD (M-REMD) methods, which have been implemented in Amber’s accelerated GPU code. Using two independent simulations, we show that M-REMD not only makes very efficient use of emerging large-scale GPU clusters, like Blue Waters at the University of Illinois, but also is critically important in generating the converged ensemble more efficiently than either T-REMD or H-REMD. With 57.6 μs aggregate sampling of a conformational ensemble with M-REMD methods, the populations can be compared to NMR data to evaluate force field reliability and further understand how putative changes to the force field may alter populations to be in more consistent agreement with experiment.

Structures and free energy landscapes of the A53T mutant-type α-synuclein protein and impact of A53T mutation on the structures of the wild-type α-synuclein protein with dynamics

The A53T genetic missense mutation of the wild-type α-synuclein (αS) protein was initially identified in Greek and Italian families with familial Parkinson's disease. Detailed understanding of the structures and the changes induced in the wild-type αS structure by the A53T mutation, as well as establishing the direct relationships between the rapid conformational changes and free energy landscapes of these intrinsically disordered fibrillogenic proteins, helps to enhance our fundamental knowledge and to gain insights into the pathogenic mechanism of Parkinson's disease. We employed extensive parallel tempering molecular dynamics simulations along with thermodynamic calculations to determine the secondary and tertiary structural properties as well as the conformational free energy surfaces of the wild-type and A53T mutant-type αS proteins in an aqueous solution medium using both implicit and explicit water models. The confined aqueous volume effect in the simulations of disordered proteins using an explicit model for water is addressed for a model disordered protein. We also assessed the stabilities of the residual secondary structure component interconversions in αS based on free energy calculations at the atomic level with dynamics using our recently developed theoretical strategy. To the best of our knowledge, this study presents the first detailed comparison of the structural properties linked directly to the conformational free energy landscapes of the monomeric wild-type and A53T mutant-type α-synuclein proteins in an aqueous solution environment. Results demonstrate that the β-sheet structure is significantly more altered than the helical structure upon A53T mutation of the monomeric wild-type αS protein in aqueous solution. The β-sheet content close to the mutation site in the N-terminal region is more abundant while the non-amyloid-β component (NAC) and C-terminal regions show a decrease in β-sheet abundance upon A53T mutation. Obtained results utilizing our new theoretical strategy show that the residual secondary structure conversion stabilities resulting in α-helix formation are not significantly affected by the mutation. Interestingly, the residual secondary structure conversion stabilities show that secondary structure conversions resulting in β-sheet formation are influenced by the A53T mutation and the most stable residual transition yielding β-sheet occurs directly from the coil structure. Long-range interactions detected between the NAC region and the N- or C-terminal regions of the wild-type αS disappear upon A53T mutation. The A53T mutant-type αS structures are thermodynamically more stable than those of the wild-type αS protein structures in aqueous solution. Overall, the higher propensity of the A53T mutant-type αS protein to aggregate in comparison to the wild-type αS protein is related to the increased β-sheet formation and lack of strong intramolecular long-range interactions in the N-terminal region in comparison to its wild-type form. The specific residual secondary structure component stabilities reported herein provide information helpful for designing and synthesizing small organic molecules that can block the β-sheet forming residues, which are reactive toward aggregation.

Structures of the E46K mutant-type α-synuclein protein and impact of E46K mutation on the structures of the wild-type α-synuclein protein

The E46K genetic missense mutation of the wild-type α-synuclein protein was recently identified in a family of Spanish origin with hereditary Parkinson's disease. Detailed understanding of the structures of the monomeric E46K mutant-type α-synuclein protein as well as the impact of the E46K missense mutation on the conformations and free energy landscapes of the wild-type α-synuclein are required for gaining insights into the pathogenic mechanism of Parkinson's disease. In this study, we use extensive parallel tempering molecular dynamics simulations along with thermodynamic calculations to assess the secondary and tertiary structural properties as well as the conformational preferences of the monomeric wild-type and E46K mutant-type α-synuclein proteins in an aqueous solution environment. We also present the residual secondary structure component conversion stabilities with dynamics using a theoretical strategy, which we most recently developed. To the best of our knowledge, this study presents the first detailed comparison of the structural and thermodynamic properties of the wild-type and E46K mutant-type α-synuclein proteins in an aqueous solution environment at the atomic level with dynamics. We find that the E46K mutation results not only in local but also in long-range changes in the structural properties of the wild-type α-synuclein protein. The mutation site shows a significant decrease in helical content as well as a large increase in β-sheet structure formation upon E46K mutation. In addition, the β-sheet content of the C-terminal region increases significantly in the E46K mutant-type αS in comparison to the wild-type αS. Our theoretical strategy developed to assess the thermodynamic preference of secondary structure transitions indicates that this shift in secondary structure is the result of a decrease in the thermodynamic preference of turn to helix conversions while the coil to β-sheet preference increases for these residues. Long-range intramolecular protein interactions of the C-terminal with the N-terminal and NAC regions increase upon E46K mutation, resulting in more compact structures for the E46K mutant-type rather than wild-type αS. However, the E46K mutant-type αS structures are less stable than the wild-type αS. Overall, our results show that the E46K mutant-type αS has a higher propensity to aggregate than the wild-type αS and that the N-terminal and C-terminal regions are reactive toward fibrillization and aggregation upon E46K mutation and we explain the associated reasons based on the structural properties herein. Small molecules or drugs that can block the specific residues forming abundant β-sheet structure, which we report here, might help to reduce the reactivity of these intrinsically disordered fibrillogenic proteins toward aggregation and their toxicity.

Structures and free energy landscapes of the wild-type and A30P mutant-type α-synuclein proteins with dynamics

The genetic missense A30P mutation of the wild-type α-synuclein protein results in the replacement of the 30th amino acid residue from alanine (Ala) to proline (Pro) and was initially found in the members of a German family who developed Parkinson's disease. Even though the structures of these proteins have been measured before, detailed understanding about the structures and their relationships with free energy landscapes is lacking, which is of interest to provide insights into the pathogenic mechanism of Parkinson's disease. We report the secondary and tertiary structures and conformational free energy landscapes of the wild-type and A30P mutant-type α-synuclein proteins in an aqueous solution environment via extensive parallel tempering molecular dynamics simulations along with thermodynamic calculations. In addition, we present the residual secondary structure component transition stabilities at the atomic level with dynamics in terms of free energy change calculations using a new strategy that we reported most recently. Our studies yield new interesting results; for instance, we find that the A30P mutation has local as well as long-range effects on the structural properties of the wild-type α-synuclein protein. The helical content at Ala18-Gly31 is less prominent in comparison to the wild-type α-synuclein protein. The β-sheet structure abundance decreases in the N-terminal region upon A30P mutation of the wild-type α-synuclein, whereas the NAC and C-terminal regions possess larger tendencies for β-sheet structure formation. Long-range intramolecular protein interactions are less abundant upon A30P mutation, especially between the NAC and C-terminal regions, which is linked to the less compact and less stable structures of the A30P mutant-type rather than the wild-type α-synuclein protein. Results including the usage of our new strategy for secondary structure transition stabilities show that the A30P mutant-type α-synuclein tendency toward aggregation is higher than the wild-type α-synuclein but we also find that the C-terminal and NAC regions of the A30P mutant-type α-synuclein are reactive toward fibrillzation and aggregation based on atomic level studies with dynamics in an aqueous solution environment. Therefore, we propose that small molecules or drugs blocking the specific residues, which we report herein, located in the NAC- and C-terminal regions of the A30P mutant-type α-synuclein protein might help to reduce the toxicity of the A30P mutant-type α-synuclein protein.

Replica exchanging self-guided Langevin dynamics for efficient and accurate conformational sampling

This work presents a replica exchanging self-guided Langevin dynamics (RXSGLD) simulation method for efficient conformational searching and sampling. Unlike temperature-based replica exchanging simulations, which use high temperatures to accelerate conformational motion, this method uses self-guided Langevin dynamics (SGLD) to enhance conformational searching without the need to elevate temperatures. A RXSGLD simulation includes a series of SGLD simulations, with simulation conditions differing in the guiding effect and/or temperature. These simulation conditions are called stages and the base stage is one with no guiding effect. Replicas of a simulation system are simulated at the stages and are exchanged according to the replica exchanging probability derived from the SGLD partition function. Because SGLD causes less perturbation on conformational distribution than high temperatures, exchanges between SGLD stages have much higher probabilities than those between different temperatures. Therefore, RXSGLD simulations have higher conformational searching ability than temperature based replica exchange simulations. Through three example systems, we demonstrate that RXSGLD can generate target canonical ensemble distribution at the base stage and achieve accelerated conformational searching. Especially for large systems, RXSGLD has remarkable advantages in terms of replica exchange efficiency, conformational searching ability, and system size extensiveness.

Comparison between empirical protein force fields for the simulation of the adsorption behavior of structured LK peptides on functionalized surfaces

All-atom empirical molecular mechanics protein force fields, which have been developed to represent the energetics of peptide folding behavior in aqueous solution, have not been parameterized for protein interactions with solid material surfaces. As a result, their applicability for representing the adsorption behavior of proteins with functionalized material surfaces should not be assumed. To address this issue, we conducted replica-exchange molecular dynamics simulations of the adsorption behavior of structured peptides to functionalized surfaces using three protein force fields that are widely used for the simulation of peptide adsorption behavior: CHARMM22, AMBER94, and OPLS-AA. Simulation results for peptide structure both in solution and when adsorbed to the surfaces were compared to experimental results for similar peptide-surface systems to provide a means of evaluating and comparing the performance of these three force fields for this type of application. Substantial differences in both solution and adsorbed peptide conformations were found amongst these three force fields, with the CHARMM22 force field found to most closely match experimental results.

Effects of solvents on the intrinsic propensity of peptide backbone conformations

We investigated the effects of solvents on the intrinsic propensity of peptide backbone conformations based on molecular dynamics simulations. The results show that compared with pure water, aqueous urea decreases the helix propensity. In comparison, methanol decreases the polyproline II (PPII) propensity. Such a solvent dependence of the intrinsic propensity of the backbone conformation is correlated with the solvent dependence of the hydration of the backbone groups and the formation probability of the local intrapeptide hydrogen bonds. Aqueous urea which has low ability to stabilize the local intrapeptide hydrogen bonds disfavors the helical conformation. Whereas, methanol which has low ability to hydrate the backbone groups disfavors the polyproline II conformation. In addition, the solvent effects can be further modulated by the side chains of the peptides. The solvent effects of the intrinsic propensity of peptide backbone conformations observed in this work suggest that changing the intrinsic propensity of the protein backbone conformations can partly contribute to the solvent-induced protein structure and dynamics variations. These results will be useful in understanding the solvent dependence of the conformational distributions of the unfolded proteins or peptides (or intrinsically disordered proteins) in which the global tertiary interactions are less important than that in the well-folded proteins.

Intrinsically disordered proteins in a physics-based world

Intrinsically disordered proteins (IDPs) are a newly recognized class of functional proteins that rely on a lack of stable structure for function. They are highly prevalent in biology, play fundamental roles, and are extensively involved in human diseases. For signaling and regulation, IDPs often fold into stable structures upon binding to specific targets. The mechanisms of these coupled binding and folding processes are of significant importance because they underlie the organization of regulatory networks that dictate various aspects of cellular decision-making. This review first discusses the challenge in detailed experimental characterization of these heterogeneous and dynamics proteins and the unique and exciting opportunity for physics-based modeling to make crucial contributions, and then summarizes key lessons from recent de novo simulations of the structure and interactions of several regulatory IDPs.

Nitrile bonds as infrared probes of electrostatics in ribonuclease S

Three different nitrile-containing amino acids, p-cyanophenylalanine, m-cyanophenylalanine, and S-cyanohomocysteine, have been introduced near the active site of the semisynthetic enzyme ribonuclease S (RNase S) to serve as probes of electrostatic fields. Vibrational Stark spectra, measured directly on the probe-modified proteins, confirm the predominance of the linear Stark tuning rate in describing the sensitivity of the nitrile stretch to external electric fields, a necessary property for interpreting observed frequency shifts as a quantitative measure of local electric fields that can be compared with simulations. The X-ray structures of these nitrile-modified RNase variants and enzymatic assays demonstrate minimal perturbation to the structure and function, respectively, by the probes and provide a context for understanding the influence of the environment on the nitrile stretching frequency. We examine the ability of simulation techniques to recapitulate the spectroscopic properties of these nitriles as a means to directly test a computational electrostatic model for proteins, specifically that in the ubiquitous Amber-99 force field. Although qualitative agreement between theory and experiment is observed for the largest shifts, substantial discrepancies are observed in some cases, highlighting the ongoing need for experimental metrics to inform the development of theoretical models of electrostatic fields in proteins.