Force field bias in protein folding simulations

In silico vaccine design for Yersinia enterocolitica: A comprehensive approach to enhanced immunogenicity, efficacy and protection

Yersinia enterocolitica, a foodborne pathogen, has emerged as a significant public health concern due to its increased prevalence and multidrug resistance. This study employed reverse vaccinology to identify novel vaccine candidates against Y. enterocolitica through comprehensive in silico analyses. The core genome's conserved protein translocase subunit SecY was selected as the target, and potential B-cell, MHC class I, and MHC class II epitopes were mapped. 3B-cell epitopes, 3 MHCI and 11 MHCII epitopes were acquired. A multi-epitope vaccine construct was designed by incorporating the identified epitopes, TLR4 Agonist was used as adjuvants to enhance the immunogenic response. EAAAK, CPGPG and AYY linkers were used to form a vaccine construct, followed by extensive computational evaluations. The vaccine exhibited desirable physicochemical properties, stable secondary and tertiary structures as evaluated by PDBSum and trRosetta. Moreover, favorable interactions with the human Toll-like receptor 4 (TLR4) was observed by ClusPro. Population coverage analysis estimated the vaccine's applicability across 99.74 % in diverse populations. In addition, molecular dynamics simulations and normal mode analysis confirmed the vaccine's structural stability and dynamics in a simulated biological environment. Furthermore, codon optimization and in silico cloning facilitated the evaluation of the vaccine's expression potential in E. coli and pET-28a was used a recombinant plasmid. This study provides a promising foundation for the development of an efficacious vaccine against Y. enterocolitica infections.

DREAMweb: An online tool for graph-based modeling of NMR protein structure

The value of accurate protein structural models closely conforming to the experimental data is indisputable. DREAMweb deploys an improved DREAM algorithm, DREAMv2, that incorporates a tighter bound in the constraint set of the underlying optimization approach. This reduces the artifacts while modeling the protein structure by solving the distance-geometry problem. DREAMv2 follows a bottom-up strategy of building smaller substructures for regions with a larger concentration of experimental bounds and consolidating them before modeling the rest of the protein structure. This improves secondary structure conformance in the final models consistent with experimental data. The proposed method efficiently models regions with sparse coverage of experimental data by reducing the possibility of artifacts compared to DREAM. To balance performance and accuracy, smaller substructures (∼ 200 $\sim 200$ atoms) are solved in this regime, allowing faster builds for the other parts under relaxed conditions. DREAMweb is accessible as an internet resource. The improvements in results are showcased through benchmarks on 10 structures. DREAMv2 can be used in tandem with any NMR-based protein structure determination workflow, including an iterative framework where the NMR assignment for the NOESY spectra is incomplete or ambiguous. DREAMweb is freely available for public use at http://pallab.cds.iisc.ac.in/DREAM/ and downloadable at https://github.com/niladriranjandas/DREAMv2.git.

Simulation-driven design of stabilized SARS-CoV-2 spike S2 immunogens

The full-length prefusion-stabilized SARS-CoV-2 spike (S) is the principal antigen of COVID-19 vaccines. Vaccine efficacy has been impacted by emerging variants of concern that accumulate most of the sequence modifications in the immunodominant S1 subunit. S2, in contrast, is the most evolutionarily conserved region of the spike and can elicit broadly neutralizing and protective antibodies. Yet, S2's usage as an alternative vaccine strategy is hampered by its general instability. Here, we use a simulation-driven approach to design S2-only immunogens stabilized in a closed prefusion conformation. Molecular simulations provide a mechanistic characterization of the S2 trimer's opening, informing the design of tryptophan substitutions that impart kinetic and thermodynamic stabilization. Structural characterization via cryo-EM shows the molecular basis of S2 stabilization in the closed prefusion conformation. Informed by molecular simulations and corroborated by experiments, we report an engineered S2 immunogen that exhibits increased protein expression, superior thermostability, and preserved immunogenicity against sarbecoviruses.

100 Years of the Lennard-Jones Potential

It is now 100 years since Lennard-Jones published his first paper introducing the now famous potential that bears his name. It is therefore timely to reflect on the many achievements, as well as the limitations, of this potential in the theory of atomic and molecular interactions, where applications range from descriptions of intermolecular forces to molecules, clusters, and condensed matter.

Computational Approaches to Predict Protein-Protein Interactions in Crowded Cellular Environments

Investigating protein-protein interactions is crucial for understanding cellular biological processes because proteins often function within molecular complexes rather than in isolation. While experimental and computational methods have provided valuable insights into these interactions, they often overlook a critical factor: the crowded cellular environment. This environment significantly impacts protein behavior, including structural stability, diffusion, and ultimately the nature of binding. In this review, we discuss theoretical and computational approaches that allow the modeling of biological systems to guide and complement experiments and can thus significantly advance the investigation, and possibly the predictions, of protein-protein interactions in the crowded environment of cell cytoplasm. We explore topics such as statistical mechanics for lattice simulations, hydrodynamic interactions, diffusion processes in high-viscosity environments, and several methods based on molecular dynamics simulations. By synergistically leveraging methods from biophysics and computational biology, we review the state of the art of computational methods to study the impact of molecular crowding on protein-protein interactions and discuss its potential revolutionizing effects on the characterization of the human interactome.

Intrinsic disorder and conformational coexistence in auxin coreceptors

AUXIN/INDOLE 3-ACETIC ACID (Aux/IAA) transcriptional repressor proteins and the TRANSPORT INHIBITOR RESISTANT 1/AUXIN SIGNALING F-BOX (TIR1/AFB) proteins to which they bind act as auxin coreceptors. While the structure of TIR1 has been solved, structural characterization of the regions of the Aux/IAA protein responsible for auxin perception has been complicated by their predicted disorder. Here, we use NMR, CD and molecular dynamics simulation to investigate the N-terminal domains of the Aux/IAA protein IAA17/AXR3. We show that despite the conformational flexibility of the region, a critical W-P bond in the core of the Aux/IAA degron motif occurs at a strikingly high (1:1) ratio of cis to trans isomers, consistent with the requirement of the cis conformer for the formation of the fully-docked receptor complex. We show that the N-terminal half of AXR3 is a mixture of multiple transiently structured conformations with a propensity for two predominant and distinct conformational subpopulations within the overall ensemble. These two states were modeled together with the C-terminal PB1 domain to provide the first complete simulation of an Aux/IAA. Using MD to recreate the assembly of each complex in the presence of auxin, both structural arrangements were shown to engage with the TIR1 receptor, and contact maps from the simulations match closely observations of NMR signal-decreases. Together, our results and approach provide a platform for exploring the functional significance of variation in the Aux/IAA coreceptor family and for understanding the role of intrinsic disorder in auxin signal transduction and other signaling systems.

Gaussian Process Regression Models for Predicting Atomic Energies and Multipole Moments

Developing a force field is a difficult task because its design is typically pulled in opposite directions by speed and accuracy. FFLUX breaks this trend by utilizing Gaussian process regression (GPR) to predict, at ab initio accuracy, atomic energies and multipole moments as obtained from the quantum theory of atoms in molecules (QTAIM). This work demonstrates that the in-house FFLUX training pipeline can generate successful GPR models for six representative molecules: peptide-capped glycine and alanine, glucose, paracetamol, aspirin, and ibuprofen. The molecules were sufficiently distorted to represent configurations from an AMBER-GAFF2 molecular dynamics run. All internal degrees of freedom were covered corresponding to 93 dimensions in the case of the largest molecule ibuprofen (33 atoms). Benefiting from active learning, the GPR models contain only about 2000 training points and return largely sub-kcal mol^-1 prediction errors for the validation sets. A proof of concept has been reached for transferring the model produced through active learning on one atomic property to that of the remaining atomic properties. The prediction of electrostatic interaction can be assessed at the intermolecular level, and the vast majority of interactions have a root-mean-square error of less than 0.1 kJ mol^-1 with a maximum value of ∼1 kJ mol^-1 for a glycine and paracetamol dimer.

Computational Models for the Study of Protein Aggregation

Protein aggregation has been studied by many groups around the world for many years because it can be the cause of a number of neurodegenerative diseases that have no effective treatment. Obtaining the structure of related fibrils and toxic oligomers, as well as describing the pathways and main factors that govern the self-organization process, is of paramount importance, but it is also very difficult. To solve this problem, experimental and computational methods are often combined to get the most out of each method. The effectiveness of the computational approach largely depends on the construction of a reasonable molecular model. Here we discussed different versions of the four most popular all-atom force fields AMBER, CHARMM, GROMOS, and OPLS, which have been developed for folded and intrinsically disordered proteins, or both. Continuous and discrete coarse-grained models, which were mainly used to study the kinetics of aggregation, are also summarized.

A fixed multi-site interaction charge model for an accurate prediction of the QM/MM interactions

A fixed multi-site interaction charge (FMIC) model was proposed for the accurate prediction of intermolecular electrostatic interactions based on the quantum mechanical linear response of a molecule to an external electric field. In such a model, some additional off-center interaction sites were added for capturing multipole interactions for a given molecule. By multivariate least-square fitting analysis of the calculated QM/MM interactions of a given molecule with the electrostatic environment and the electrostatic potentials of the environment at the pre-defined distributed interaction sites, the FMIC of the molecule was obtained. The model system of CO in myoglobin (Mb) was utilized to demonstrate the derivation of the FMIC. The accuracy of FMIC in predicting the electrostatic interactions between CO and the Mb environment was investigated using 10 000 different Mb-CO configurations generated from the 400 ps QM/MM MD simulation. In comparison to the QM/MM calculations at the B3LYP/aug-cc-pVTZ/ff99SB level, the mean unsigned error (MUE) of the results based on the FMIC model was merely 0.10 kcal mol^-1, and the root mean square error (RMSE) was only 0.13 kcal mol^-1, which are significantly lower than the results predicted by the ESP charge model (MUE = 1.45 kcal mol^-1, and RMSE = 1.7 kcal mol^-1, respectively). The transferability of FMIC was tested by applying the obtained FMIC in the wild type Mb-CO system to the mutants of V68F and H64L Mb-CO systems. The MUEs of the obtained results for 10 000 different configurations are both smaller than 0.2 kcal mol^-1 for the V68F and H64L Mb-CO systems in comparison to the B3LYP/aug-cc-pVTZ/ff99SB calculations, and the RMSEs are also lower than 0.2 kcal mol^-1 for both mutants. The applications of FMIC were extended to model the electrostatic interactions between a hydrogen fluoride molecule and 492 waters in a truncated octahedron box; our study showed that the FMIC could give satisfactory results with a MUE of 0.12 kcal mol^-1 and a RMSE of 0.16 kcal mol^-1 in comparison to the B3LYP/aug-cc-pVDZ/TIP3P calculations for 10 000 different configurations generated using the 10 ns classical MD simulation. Therefore, the FMIC method provides an accurate and efficient tool for predicting intermolecular electrostatic interactions, which can be utilized in the future development of molecular force fields.

Site-Specific Backbone and Side-Chain Contributions to Thermodynamic Stabilizing Forces of the WW Domain

The native structure of a protein is stabilized by a number of interactions such as main-chain hydrogen bonds and side-chain hydrophobic contacts. However, it has been challenging to determine how these interactions contribute to protein stability at single amino acid resolution. Here, we quantified site-specific thermodynamic stability at the molecular level to extend our understanding of the stabilizing forces in protein folding. We derived the free energy components of individual amino acid residues separately for the folding of the human Pin WW domain based on simulated structures. A further decomposition of the thermodynamic properties into contributions from backbone and side-chain groups enabled us to identify the critical residues in the secondary structure and hydrophobic core formation, without introducing physical modifications to the system as in site-directed mutagenesis methods. By relating the structural and thermodynamic changes upon folding for each residue, we find that the simultaneous formation of the backbone hydrogen bonds and side-chain contacts cooperatively stabilizes the folded structure. The identification of stabilizing interactions in a folding protein at atomic resolution will provide molecular insights into understanding the origin of the protein structure and into engineering a more stable protein.

Polarizable and non-polarizable force fields: Protein folding, unfolding, and misfolding

Molecular dynamics simulations are an invaluable tool to characterize the dynamic motions of proteins in atomistic detail. However, the accuracy of models derived from simulations inevitably relies on the quality of the underlying force field. Here, we present an evaluation of current non-polarizable and polarizable force fields (AMBER ff14SB, CHARMM 36m, GROMOS 54A7, and Drude 2013) based on the long-standing biophysical challenge of protein folding. We quantify the thermodynamics and kinetics of the β-hairpin formation using Markov state models of the fast-folding mini-protein CLN025. Furthermore, we study the (partial) folding dynamics of two more complex systems, a villin headpiece variant and a WW domain. Surprisingly, the polarizable force field in our set, Drude 2013, consistently leads to destabilization of the native state, regardless of the secondary structure element present. All non-polarizable force fields, on the other hand, stably characterize the native state ensembles in most cases even when starting from a partially unfolded conformation. Focusing on CLN025, we find that the conformational space captured with AMBER ff14SB and CHARMM 36m is comparable, but the ensembles from CHARMM 36m simulations are clearly shifted toward disordered conformations. While the AMBER ff14SB ensemble overstabilizes the native fold, CHARMM 36m and GROMOS 54A7 ensembles both agree remarkably well with experimental state populations. In addition, GROMOS 54A7 also reproduces experimental folding times most accurately. Our results further indicate an over-stabilization of helical structures with AMBER ff14SB. Nevertheless, the presented investigations strongly imply that reliable (un)folding dynamics of small proteins can be captured in feasible computational time with current additive force fields.

Length Dependent Folding Kinetics of Alanine-Based Helical Peptides from Optimal Dimensionality Reduction

We present a computer simulation study of helix folding in alanine homopeptides (ALA)n of length n = 5, 8, 15, and 21 residues. Based on multi-microsecond molecular dynamics simulations at room temperature, we found helix populations and relaxation times increasing from about 6% and ~2 ns for ALA5 to about 60% and ~500 ns for ALA21, and folding free energies decreasing linearly with the increasing number of residues. The helix folding was analyzed with the Optimal Dimensionality Reduction method, yielding coarse-grained kinetic models that provided a detailed representation of the folding process. The shorter peptides, ALA5 and ALA8, tended to convert directly from coil to helix, while ALA15 and ALA21 traveled through several intermediates. Coarse-grained aggregate states representing the helix, coil, and intermediates were heterogeneous, encompassing multiple peptide conformations. The folding involved multiple pathways and interesting intermediate states were present on the folding paths, with partially formed helices, turns, and compact coils. Statistically, helix initiation was favored at both termini, and the helix was most stable in the central region. Importantly, we found the presence of underlying universal local dynamics in helical peptides with correlated transitions for neighboring hydrogen bonds. Overall, the structural and dynamical parameters extracted from the trajectories are in good agreement with experimental observables, providing microscopic insights into the complex helix folding kinetics.

Process, Outcomes and Possible Elimination of Aggregation with Special Reference to Heme Proteins; Likely Remediations of Proteinopathies

Protein folding is a natural phenomenon through which a linear polypeptide possessing necessary information attains three-dimension functionally active conformation. This is a complex and multistep process and therefore, the presence of several intermediary structures could be speculated as a result of protein folding. In in vivo, this folding process is governed by the assistance of other proteins called molecular chaperones and heat shock proteins. Due to the mechanism of protein folding, these intermediary structures remain major challenge for modern biology. Mutation in gene encoding amino acid can cause adverse environmental conditions which may result in misfolding of the linear polypeptide followed by the formation of aggregates and amyloidosis. Aggregation contributes to the pathophysiology of several maladies including diabetes mellitus, Huntington's and Alzheimer's disease. The propensity of native structure to form aggregated and fibrillar assemblies is a hallmark of amyloidosis. During aggregation of a protein, transition from α helix to β sheet is observed, and mainly β sheeted structure is visualised in a mature fibril. Heme proteins are very crucial for major life activities like transport of oxygen and carbon dioxide, synthesis of ATP, role in electron transport chain, and detoxification of free radicals formed during biochemical reactions. Any structural variation in the heme proteins may lead to a fatal response. Hence characterization of the folding intermediates becomes crucial. The characterization has been deciphered with the help of strong denaturants like acetonitrile and TFE. Moreover, possible role of elimination of these aggregates and prevention of protein denaturation is also discussed. Current review deals with the basic process and mechanism of the protein folding in general and the ultimate outcomes of the protein misfolding. Since Native conformation of heme proteins is essential for some vital activities as listed above, we have discussed possible prevention of denaturation and aggregation of heme proteins such as Hb, cyt c, catalase & peroxidase.

DBAASP v3: database of antimicrobial/cytotoxic activity and structure of peptides as a resource for development of new therapeutics

The Database of Antimicrobial Activity and Structure of Peptides (DBAASP) is an open-access, comprehensive database containing information on amino acid sequences, chemical modifications, 3D structures, bioactivities and toxicities of peptides that possess antimicrobial properties. DBAASP is updated continuously, and at present, version 3.0 (DBAASP v3) contains >15 700 entries (8000 more than the previous version), including >14 500 monomers and nearly 400 homo- and hetero-multimers. Of the monomeric antimicrobial peptides (AMPs), >12 000 are synthetic, about 2700 are ribosomally synthesized, and about 170 are non-ribosomally synthesized. Approximately 3/4 of the entries were added after the initial release of the database in 2014 reflecting the recent sharp increase in interest in AMPs. Despite the increased interest, adoption of peptide antimicrobials in clinical practice is still limited as a consequence of several factors including side effects, problems with bioavailability and high production costs. To assist in developing and optimizing de novo peptides with desired biological activities, DBAASP offers several tools including a sophisticated multifactor analysis of relevant physicochemical properties. Furthermore, DBAASP has implemented a structure modelling pipeline that automates the setup, execution and upload of molecular dynamics (MD) simulations of database peptides. At present, >3200 peptides have been populated with MD trajectories and related analyses that are both viewable within the web browser and available for download. More than 400 DBAASP entries also have links to experimentally determined structures in the Protein Data Bank. DBAASP v3 is freely accessible at http://dbaasp.org.

EmhR is an indole-sensing transcriptional regulator responsible for the indole-induced antibiotic tolerance in Pseudomonas fluorescens

Indole is well known as an interspecies signalling molecule to modulate bacterial physiology; however, it is not clear how the indole signal is perceived and responded to by plant growth promoting rhizobacteria (PGPR) in the rhizosphere. Here, we demonstrated that indole enhanced the antibiotic tolerance of Pseudomonas fluorescens 2P24, a PGPR well known for its biocontrol capacity. Proteomic analysis revealed that indole influenced the expression of multiple genes including the emhABC operon encoding a major multidrug efflux pump. The expression of emhABC was regulated by a TetR-family transcription factor EmhR, which was demonstrated to be an indole-responsive regulator. Molecular dynamics simulation showed that indole allosterically affected the distance between the two DNA-recognizing helices within the EmhR dimer, leading to diminished EmhR-DNA interaction. It was further revealed the EmhR ortholog in Pseudomonas syringae was also responsible for indole-induced antibiotic tolerance, suggesting this EmhR-dependent, indole-induced antibiotic tolerance is likely to be conserved among Pseudomonas species. Taken together, our results elucidated the molecular mechanism of indole-induced antibiotic tolerance in Pseudomonas species and had important implications on how rhizobacteria sense and respond to indole in the rhizosphere.

Molecular dynamics simulations identify the regions of compromised thermostability in SazCA

The present study examined the structure and dynamics of the most active and thermostable carbonic anhydrase, SazCA, probed using molecular dynamics simulations. The molecular system was described by widely used biological force-fields (AMBER, CHARMM22, CHARMM36, and OPLS-AA) in conjunction with TIP3P water model. The comparison of molecular dynamics simulation results suggested AMBER to be a suitable choice to describe the structure and dynamics of SazCA. In addition to this, we also addressed the effect of temperature on the stability of SazCA. We performed molecular dynamics simulations at 313, 333, 353, 373, and 393 K to study the relationship between thermostability and flexibility in SazCA. The amino acid residues VAL98, ASN99, GLY100, LYS101, GLU145, and HIS207 were identified as the most flexible residues from root-mean-square fluctuations. The salt bridge analysis showed that ion-pairs ASP113-LYS81, ASP115-LYS81, ASP115-LYS114, GLU144-LYS143, and GLU144-LYS206, were responsible for the compromised thermal stability of SazCA.

Comparison and Evaluation of Force Fields for Intrinsically Disordered Proteins

Molecular dynamics (MD) simulations of six upgraded empirical force fields were compared and evaluated with short peptides, intrinsically disordered proteins, and folded proteins using trajectories of 1, 1.5, 5, or 10 μs (five replicates of 200 ns, 300 ns, 1 μs, or 2 μs) for each system. Previous studies have shown that different force fields, water models, simulation methods, and parameters can affect simulation outcomes. Here, the MD simulations were done in an explicit solvent with RS-peptide, HEWL19, HIV-rev, β amyloid (Aβ)-40, Aβ-42, phosphodiesterase-γ, CspTm, and ubiquitin using ff99IDPs, ff14IDPs, ff14IDPSFF, ff03w, CHARMM36m, and CHARMM22* force fields. The IDP ensembles generated by six all-atom empirical force fields were compared against NMR data. Despite using identical starting structures and simulation parameters, ensembles obtained with different force fields exhibit significant differences in NMR RMDs, secondary structure contents, and global properties such as the radius of gyration. The intrinsically disordered protein (IDP)-specific force fields could substantially reproduce the experimental observables in force field comparison: they have the lowest error in chemical shifts and J-couplings for short peptides/proteins, reasonably well for large IDPs and reasonably well with the radius of gyration. A high population of disorderness was observed in the IDP-specific force field for the IDP ensemble with a fraction of β sheets for β-amyloids. CHARMM22* performs better for many observables; however, it still has a preference toward the helicity for short peptides. The results of β-amyloid 42 starting from two different initial structures (Aβ42^1Z0Q and Aβ42^model) were also compared with DSSP and NMR data. The results obtained with IDP-specific force fields within 2 μs simulation time are similar, even though starting from different structures. The current force fields perform equally well for folded proteins. The results of currently developed or modified force fields for IDPs are capable of enlightening the overall performance of the force field for disordered as well as folded proteins, thereby contributing to force field development.

Crowding, Sticking, and Partial Folding of GTT WW Domain in a Small Cytoplasm Model

Recent experimental data has shown that protein folding in the cytoplasm can differ from in vitro folding with respect to speed, stability, and residual structure. Here we investigate the all-atom molecular dynamics (MD) simulations of 9 copies of the model protein GTT WW domain in a small bacterial cytoplasm model using three force fields. GTT has been well-studied by MD in aqueous solution for comparison. We find that folded copies remain folded for up 25 μs, whereas unfolded copies do not fold for up to 190 μs. Unfolded GTT in our cytoplasm model does populate partly folded intermediates with one of the two hairpins formed. Relative to aqueous solution, GTT gets stuck in metastable states with a small RMSD and radius of gyration and extensive burial of surface area against other macromolecules. In particular, GTT is even able to form transient intermolecular β-sheets with other proteins, resulting in a "chimeric structure" that could be a precursor to oligomeric β-aggregates. We conclude that sticking, enhanced by the non-native mutations of GTT, is largely responsible, and we propose, on the basis of our result as well as recent experiments, that coevolution of protein surfaces with their solvation environment (including chaperones) is important for folding and diffusion of proteins in the cytoplasm.

Pairwise-additive and polarizable atomistic force fields for molecular dynamics simulations of proteins

Protein force fields have been undergoing continual development since the first complete parameter sets were introduced nearly four decades ago. The functional forms that underlie these models have many common elements for the treatment of bonded and nonbonded forces, which are reviewed here. The most widely used force fields to date use a fixed-charge convention in which electronic polarization effects are treated via a mean-field approximation during partial charge assignment. Despite success in modeling folded proteins over many years, the fixed-charge assumption has limitations that cannot necessarily be overcome within their potential energy equations. To overcome these limitations, several force fields have recently been derived that explicitly treat electronic polarization effects with straightforward extensions of the potential energy functions used by nonpolarizable force fields. Here, we review the history of the most popular nonpolarizable force fields (AMBER, CHARMM, OPLS, and GROMOS) as well as studies that have validated them and applied them to studies of protein folding and misfolding. Building upon these force fields are more recent polarizable interaction potentials, including fluctuating charge models, POSSIM, AMOEBA, and the classical Drude oscillator. These force fields differ in their implementations but all attempt to model electronic polarization in a computationally tractable manner. Despite their recent emergence in the field of protein folding, several studies have already applied these polarizable models to challenging problems in this domain, including the role of polarization in folding free energies and sequence-specific effects on the stability of α-helical structures.

Exploration of the Misfolding Mechanism of Transthyretin Monomer: Insights from Hybrid-Resolution Simulations and Markov State Model Analysis

Misfolding and aggregation of transthyretin (TTR) is widely known to be responsible for a progressive systemic disorder called amyloid transthyretin (ATTR) amyloidosis. Studies suggest that TTR aggregation is initiated by a rate-limiting dissociation of the homo-tetramer into its monomers, which can rapidly misfold and self-assemble into amyloid fibril. Thus, exploring conformational change involved in TTR monomer misfolding is of vital importance for understanding the pathogenesis of ATTR amyloidosis. In this work, microsecond timescale hybrid-resolution molecular dynamics (MD) simulations combined with Markov state model (MSM) analysis were performed to investigate the misfolding mechanism of the TTR monomer. The results indicate that a macrostate with partially unfolded conformations may serve as the misfolded state of the TTR monomer. This misfolded state was extremely stable with a very large equilibrium probability of about 85.28%. With secondary structure analysis, we found the DAGH sheet in this state to be significantly destroyed. The CBEF sheet was relatively stable and sheet structure was maintained. However, the F-strand in this sheet was likely to move away from E-strand and reform a new β-sheet with the H-strand. This observation is consistent with experimental finding that F and H strands in the outer edge drive the misfolding of TTR. Finally, transition pathways from a near native state to this misfolded macrostate showed that the conformational transition can occur either through a native-like β-sheet intermediates or through partially unfolded intermediates, while the later appears to be the main pathway. As a whole, we identified a potential misfolded state of the TTR monomer and elucidated the misfolding pathway for its conformational transition. This work can provide a valuable theoretical basis for understanding of TTR aggregation and the pathogenesis of ATTR amyloidosis at the atomic level.

Accurate Biomolecular Simulations Account for Electronic Polarization

In this perspective, we discuss where and how accounting for electronic many-body polarization affects the accuracy of classical molecular dynamics simulations of biomolecules. While the effects of electronic polarization are highly pronounced for molecules with an opposite total charge, they are also non-negligible for interactions with overall neutral molecules. For instance, neglecting these effects in important biomolecules like amino acids and phospholipids affects the structure of proteins and membranes having a large impact on interpreting experimental data as well as building coarse grained models. With the combined advances in theory, algorithms and computational power it is currently realistic to perform simulations with explicit polarizable dipoles on systems with relevant sizes and complexity. Alternatively, the effects of electronic polarization can also be included at zero additional computational cost compared to standard fixed-charge force fields using the electronic continuum correction, as was recently demonstrated for several classes of biomolecules.

Planarity and out-of-plane vibrational modes of tryptophan and tyrosine in biomolecular modeling

Tryptophan and tyrosine are amino acids that play significant roles in the folding processes of proteins at water-membrane interfaces because of their amphipathic heteroaromatic rings. Employing appropriate heteroaromatic molecular structures is essential for obtaining accurate dynamics and predictive capabilities in molecular simulations of these amino acids. In this study, molecular dynamics simulations that applied the most recent version of the CHARMM36 force field were conducted on aqueous solutions of tryptophan and of tyrosine. Geometric analysis and dynamics quantified how aromatic rings deviated from planar structures and exhibited out-of-plane fluctuations. Radial distribution functions showed possible biological significance because the extent of ring planarity slightly affected local water concentrations near aromatic rings. Instantaneous all-atom normal mode analysis (NMA) and Fourier transformation of time autocorrelation functions of out-of-plane displacements were applied to study out-of-plane vibrations of atoms in these rings. The NMA started with minimum energy configurations and then averaged over fluctuations in aqueous solution. The frequencies and frequency patterns that were obtained for tryptophan and tyrosine with CHARMM36 differed from literature reports of Raman spectra, infrared spectra, and frequencies calculated using quantum mechanics, with some out-of-plane modes found at higher frequencies. Effects of imposing improper torsion potentials and changing torsion angle force constants were investigated for all atoms in the rings of tryptophan and tyrosine. Results show that these coarse force field variations only affect planarity and out-of-plane vibrations of atoms within the rings, and not other vibrations. Although increasing improper torsion force constants reduced deviations from aromatic ring planarity significantly, it increased out-of-plane mode frequencies. Reducing torsion angle force constants (with and without improper torsions) shifted modes to lower frequencies. A combination of decreasing most torsion angle force constants for ring atoms in both amino acids and including improper torsion forces attained frequencies and frequency patterns for out-of-plane normal modes that were more similar to the literature spectra. These force field variations decreased the extents of out-of-plane vibrations within the heteroaromatic rings of tryptophan, especially around the nitrogen atom in the ring, but not within the heteroaromatic ring of tyrosine. Conclusions were unaffected by the peptide endgroup, water, or simulation ensemble.

Applications of Molecular Dynamics Simulation in Structure Prediction of Peptides and Proteins

Compared with rapid accumulation of protein sequences from high-throughput DNA sequencing, obtaining experimental 3D structures of proteins is still much more difficult, making protein structure prediction (PSP) potentially very useful. Currently, a vast majority of PSP efforts are based on data mining of known sequences, structures and their relationships (informatics-based). However, if closely related template is not available, these methods are usually much less reliable than experiments. They may also be problematic in predicting the structures of naturally occurring or designed peptides. On the other hand, physics-based methods including molecular dynamics (MD) can utilize our understanding of detailed atomic interactions determining biomolecular structures. In this mini-review, we show that all-atom MD can predict structures of cyclic peptides and other peptide foldamers with accuracy similar to experiments. Then, some notable successes in reproducing experimental 3D structures of small proteins through MD simulations (some with replica-exchange) of the folding were summarized. We also describe advancements of MD-based refinement of structure models, and the integration of limited experimental or bioinformatics data into MD-based structure modeling.

Factors Stabilizing β-Sheets in Protein Structures from a Quantum-Chemical Perspective

Protein folds are determined by the interplay between various (de)stabilizing forces, which can be broadly divided into a local strain of the protein chain and intramolecular interactions. In contrast to the α-helix, the β-sheet secondary protein structure is significantly stabilized by long-range interactions between the individual β-strands. It has been observed that quite diverse amino acid sequences can form a very similar small β-sheet fold, such as in the three-β-strand WW domain. Employing "calibrated" quantum-chemical methods, we show herein on two sequentially diverse examples of the WW domain that the internal strain energy is higher in the β-strands and lower in the loops, while the interaction energy has an opposite trend. Low strain energy computed for peptide sequences in the loop 1 correlates with its postulated early formation in the folding process. The relatively high strain energy within the β-strands (up to 8 kcal mol^-1 per amino acid residue) is compensated by even higher intramolecular interaction energy (up to 15 kcal mol^-1 per residue). It is shown in a quantitative way that the most conserved residues across the structural family of WW domains have the highest contributions to the intramolecular interaction energy. On the other hand, the residues in the regions with the lowest strain are not conserved. We conclude that the internal interaction energy is the physical quantity tuned by evolution to define the β-sheet protein fold.

Computational Estimation of Microsecond to Second Atomistic Folding Times

Despite the development of massively parallel computing hardware including inexpensive graphics processing units (GPUs), it has remained infeasible to simulate the folding of atomistic proteins at room temperature using conventional molecular dynamics (MD) beyond the microsecond scale. Here, we report the folding of atomistic, implicitly solvated protein systems with folding times τ ranging from ∼10 μs to ∼100 ms using the weighted ensemble (WE) strategy in combination with GPU computing. Starting from an initial structure or set of structures, WE organizes an ensemble of GPU-accelerated MD trajectory segments via intermittent pruning and replication events to generate statistically unbiased estimates of rate constants for rare events such as folding; no biasing forces are used. Although the variance among atomistic WE folding runs is significant, multiple independent runs are used to reduce and quantify statistical uncertainty. Folding times are estimated directly from WE probability flux and from history-augmented Markov analysis of the WE data. Three systems were examined: NTL9 at low solvent viscosity (yielding τf = 0.8-9 μs), NTL9 at water-like viscosity (τf = 0.2-2 ms), and Protein G at low viscosity (τf = 3-200 ms). In all cases, the folding time, uncertainty, and ensemble properties could be estimated from WE simulation; for Protein G, this characterization required significantly less overall computing than would be required to observe a single folding event with conventional MD simulations. Our results suggest that the use and calibration of force fields and solvent models for precise estimation of kinetic quantities is becoming feasible.

Rotamer Dynamics: Analysis of Rotamers in Molecular Dynamics Simulations of Proteins

Given by χ torsional angles, rotamers describe the side-chain conformations of amino acid residues in a protein based on the rotational isomers (hence the word rotamer). Constructed rotamer libraries, based on either protein crystal structures or dynamics studies, are the tools for classifying rotamers (torsional angles) in a way that reflect their frequency in nature. Rotamer libraries are routinely used in structure modeling and evaluation. In this perspective article, we would like to encourage researchers to apply rotamer analyses beyond their traditional use. Molecular dynamics (MD) of proteins highlight the in silico behavior of molecules in solution and thus can identify favorable side-chain conformations. In this article, we used simple computational tools to study rotamer dynamics (RD) in MD simulations. First, we isolated each frame in the MD trajectories in separate Protein Data Bank files via the cpptraj module in AMBER. Then, we extracted torsional angles via the Bio3D module in R language. The classification of torsional angles was also done in R according to the penultimate rotamer library. RD analysis is useful for various applications such as protein folding, study of rotamer-rotamer relationship in protein-protein interaction, real-time correlation between secondary structures and rotamers, study of flexibility of side chains in binding site for molecular docking preparations, use of RD as guide in functional analysis and study of structural changes caused by mutations, providing parameters for improving coarse-grained MD accuracy and speed, and many others. Major challenges facing RD to emerge as a new scientific field involve the validation of results via easy, inexpensive wet-lab methods. This realm is yet to be explored.

Atomistic Force Fields for Proteins

All-atom, classical force fields for protein molecular dynamics (MD) simulations currently occupy a sweet spot in the universe of computational models, sufficiently detailed to be of predictive value in many cases, yet also simple enough that some biologically relevant time scales (microseconds or more) can now be sampled via specialized hardware or enhanced sampling methods. However, due to their long evolutionary history, there is now a myriad of force field branches in current use, which can make it hard for those entering the simulation field to know which would be the best set of parameters for a given application. In this chapter, I try to give an overview of the historical motivation for the different force fields available, suggestions for how to determine the most appropriate model and what to do if the results are in conflict with experimental evidence.

Folding Molecular Dynamics Simulation of a gp41-Derived Peptide Reconcile Divergent Structure Determinations

T-20 peptide is the first FDA-approved fusion inhibitor against AIDS/HIV-1 gp41 protein. Part of it, the gp41[659-671] peptide, that contains the complete epitope for the neutralizing 2F5 monoclonal antibody, has been found experimentally in a number of divergent structures. Herein, we attempt to reconcile them by using unbiased large-scale all-atom molecular dynamics folding simulations. We show that our approach can successfully capture the peptide's heterogeneity and reach each and every experimentally determined conformation in sub-angstrom accuracy, whilst preserving the peptide's disordered nature. Our analysis also unveils that the minor refinements within the AMBER99SB family of force fields can lead to appreciable differences in the predicted conformational stability arising from subtle differences in the helical- and β-region of the Ramachandran plot. Our work underlines the contribution of molecular dynamics simulation in structurally characterizing pharmacologically important peptides of ambiguous structure.

A spin-1 representation for dual-funnel energy landscapes

The interconversion between the left- and right-handed helical folds of a polypeptide defines a dual-funneled free energy landscape. In this context, the funnel minima are connected through a continuum of unfolded conformations, evocative of the classical helix-coil transition. Physical intuition and recent conjectures suggest that this landscape can be mapped by assigning a left- or right-handed helical state to each residue. We explore this possibility using all-atom replica exchange molecular dynamics and an Ising-like model, demonstrating that the energy landscape architecture is at odds with a two-state picture. A three-state model-left, right, and unstructured-can account for most key intermediates during chiral interconversion. Competing folds and excited conformational states still impose limitations on the scope of this approach. However, the improvement is stark: Moving from a two-state to a three-state model decreases the fit error from 1.6 kBT to 0.3 kBT along the left-to-right interconversion pathway.

Experimental and Theoretical Insights into the Inhibition Mechanism of Prion Fibrillation by Resveratrol and its Derivatives

Resveratrol and its derivatives have been shown to display beneficial effects to neurodegenerative diseases. However, the molecular mechanism of resveratrol and its derivatives on prion conformational conversion is poorly understood. In this work, the interaction mechanism between prion and resveratrol as well as its derivatives was investigated using steady-state fluorescence quenching, Thioflavin T binding assay, Western blotting, and molecular dynamics simulation. Protein fluorescence quenching method and Thioflavin T assay revealed that resveratrol and its derivatives could interact with prion and interrupt prion fibril formation. Molecular dynamics simulation results indicated that resveratrol can stabilize the PrP_127-147 peptide mainly through π-π stacking interactions between resveratrol and Tyr128. The hydrogen bonds interactions between resveratrol and the PrP_127-147 peptide could further reduce the flexibility and the propensity to aggregate. The results of this study not only can provide useful information about the interaction mechanism between resveratrol and prion, but also can provide useful clues for further design of new inhibitors inhibiting prion aggregation.

Deficiencies in Molecular Dynamics Simulation-Based Prediction of Protein-DNA Binding Free Energy Landscapes

Transcriptional regulation allows cells to match their gene expression profiles to their current requirements based on environment, cellular physiological state, and extracellular signals. DNA binding transcription factors are major agents of transcriptional regulation, and bind to DNA with a factor-specific sequence preference to exert regulatory effects. A crucial step in unraveling the logic of a regulatory network is determining the sequence-specific binding affinity landscapes for the transcription factors in it. While such landscapes can be measured experimentally, the ability to predict them computationally would both reduce the effort required to obtain the needed data and provide additional insight into the key interactions shaping protein-DNA interactions. Here we apply free energy calculations based on all-atom molecular dynamics simulations to predict the changes in binding free energy for all single base pair perturbations of the binding sites for four eukaryotic transcription factors for which high-quality experimental data exist. We find that the simulated results both vastly overestimate the magnitude of changes in binding free energy, and frequently predict the incorrect signs. These simulations will nevertheless serve as a jumping-off point for refining our current representation of protein-DNA interactions to allow quantitative reproduction of experimental data on such systems in the future.

Computer Simulations of Intrinsically Disordered Proteins

The investigation of intrinsically disordered proteins (IDPs) is a new frontier in structural and molecular biology that requires a new paradigm to connect structural disorder to function. Molecular dynamics simulations and statistical thermodynamics potentially offer ideal tools for atomic-level characterizations and thermodynamic descriptions of this fascinating class of proteins that will complement experimental studies. However, IDPs display sensitivity to inaccuracies in the underlying molecular mechanics force fields. Thus, achieving an accurate structural characterization of IDPs via simulations is a challenge. It is also daunting to perform a configuration-space integration over heterogeneous structural ensembles sampled by IDPs to extract, in particular, protein configurational entropy. In this review, we summarize recent efforts devoted to the development of force fields and the critical evaluations of their performance when applied to IDPs. We also survey recent advances in computational methods for protein configurational entropy that aim to provide a thermodynamic link between structural disorder and protein activity.

Direct folding simulation of helical proteins using an effective polarizable bond force field

Snapshots of the intermediate conformation of Trp-cage at various simulation times using AMBER03, EPB03, AMBER12SB, and EPB12SB. Here, the N terminal is always on the top.

Effect of electrostatic polarization and bridging water on CDK2–ligand binding affinities calculated using a highly efficient interaction entropy method

A new highly efficient interaction entropy (IE) method combined with the polarized protein-specific charge (PPC) force field is employed to investigate the interaction mechanism of CDK2–ligand binding and the effect of the bridging water.

Conformations of a Metastable SH3 Domain Characterized by smFRET and an Excluded-Volume Polymer Model

Conformational states of the metastable drkN SH3 domain were characterized using single-molecule fluorescence techniques. Under nondenaturing conditions, two Förster resonance energy transfer (FRET) populations were observed that corresponded to a folded and an unfolded state. FRET-estimated radii of gyration and hydrodynamic radii estimated by fluorescence correlation spectroscopy of the two coexisting conformations are in agreement with previous ensemble x-ray scattering and NMR measurements. Surprisingly, when exposed to high concentrations of urea and GdmCl denaturants, the protein still exhibits two distinct FRET populations. The dominant conformation is expanded, showing a low FRET efficiency, consistent with the expected behavior of a random chain with excluded volume. However, approximately one-third of the drkN SH3 conformations showed high, nearly 100%, FRET efficiency, which is shown to correspond to denaturation-induced looped conformations that remain stable on a timescale of at least 100 μs. These loops may contain interconverting conformations that are more globally collapsed, hairpin-like, or circular, giving rise to the observed heterogeneous broadening of this population. Although the underlying mechanism of chain looping remains elusive, FRET experiments in formamide and dimethyl sulfoxide suggest that interactions between hydrophobic groups in the distal regions may play a significant role in the formation of the looped state.

A multi-resolution model to capture both global fluctuations of an enzyme and molecular recognition in the ligand-binding site

In multi-resolution simulations, different system components are simultaneously modeled at different levels of resolution, these being smoothly coupled together. In the case of enzyme systems, computationally expensive atomistic detail is needed in the active site to capture the chemistry of ligand binding. Global properties of the rest of the protein also play an essential role, determining the structure and fluctuations of the binding site; however, these can be modeled on a coarser level. Similarly, in the most computationally efficient scheme only the solvent hydrating the active site requires atomistic detail. We present a methodology to couple atomistic and coarse-grained protein models, while solvating the atomistic part of the protein in atomistic water. This allows a free choice of which protein and solvent degrees of freedom to include atomistically. This multi-resolution methodology can successfully model stable ligand binding, and we further confirm its validity by exploring the reproduction of system properties relevant to enzymatic function. In addition to a computational speedup, such an approach can allow the identification of the essential degrees of freedom playing a role in a given process, potentially yielding new insights into biomolecular function. Proteins 2016; 84:1902-1913. © 2016 Wiley Periodicals, Inc.

A Distal Disulfide Bridge in OXA-1 β-Lactamase Stabilizes the Catalytic Center and Alters the Dynamics of the Specificity Determining Ω Loop

Widespread antibiotic resistance, particularly when mediated by broad-spectrum β-lactamases, has major implications for public health. Substitutions in the active site often allow broad-spectrum enzymes to accommodate diverse types of β-lactams. Substitutions observed outside the active site are thought to compensate for the loss of thermal stability. The OXA-1 clade of class D β-lactamases contains a pair of conserved cysteines located outside the active site that forms a disulfide bond in the periplasm. Here, the effect of the distal disulfide bond on the structure and dynamics of OXA-1 was investigated via 4 μs molecular dynamics simulations. The results reveal that the disulfide promotes the preorganized orientation of the catalytic residues and affects the conformation of the functionally important Ω loop. Furthermore, principal component analysis reveals differences in the global dynamics between the oxidized and reduced forms, especially in the motions involving the Ω loop. A dynamical network analysis indicates that, in the oxidized form, in addition to its role in ligand binding, the KTG family motif is a central hub of the global dynamics. As activity of OXA-1 has been measured only in the reduced form, we suggest that accurate assessment of its functional profile would require oxidative conditions mimicking periplasm.

The inhibitory mechanism of a fullerene derivative against amyloid-β peptide aggregation: an atomistic simulation study

Alzheimer's disease (AD) is associated with the pathological self-assembly of amyloid-β (Aβ) peptides into β-sheet enriched fibrillar aggregates. Aβ dimers formed in the initial step of Aβ aggregation were reported to be the smallest toxic species. Inhibiting the formation of β-sheet-rich oligomers and fibrils is considered as the primary therapeutic strategy for AD. Previous studies reported that fullerene derivatives strongly inhibit Aβ fibrillation. However, the underlying inhibitory mechanism remains elusive. As a first step to understand fullerene-modulated full-length Aβ aggregation, we investigated the conformational ensemble of the Aβ1-42 dimer with and without 1,2-(dimethoxymethano)fullerene (DMF) - a more water-soluble fullerene derivative - by performing a 340 ns explicit-solvent replica exchange molecular dynamics simulation. Our simulations show that although disordered states are the most abundant conformations of the Aβ1-42 dimer, conformations containing diverse extended β-hairpins are also populated. The first most-populated β-hairpins involving residues L17-D23 and A30-V36 strongly resemble the engineered β-hairpin which is a building block of toxic Aβ oligomers. We find that the interaction of DMFs with Aβ peptides greatly impedes the formation of such β-hairpins and inter-peptide β-sheets. Binding energy analyses demonstrate that DMF preferentially binds not only to the central hydrophobic motif LVFFA of the Aβ peptide as suggested experimentally, but also to the aromatic residues including F4 and Y10 and the C-terminal hydrophobic region I31-V40. This study reveals a complete picture of the inhibitory mechanism of full-length Aβ1-42 aggregation by fullerenes, providing theoretical insights into the development of drug candidates against AD.

Principles and Overview of Sampling Methods for Modeling Macromolecular Structure and Dynamics

Investigation of macromolecular structure and dynamics is fundamental to understanding how macromolecules carry out their functions in the cell. Significant advances have been made toward this end in silico, with a growing number of computational methods proposed yearly to study and simulate various aspects of macromolecular structure and dynamics. This review aims to provide an overview of recent advances, focusing primarily on methods proposed for exploring the structure space of macromolecules in isolation and in assemblies for the purpose of characterizing equilibrium structure and dynamics. In addition to surveying recent applications that showcase current capabilities of computational methods, this review highlights state-of-the-art algorithmic techniques proposed to overcome challenges posed in silico by the disparate spatial and time scales accessed by dynamic macromolecules. This review is not meant to be exhaustive, as such an endeavor is impossible, but rather aims to balance breadth and depth of strategies for modeling macromolecular structure and dynamics for a broad audience of novices and experts.