Folding Thermodynamics and Mechanism of Five Trp-Cage Variants from Replica-Exchange MD Simulations with RSFF2 Force Field

Effect of Main and Side Chains on the Folding Mechanism of the Trp-Cage Miniprotein

Proteins that do not fold into their functional native state have been linked to diseases. In this study, the influence of the main and side chains of individual amino acids on the folding of the tryptophan cage (Trp-cage), a designed 20-residue miniprotein, was analyzed. For this purpose, we calculated the solvation free energy (SFE) contributions of individual atoms by using the 3D-reference interaction site model with the atomic decomposition method. The mechanism by which the Trp-cage is stabilized during the folding process was examined by calculating the total energy, which is the sum of the conformational energy and SFE. The folding process of the Trp-cage resulted in a stable native state, with a total energy that was 62.4 kcal/mol lower than that of the unfolded state. The solvation entropy, which is considered to be responsible for the hydrophobic effect, contributed 31.3 kcal/mol to structural stabilization. In other words, the contribution of the solvation entropy accounted for approximately half of the total contribution to Trp-cage folding. The hydrophobic core centered on Trp6 contributed 15.6 kcal/mol to the total energy, whereas the solvation entropy contribution was 6.3 kcal/mol. The salt bridge formed by the hydrophilic side chains of Asp9 and Arg16 contributed 10.9 and 5.0 kcal/mol, respectively. This indicates that not only the hydrophobic core but also the salt bridge of the hydrophilic side chains gain solvation entropy and contribute to stabilizing the native structure of the Trp-cage.

Computational Modeling of the Thermodynamics of the Mesophilic and Thermophilic Mutants of Trp-Cage Miniprotein

We characterize the folding-unfolding thermodynamics of two mutants of the miniprotein Trp-cage by combining extended molecular dynamics simulations and an advanced statistical-mechanical-based approach. From a set of molecular dynamics simulations in an explicit solvent performed along a reference isobar, we evaluated the structural and thermodynamic behaviors of a mesophilic and a thermophilic mutant of the Trp-cage and their temperature dependence. In the case of the thermophilic mutant, computational data confirm that our theoretical-computational approach is able to reproduce the available experimental estimate with rather good accuracy. On the other hand, the mesophilic mutant does not show a clear two-state (folded and unfolded) behavior, preventing us from reconstructing its thermodynamics; thus, an analysis of its structural behavior along a reference isobar is presented. Our results show that an extended sampling of these kinds of systems coupled to an advanced statistical-mechanical-based treatment of the data can provide an accurate description of the folding-unfolding thermodynamics along a reference isobar, rationalizing the discrepancies between the simulated and experimental systems.

Cyclic peptides: backbone rigidification and capability of mimicking motifs at protein-protein interfaces

Cyclization is commonly employed in efforts to improve the target binding affinity of peptide-based probes and therapeutics. Many structural motifs have been identified at protein-protein interfaces and provide promising targets for inhibitor design using cyclic peptides. Cyclized peptides are generally assumed to be rigidified relative to their linear counterparts. This rigidification potentially pre-organizes the molecules to interact properly with their targets. However, the actual impact of cyclization on, for example, peptide configurational entropy, is currently poorly understood in terms of both its magnitude and molecular-level origins. Moreover, even with thousands of desired structural motifs at hand, it is currently not possible to a priori identify the ones that are most promising to mimic using cyclic peptides nor to select the ideal linker length. Instead, labor-intensive chemical synthesis and experimental characterization of various cyclic peptide designs are required, in hopes of finding one with improved target affinity. Herein, using molecular dynamics simulations of polyglycines, we elucidated how head-to-tail cyclization impacts peptide backbone dihedral entropy and developed a simple strategy to rapidly screen for structures that can be reliably mimicked by preorganized cyclic peptides. As expected, cyclization generally led to a reduction in backbone dihedral entropy; notably, however, this effect was minimal when the length of polyglycines was >9 residues. We also found that the reduction in backbone dihedral entropy upon cyclization of small polyglycine peptides does not result from more restricted distributions of the dihedrals; rather, it was the correlations between specific dihedrals that caused the decrease in configurational entropy in the cyclic peptides. Using our comprehensive cyclo-Gn structural ensembles, we obtained a holistic picture of what conformations are accessible to cyclic peptides. Using "hot loops" recently identified at protein-protein interfaces as an example, we provide clear guidelines for choosing the "easiest" hot loops for cyclic peptides to mimic and for identifying appropriate cyclic peptide lengths. In conclusion, our results provide an understanding of the thermodynamics and structures of this interesting class of molecules. This information should prove particularly useful for designing cyclic peptide inhibitors of protein-protein interactions.

Increasing the Sampling Efficiency of Protein Conformational Change by Combining a Modified Replica Exchange Molecular Dynamics and Normal Mode Analysis

Understanding conformational change at an atomic level is significant when determining a protein functional mechanism. Replica exchange molecular dynamics (REMD) is a widely used enhanced sampling method to explore protein conformational space. However, REMD with an explicit solvent model requires huge computational resources, immensely limiting its application. In this study, a variation of parallel tempering metadynamics (PTMetaD) with the omission of solvent-solvent interactions in exchange attempts and the use of low-frequency modes calculated by normal-mode analysis (NMA) as collective variables (CVs), namely ossPTMetaD, is proposed with the aim to accelerate MD simulations simultaneously in temperature and geometrical spaces. For testing the performance of ossPTMetaD, five protein systems with diverse biological functions and motion patterns were selected, including large-scale domain motion (AdK), flap movement (HIV-1 protease and BACE1), and DFG-motif flip in kinases (p38α and c-Abl). The simulation results showed that ossPTMetaD requires much fewer numbers of replicas than temperature REMD (T-REMD) with a reduction of ∼70% to achieve a similar exchange ratio. Although it does not obey the detailed balance condition, ossPTMetaD provides consistent results with T-REMD and experimental data. The high accessibility of the large conformational change of protein systems by ossPTMetaD, especially in simulating the very challenging DFG-motif flip of protein kinases, demonstrated its high efficiency and robustness in the characterization of the large-scale protein conformational change pathway and associated free energy profile.

Unperturbed Detection of the Dynamic Structure in the Hydrophobic Core of Trp-Cage via Two-Dimensional Infrared Spectroscopy

The tyrosine ring mode is an intrinsic non-perturbing site-specific infrared reporter for conformational dynamics within protein systems. This transition is influenced by direct and indirect interactions associated with the electron-donating ability and the hydrophobicity of the surrounding molecules. Utilizing an intrinsic tyrosine moiety, two-dimensional infrared spectra of Trp-cage, often called the "hydrogen atom" of protein folding, were measured in the folded and denatured states to uncover the dynamics of the hydrophobic core. The vibrational lifetimes and the correlation decays of the tyrosine ring mode showed significant changes upon both temperature and chemical denaturation of the Trp-cage miniprotein, indicating important structural features of the hydrophobic core and its dynamics. The observed Trp6-Tyr3 interactions are in good agreement with the prior studies of the folded state, but they reach beyond the static structure. These stacking interactions and orientations fluctuate on the picosecond time scale as measured through the spectral dephasing within a dehydrated environment.

Exploring Conformational Change of Adenylate Kinase by Replica Exchange Molecular Dynamic Simulation

Replica exchange molecular dynamics (REMD) simulation is a popular enhanced sampling method that is widely used for exploring the atomic mechanism of protein conformational change. However, the requirement of huge computational resources for REMD, especially with the explicit solvent model, largely limits its application. In this study, the availability and efficiency of a variant of velocity-scaling REMD (vsREMD) was assessed with adenylate kinase as an example. Although vsREMD achieved results consistent with those from conventional REMD and experimental studies, the number of replicas required for vsREMD (30) was much less than that for conventional REMD (80) to achieve a similar acceptance rate (∼0.2), demonstrating high efficiency of vsREMD to characterize the protein conformational change and associated free-energy profile. Thus, vsREMD is a highly efficient approach for studying the large-scale conformational change of protein systems.

Theoretical-computational modelling of the temperature dependence of the folding-unfolding thermodynamics and kinetics: the case of a Trp-cage

Here we present a theoretical-computational study of the thermodynamics and kinetics of an aqueous Trp-cage, a 20-residue long miniprotein. The combined use of accurate molecular dynamics simulations rigorously reconstructing the proper isobar of the system and a sound statistical-mechanical model provides a quantitative description of the temperature dependence of the relevant physical-chemical properties and insights into the detailed mechanisms regulating the folding-unfolding properties.

Accelerated Molecular Dynamics Simulation for Helical Proteins Folding in Explicit Water

In this study, we examined the folding processes of eight helical proteins (2I9M, TC5B, 1WN8, 1V4Z, 1HO2, 1HLL, 2KFE, and 1YYB) at room temperature using the explicit solvent model under the AMBER14SB force field with the accelerated molecular dynamics (AMD) and traditional molecular dynamics (MD), respectively. We analyzed and compared the simulation results obtained by these two methods based on several aspects, such as root mean square deviation (RMSD), native contacts, cluster analysis, folding snapshots, free energy landscape, and the evolution of the radius of gyration, which showed that these eight proteins were successfully and consistently folded into the corresponding native structures by AMD simulations carried out at room temperature. In addition, the folding occurred in the range of 40~180 ns after starting from the linear structures of the eight proteins at 300 K. By contrast, these stable folding structures were not found when the traditional molecular dynamics (MD) simulation was used. At the same time, the influence of high temperatures (350, 400, and 450 K) is also further investigated. Study found that the simulation efficiency of AMD is higher than that of MD simulations, regardless of the temperature. Of these temperatures, 300 K is the most suitable temperature for protein folding for all systems. To further investigate the efficiency of AMD, another trajectory was simulated for eight proteins with the same linear structure but different random seeds at 300 K. Both AMD trajectories reached the correct folded structures. Our result clearly shows that AMD simulation are a highly efficient and reliable method for the study of protein folding.

Long Time-Scale Atomistic Simulations of the Structure and Dynamics of Transcription Factor-DNA Recognition

Recent years have witnessed an explosion of interest in computational studies of DNA binding proteins, including both coarse-grained and atomistic simulations of transcription factor-DNA recognition, to understand how these transcription factors recognize their binding sites on the DNA with such exquisite specificity. The present study performs microsecond time scale all-atom simulations of the dimeric form of the lactose repressor (LacI), both in the absence of any DNA and in the presence of both specific and nonspecific complexes, considering three different DNA sequences. We examine, specifically, the conformational differences between specific and nonspecific protein-DNA interactions, as well as the behavior of the helix-turn-helix motif of LacI when interacting with the DNA. Our simulations suggest that stable LacI binding occurs primarily to bent A-form DNA, with a loss of LacI conformational entropy and optimization of correlated conformational equilibria across the protein. In addition, binding to the specific operator sequence involves a slightly larger number of stabilizing DNA-protein hydrogen bonds (in comparison to nonspecific complexes), which may account for the experimentally observed specificity for this operator. In doing so, our simulations provide a detailed atomistic description of potential structural drivers for LacI selectivity.

Hybrid Cascade-Type Molecular Dynamics with a Markov State Model for Efficient Free Energy Calculations

A protocol for calculating free energy landscapes (FELs) is proposed based on a combination of two cascade-type molecular dynamics (MD) methods, parallel cascade selection MD (PaCS-MD) and outlier flooding method (OFLOOD), with the help of a Markov state model (MSM). The former rapidly generates approximated transition paths directly connecting reactants with products, and the latter complementary resamples marginal conformational subspaces. Trajectories obtained by them give reliable microstates in MSM providing accurate FEL with low computational costs. As a demonstration, the present method was applied to a miniprotein (Chignolin and Trp-cage) in explicit water and successfully elucidated multiple folding paths on their free energy landscapes. Our method could be applicable to a wide variety of biological systems to estimate their free energy profiles.

Flexible selection of the solute region in replica exchange with solute tempering: Application to protein-folding simulations

Replica-exchange molecular dynamics (REMD) and their variants have been widely used in simulations of the biomolecular structure and dynamics. Replica exchange with solute tempering (REST) is one of the methods where temperature of a pre-defined solute molecule is exchanged between replicas, while solvent temperatures in all the replicas are kept constant. REST greatly reduces the number of replicas compared to the temperature REMD, while replicas at low temperatures are often trapped under their conditions, interfering with the conformational sampling. Here, we introduce a new scheme of REST, referred to as generalized REST (gREST), where the solute region is defined as a part of a molecule or a part of the potential energy terms, such as the dihedral-angle energy term or Lennard-Jones energy term. We applied this new method to folding simulations of a β-hairpin (16 residues) and a Trp-cage (20 residues) in explicit water. The protein dihedral-angle energy term is chosen as the solute region in the simulations. gREST reduces the number of replicas necessary for good random walks in the solute-temperature space and covers a wider conformational space compared to the conventional REST2. Considering the general applicability, gREST should become a promising tool for the simulations of protein folding, conformational dynamics, and an in silico drug design.

Protein dynamics and structural waters in bromodomains

Bromodomains are epigenetic readers of acetylated lysines that are integral parts of histone tails. The 61 bromodomains in humans are structurally highly conserved but specifically bind to widely varying recognition motifs, suggesting that dynamic rather than static factors are responsible for recognition selectivity. To test this hypothesis, the dynamics of the binding sites and structural water molecules of four bromodomains (ATAD2, BAZ2B, BRD2(1) and CREBBP) representing four different subtypes is studied with 1 μs MD simulations using the RSFF2 force field. The different dynamics of the ZA-loops and BC-loops between the four bromodomains leads to distinct patterns for the opening and closing of the binding pocket. This in turn determines the structural and energetic properties of the structural waters in the binding pocket, suggesting that these waters are not only important for the recognition itself, as has been proposed previously, but also contribute to the selectivity of different bromodomains.

Temperature evolution of Trp-cage folding pathways: An analysis by dividing the probability flux field into stream tubes

Owing to its small size and very fast folding rate, the Trp-cage miniprotein has become a benchmark system to study protein folding. Two folding pathways were found to be characteristic of this protein: pathway I, in which the hydrophobic collapse precedes the formation of α-helix, and pathway II, in which the events occur in the reverse order. At the same time, the relative contribution of these pathways at different temperatures as well as the nature of transition from one pathway to the other remain unclear. To gain insight into this issue, we employ a recently proposed hydrodynamic description of protein folding, in which the process of folding is considered as a motion of a "folding fluid" (Chekmarev et al., Phys. Rev. Lett. 100(1), 018107 2008). Using molecular dynamics simulations, we determine the field of probability fluxes of transitions in a space of collective variables and divide it into stream tubes. Each tube contains a definite fraction of the total folding flow and can be associated with a certain pathway. Specifically, three temperatures were considered, T = 285K, T = 315K, and T = 325K. We have found that as the temperature increases, the contribution of pathway I, which is approximately 90% of the total folding flow at T = 285K, decreases to approximately 10% at T = 325K, i.e., pathway II becomes dominant. At T = 315K, both pathways contribute approximately equally. All these temperatures are found below the calculated melting point, which suggests that the Trp-cage folding mechanism is determined by kinetic factors rather than thermodynamics.

Thermodynamics of Protein Folding Studied by Umbrella Sampling along a Reaction Coordinate of Native Contacts

Spontaneous transitions between the native and non-native protein conformations are normally rare events that hardly take place in typical unbiased molecular dynamics simulations. It was recently demonstrated that such transitions can be well described by a reaction coordinate, Q, that represents the collective fraction of the native contacts between the protein atoms. Here we attempt to use this reaction coordinate to enhance the conformational sampling. We perform umbrella sampling simulations with biasing potentials on Q for two model proteins, Trp-Cage and BBA, using the CHARMM force field. Hamiltonian replica exchange is implemented in these simulations to further facilitate the sampling. The simulations appear to have reached satisfactory convergence, resulting in unbiased free energies as a function of Q. In addition to the native structure, multiple folded conformations are identified in the reconstructed equilibrium ensemble. Some conformations without any native contacts nonetheless have rather compact geometries and are stabilized by hydrogen bonds not present in the native structure. Whereas the enhanced sampling along Q reasonably reproduces the equilibrium conformational space, we also find that the folding of an α-helix in Trp-Cage is a slow degree of freedom orthogonal to Q and therefore cannot be accelerated by biasing the reaction coordinate. Overall, we conclude that whereas Q is an excellent parameter to analyze the simulations, it is not necessarily a perfect reaction coordinate for enhanced sampling, and better incorporation of other slow degrees of freedom may further improve this reaction coordinate.

Mechanism of Phosphorylation-Induced Folding of 4E-BP2 Revealed by Molecular Dynamics Simulations

Site-specific phosphorylation of an intrinsically disordered protein, eIF4E-binding protein isoform 2 (4E-BP2), can suppress its native function by folding it into a four-stranded β-sheet, but the mechanism of this phosphorylation-induced folding is unclear. In this work, we use all-atom molecular dynamics simulations to investigate both the folded and unfolded states of 4E-BP2 under different phosphorylation states of T37 and T46. The results show that the phosphorylated forms of both T37 and T46 play important roles in stabilizing the folded structure, especially for the β-turns and the sequestered binding motif. The phosphorylated residues not only guide the folding of the protein through several intermediate states but also affect the conformational distribution of the unfolded ensemble. Significantly, the phosphorylated residues can function as nucleation sites for the folding of the protein by forming certain local structures that are stabilized by hydrogen bonding involving the phosphate group. The region around phosphorylated T46 appears to fold before that around phosphorylated T37. These findings provide new insight into the intricate effects of protein phosphorylation.

The Folding of de Novo Designed Protein DS119 via Molecular Dynamics Simulations

As they are not subjected to natural selection process, de novo designed proteins usually fold in a manner different from natural proteins. Recently, a de novo designed mini-protein DS119, with a βαβ motif and 36 amino acids, has folded unusually slowly in experiments, and transient dimers have been detected in the folding process. Here, by means of all-atom replica exchange molecular dynamics (REMD) simulations, several comparably stable intermediate states were observed on the folding free-energy landscape of DS119. Conventional molecular dynamics (CMD) simulations showed that when two unfolded DS119 proteins bound together, most binding sites of dimeric aggregates were located at the N-terminal segment, especially residues 5-10, which were supposed to form β-sheet with its own C-terminal segment. Furthermore, a large percentage of individual proteins in the dimeric aggregates adopted conformations similar to those in the intermediate states observed in REMD simulations. These results indicate that, during the folding process, DS119 can easily become trapped in intermediate states. Then, with diffusion, a transient dimer would be formed and stabilized with the binding interface located at N-terminals. This means that it could not quickly fold to the native structure. The complicated folding manner of DS119 implies the important influence of natural selection on protein-folding kinetics, and more improvement should be achieved in rational protein design.

Unfolding mechanism of thrombin-binding aptamer revealed by molecular dynamics simulation and Markov State Model

Thrombin-binding aptamer (TBA) with the sequence 5′GGTTGGTGTGGTTGG3′ could fold into G-quadruplex, which correlates with functionally important genomic regionsis. However, unfolding mechanism involved in the structural stability of G-quadruplex has not been satisfactorily elucidated on experiments so far. Herein, we studied the unfolding pathway of TBA by a combination of molecular dynamics simulation (MD) and Markov State Model (MSM). Our results revealed that the unfolding of TBA is not a simple two-state process but proceeds along multiple pathways with multistate intermediates. One high flux confirms some observations from NMR experiment. Another high flux exhibits a different and simpler unfolding pathway with less intermediates. Two important intermediate states were identified. One is similar to the G-triplex reported in the folding of G-quadruplex, but lack of H-bonding between guanines in the upper plane. More importantly, another intermediate state acting as a connector to link the folding region and the unfolding one, was the first time identified, which exhibits higher population and stability than the G-triplex-like intermediate. These results will provide valuable information for extending our understanding the folding landscape of G-quadruplex formation.

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.