Folding dynamics of the Trp-cage miniprotein: evidence for a native-like intermediate from combined time-resolved vibrational spectroscopy and molecular dynamics simulations

Unnatural Amino Acids for Biological Spectroscopy and Microscopy

Due to advances in methods for site-specific incorporation of unnatural amino acids (UAAs) into proteins, a large number of UAAs with tailored chemical and/or physical properties have been developed and used in a wide array of biological applications. In particular, UAAs with specific spectroscopic characteristics can be used as external reporters to produce additional signals, hence increasing the information content obtainable in protein spectroscopic and/or imaging measurements. In this Review, we summarize the progress in the past two decades in the development of such UAAs and their applications in biological spectroscopy and microscopy, with a focus on UAAs that can be used as site-specific vibrational, fluorescence, electron paramagnetic resonance (EPR), or nuclear magnetic resonance (NMR) probes. Wherever applicable, we also discuss future directions.

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.

An interoperable implementation of collective-variable based enhanced sampling methods in extended phase space within the OpenMM package

Collective variable (CV)-based enhanced sampling techniques are widely used today for accelerating barrier-crossing events in molecular simulations. A class of these methods, which includes temperature accelerated molecular dynamics (TAMD)/driven-adiabatic free energy dynamics (d-AFED), unified free energy dynamics (UFED), and temperature accelerated sliced sampling (TASS), uses an extended variable formalism to achieve quick exploration of conformational space. These techniques are powerful, as they enhance the sampling of a large number of CVs simultaneously compared to other techniques. Extended variables are kept at a much higher temperature than the physical temperature by ensuring adiabatic separation between the extended and physical subsystems and employing rigorous thermostatting. In this work, we present a computational platform to perform extended phase space enhanced sampling simulations using the open-source molecular dynamics engine OpenMM. The implementation allows users to have interoperability of sampling techniques, as well as employ state-of-the-art thermostats and multiple time-stepping. This work also presents protocols for determining the critical parameters and procedures for reconstructing high-dimensional free energy surfaces. As a demonstration, we present simulation results on the high dimensional conformational landscapes of the alanine tripeptide in vacuo, tetra-N-methylglycine (tetra-sarcosine) peptoid in implicit solvent, and the Trp-cage mini protein in explicit water.

The "Beacon" Structural Model of Protein Folding: Application for Trp-Cage in Water

Protein folding is a process in which a polypeptide must undergo folding process to obtain its three-dimensional structure. Thermodynamically, it is a process of enthalpy to overcome the loss of conformational entropy in folding. Folding is primarily related to hydrophobic interactions and intramolecular hydrogen bondings. During folding, hydrophobic interactions are regarded to be the driving forces, especially in the initial structural collapse of a protein. Additionally, folding is guided by the strong interactions within proteins, such as intramolecular hydrogen bondings related to the α-helices and β-sheets of proteins. Therefore, a protein is divided into the folding key (FK) regions related to intramolecular hydrogen bondings and the non-folding key (non-FK) regions. Various conformations are expected for FK and non-FK regions. Different from non-FK regions, it is necessary for FK regions to form the specific conformations in folding, which are regarded as the necessary folding pathways (or "beacons"). Additionally, sequential folding is expected for the FK regions, and the intermediate state is found during folding. They are reflected on the local basins in the free energy landscape (FEL) of folding. To demonstrate the structural model, molecular dynamics (MD) simulations are conducted on the folding pathway of the TRP-cage in water.

Folding-unfolding asymmetry and a RetroFold computational algorithm

We treat protein folding as molecular self-assembly, while unfolding is viewed as disassembly. Fracture is typically a much faster process than self-assembly. Self-assembly is often an exponentially decaying process, since energy relaxes due to dissipation, while fracture is a constant-rate process as the driving force is opposed by damping. Protein folding takes two orders of magnitude longer than unfolding. We suggest a mathematical transformation of variables, which makes it possible to view self-assembly as time-reversed disassembly, thus folding can be studied as reversed unfolding. We investigate the molecular dynamics modelling of folding and unfolding of the short Trp-cage protein. Folding time constitutes about 800 ns, while unfolding (denaturation) takes only about 5.0 ns and, therefore, fewer computational resources are needed for its simulation. This RetroFold approach can be used for the design of a novel computation algorithm, which, while approximate, is less time-consuming than traditional folding algorithms.

Visualizing the Residue Interaction Landscape of Proteins by Temporal Network Embedding

Characterizing the structural dynamics of proteins with heterogeneous conformational landscapes is crucial to understanding complex biomolecular processes. To this end, dimensionality reduction algorithms are used to produce low-dimensional embeddings of the high-dimensional conformational phase space. However, identifying a compact and informative set of input features for the embedding remains an ongoing challenge. Here, we propose to harness the power of Residue Interaction Networks (RINs) and their centrality measures, established tools to provide a graph theoretical view on molecular structure. Specifically, we combine the closeness centrality, which captures global features of the protein conformation at residue-wise resolution, with EncoderMap, a hybrid neural-network autoencoder/multidimensional-scaling like dimensionality reduction algorithm. We find that the resulting low-dimensional embedding is a meaningful visualization of the residue interaction landscape that resolves structural details of the protein behavior while retaining global interpretability. This feature-based graph embedding of temporal protein graphs makes it possible to apply the general descriptive power of RIN formalisms to the analysis of protein simulations of complex processes such as protein folding and multidomain interactions requiring no protein-specific input. We demonstrate this on simulations of the fast folding protein Trp-Cage and the multidomain signaling protein FAT10. Due to its generality and modularity, the presented approach can easily be transferred to other protein systems.

The Role of Transient Intermediate Structures in the Unfolding of the Trp-Cage Fast-Folding Protein: Generating Ensembles from Time-Resolved X-ray Solution Scattering with Genetic Algorithms

The Trp-cage miniprotein is one of the smallest systems to exhibit a stable secondary structure and fast-folding dynamics, serving as an apt model system to study transient intermediates with both experimental and computational analyses. Previous spectroscopic characterizations that have been done on Trp-cage have inferred a single stable intermediate on a pathway from folded to unfolded basins. We aim to bridge the understanding of Trp-cage structural folding dynamics on microsecond-time scales, by utilizing time-resolved X-ray solution scattering to probe the temperature-induced unfolding pathway. Our results indicate the formation of a conformationally extended intermediate on the time scale of 1 μs, which undergoes complete unfolding within 5 μs. We further investigated the atomistic structural details of the unfolding pathway using a genetic algorithm to generate ensemble model fits to the scattering profiles. This analysis paves the way for direct benchmarking of theoretical models of protein folding ensembles produced with molecular dynamics simulations.

Topological bio-scaling analysis as a universal measure of protein folding

Scaling relationships for polymeric molecules establish power law dependencies between the number of molecular segments and linear dimensions, such as the radius of gyration. They also establish spatial topological properties of the chains, such as their dimensionality. In the spatial domain, power exponents α = 1 (linear stretched molecule), α = 0.5 (the ideal chain) and α = 0.333 (compact globule) are significant. During folding, the molecule undergoes the transition from the one-dimensional linear to the three-dimensional globular state within a very short time. However, intermediate states with fractional dimensions can be stabilized by modifying the solubility (e.g. by changing the solution temperature). Topological properties, such as dimension, correlate with the interaction energy, and thus by tuning the solubility one can control molecular interaction. We investigate these correlations using the example of a well-studied short model of Trp-cage protein. The radius of gyration is used to estimate the fractal dimension of the chain at different stages of folding. It is expected that the same principle is applicable to much larger molecules and that topological (dimensional) characteristics can provide insights into molecular folding and interactions.

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.

GraphVAMPNet, using graph neural networks and variational approach to Markov processes for dynamical modeling of biomolecules

Finding low dimensional representation of data from long-timescale trajectories of biomolecular processes such as protein-folding or ligand-receptor binding is of fundamental importance and kinetic models such as Markov modeling have proven useful in describing the kinetics of these systems. Recently, an unsupervised machine learning technique called VAMPNet was introduced to learn the low dimensional representation and linear dynamical model in an end-to-end manner. VAMPNet is based on variational approach to Markov processes (VAMP) and relies on neural networks to learn the coarse-grained dynamics. In this contribution, we combine VAMPNet and graph neural networks to generate an end-to-end framework to efficiently learn high-level dynamics and metastable states from the long-timescale molecular dynamics trajectories. This method bears the advantages of graph representation learning and uses graph message passing operations to generate an embedding for each datapoint which is used in the VAMPNet to generate a coarse-grained representation. This type of molecular representation results in a higher resolution and more interpretable Markov model than the standard VAMPNet enabling a more detailed kinetic study of the biomolecular processes. Our GraphVAMPNet approach is also enhanced with an attention mechanism to find the important residues for classification into different metastable states.

AIM: A Mapping Program for Infrared Spectroscopy of Proteins

Here, we present a new analysis program, AIM, that allows extracting the vibrational amide-I Hamiltonian using molecular dynamics trajectories for protein infrared spectroscopy modeling. The constructed Hamiltonians can be used as input for spectral calculations allowing the calculation of infrared absorption spectra, vibrational circular dichroism, and two-dimensional infrared spectra. These spectroscopies allow the study of the structure and dynamics of proteins. We will explain the essence of how AIM works and give examples of the information and spectra that can be obtained with the program using the Trypsin Inhibitor as an example. AIM is freely available from GitHub, and the package contains a demonstration allowing easy introduction to the use of the program.

Boosting the conformational sampling by combining replica exchange with solute tempering and well-sliced metadynamics

Methods that combine collective variable (CV) based enhanced sampling and global tempering approaches are used in speeding-up the conformational sampling and free energy calculation of large and soft systems with a plethora of energy minima. In this paper, a new method of this kind is proposed in which the well-sliced metadynamics approach (WSMTD) is united with replica exchange with solute tempering (REST2) method. WSMTD employs a divide-and-conquer strategy wherein high-dimensional slices of a free energy surface are independently sampled and combined. The method enables one to accomplish a controlled exploration of the CV-space with a restraining bias as in umbrella sampling, and enhance-sampling of one or more orthogonal CVs using a metadynamics like bias. The new hybrid method proposed here enables boosting the sampling of more slow degrees of freedom in WSMTD simulations, without the need to specify associated CVs, through a replica exchange scheme within the framework of REST2. The high-dimensional slices of the probability distributions of CVs computed from the united WSMTD and REST2 simulations are subsequently combined using the weighted histogram analysis method to obtain the free energy surface. We show that the new method proposed here is accurate, improves the conformational sampling, and achieves quick convergence in free energy estimates. We demonstrate this by computing the conformational free energy landscapes of solvated alanine tripeptide and Trp-cage mini protein in explicit water.

Variational embedding of protein folding simulations using Gaussian mixture variational autoencoders

Conformational sampling of biomolecules using molecular dynamics simulations often produces a large amount of high dimensional data that makes it difficult to interpret using conventional analysis techniques. Dimensionality reduction methods are thus required to extract useful and relevant information. Here, we devise a machine learning method, Gaussian mixture variational autoencoder (GMVAE), that can simultaneously perform dimensionality reduction and clustering of biomolecular conformations in an unsupervised way. We show that GMVAE can learn a reduced representation of the free energy landscape of protein folding with highly separated clusters that correspond to the metastable states during folding. Since GMVAE uses a mixture of Gaussians as its prior, it can directly acknowledge the multi-basin nature of the protein folding free energy landscape. To make the model end-to-end differentiable, we use a Gumbel-softmax distribution. We test the model on three long-timescale protein folding trajectories and show that GMVAE embedding resembles the folding funnel with folded states down the funnel and unfolded states outside the funnel path. Additionally, we show that the latent space of GMVAE can be used for kinetic analysis and Markov state models built on this embedding produce folding and unfolding timescales that are in close agreement with other rigorous dynamical embeddings such as time independent component analysis.

Resolving Dynamics in the Ensemble: Finding Paths through Intermediate States and Disordered Protein Structures

Proteins have been found to inhabit a diverse set of three-dimensional structures. The dynamics that govern protein interconversion between structures happen over a wide range of time scales─picoseconds to seconds. Our understanding of protein functions and dynamics is largely reliant upon our ability to elucidate physically populated structures. From an experimental structural characterization perspective, we are often limited to measuring the ensemble-averaged structure both in the steady-state and time-resolved regimes. Generating kinetic models and understanding protein structure-function relationships require atomistic knowledge of the populated states in the ensemble. In this Perspective, we present ensemble refinement methodologies that integrate time-resolved experimental signals with molecular dynamics models. We first discuss integration of experimental structural restraints to molecular models in disordered protein systems that adhere to the principle of maximum entropy for creating a complete set of ensemble structures. We then propose strategies to find kinetic pathways between the refined structures, using time-resolved inputs to guide molecular dynamics trajectories and the use of inference to generate tailored stimuli to prepare a desired ensemble of protein states.

CORE-MD II: A fast, adaptive, and accurate enhanced sampling method

In this paper, we present a fast and adaptive correlation guided enhanced sampling method (CORE-MD II). The CORE-MD II technique relies, in part, on partitioning of the entire pathway into short trajectories that we refer to as instances. The sampling within each instance is accelerated by adaptive path-dependent metadynamics simulations. The second part of this approach involves kinetic Monte Carlo (kMC) sampling between the different states that have been accessed during each instance. Through the combination of the partition of the total simulation into short non-equilibrium simulations and the kMC sampling, the CORE-MD II method is capable of sampling protein folding without any a priori definitions of reaction pathways and additional parameters. In the validation simulations, we applied the CORE-MD II on the dialanine peptide and the folding of two peptides: TrpCage and TrpZip2. In a comparison with long time equilibrium Molecular Dynamics (MD), 1 µs replica exchange MD (REMD), and CORE-MD I simulations, we find that the level of convergence of the CORE-MD II method is improved by a factor of 8.8, while the CORE-MD II method reaches acceleration factors of ∼120. In the CORE-MD II simulation of TrpZip2, we observe the formation of the native state in contrast to the REMD and the CORE-MD I simulations. The method is broadly applicable for MD simulations and is not restricted to simulations of protein folding or even biomolecules but also applicable to simulations of protein aggregation, protein signaling, or even materials science simulations.

Examination of Interactions of Hexafluoro-2-propanol with Trp-Cage in Hexafluoro-2-propanol-Water by MD Simulations and Intermolecular Nuclear Overhauser Effects

All-atom molecular dynamic simulations of the peptide Trp-cage in 30% hexafluoro-2-propanol- water (V/V) at 278 K have been carried out with the goal of exploring peptide hydrogen-solvent fluorine nuclear spin cross relaxation. Force field parameters for HFIP reported by Fioroni et al. along with the fluorine parameters of the TFE5 model reported by this lab were used. Water was represented by the TIP5P-Ew model. Peptide modeling used the AMBER99SB-ILDN force field. Translational diffusion coefficients of solution components at 278 K were predicted to within 35% of experimental values using these parameter sets. The simulations indicate that the solvent mixture is not homogeneous, with HFIP molecules clustered into aggregates as large as 53 fluoroalcohol molecules. The solvent environment of surface atoms of Trp-cage fluctuates between being HFIP-rich and more water-rich about every 10 ns. In accord with previous studies by other groups, the average concentration of HFIP near the surface of the peptide is significantly enhanced over the concentration of HFIP in the bulk solvent. In the simulations, ∼7% of the initial contacts between HFIP molecules and Trp-cage develop into peptide-fluoroalcohol interactions that persist for times as long as 8 ns. Most of the available experimental nuclear spin cross-relaxation rates (Σ_HF) for hydrogens of the Trp-cage in 30% HFIP-water are reproduced from the MD trajectories to within uncertainties of the experimental data and the simulations. However, a few calculated Σ_HF values for hydrogens of the Trp-cage do not agree with experiment. These tend to be situations where long-lived peptide-HFIP interactions are predicted. The disagreements between observed and calculated Σ_HF in these instances signal defects in the modeling parameters and procedures that are presently unrecognized.

Visible-light-driven photocontrol of the Trp-cage protein fold by a diazocine cross-linker

Diazocines are characterized by extraordinary photochemical properties rendering them of particular interest for switching the conformation of biomolecules with visible light. Current developments afford synthetic access to unprecedented diazocine derivatives promising particular opportunities in photocontrol of proteins and biological systems. In this work, the well-established approach of photocontrolling the secondary structure of α-helices was exploited using a diazocine to reversibly fold and unfold the tertiary structure of a small protein. The protein of choice was the globulary folded Trp-cage, a widely used model system for the elucidation of protein folding pathways. A specifically designed, short and rigid dicarboxy-functionalized diazocine-based cross-linker was attached to two solvent-exposed side chains at the α-helix of the miniprotein through the use of a primary amine-selective active ester. This cross-linking strategy is orthogonal to the common cysteine-based chemistry. The cross-linked Trp-cage was successfully photoisomerized and exhibited a strong correlation between protein fold and diazocine isomeric state. As determined by NMR spectroscopy, the cis-isomer stabilized the fold, while the trans-isomer led to complete protein unfolding. The successful switching of the protein fold in principle demonstrates the ability to control protein function, as the activity depends on their structural integrity.

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.

Refinement of Peptide Conformational Ensembles by 2D IR Spectroscopy: Application to Ala‒Ala‒Ala

Characterizing ensembles of intrinsically disordered proteins is experimentally challenging because of the ill-conditioned nature of ensemble determination with limited data and the intrinsic fast dynamics of the conformational ensemble. Amide I two-dimensional infrared (2D IR) spectroscopy has picosecond time resolution to freeze structural ensembles as needed for probing disordered-protein ensembles and conformational dynamics. Also, developments in amide I computational spectroscopy now allow a quantitative and direct prediction of amide I spectra based on conformational distributions drawn from molecular dynamics simulations, providing a route to ensemble refinement against experimental spectra. We performed a Bayesian ensemble refinement method on Ala-Ala-Ala against isotope-edited Fourier-transform infrared spectroscopy and 2D IR spectroscopy and tested potential factors affecting the quality of ensemble refinements. We found that isotope-edited 2D IR spectroscopy provides a stringent constraint on Ala-Ala-Ala conformations and returns consistent conformational ensembles with the dominant ppII conformer across varying prior distributions from many molecular dynamics force fields and water models. The dominant factor influencing ensemble refinements is the systematic frequency uncertainty from spectroscopic maps. However, the uncertainty of conformer populations can be significantly reduced by incorporating 2D IR spectra in addition to traditional Fourier-transform infrared spectra. Bayesian ensemble refinement against isotope-edited 2D IR spectroscopy thus provides a route to probe equilibrium-complex protein ensembles and potentially nonequilibrium conformational dynamics.

Optimizing the fold stability of the circularly permuted Trp-cage motif

Through optimization of the linker region and key stabilizing mutations, it has been possible to improve the stability of the circularly permuted (cp) Trp-cage miniprotein. However, even the most stable Trp-cage circular permutants are still less stable than the analogous standard topology (std) Trp-cages. Extending mutational studies of Trp-cage fold stability to cp-species, including analogs lacking chain terminal charges, has uncovered and quantitated some additional stabilizing and destabilizing interactions. Upon protonation, the circular permutants are destabilized to a much greater extent than the standard topology series. End effects, particularly Coulombic interactions, appear to be more important for the cp-series while the Y10/P4 interaction in the cp-series is not as significant a stabilizing feature as the corresponding Y3/P19 in the standard topology series.

Enhancing Biomolecular Sampling with Reinforcement Learning: A Tree Search Molecular Dynamics Simulation Method

This paper proposes a novel molecular simulation method, called tree search molecular dynamics (TS-MD), to accelerate the sampling of conformational transition pathways, which require considerable computation. In TS-MD, a tree search algorithm, called upper confidence bounds for trees, which is a type of reinforcement learning algorithm, is applied to sample the transition pathway. By learning from the results of the previous simulations, TS-MD efficiently searches conformational space and avoids being trapped in local stable structures. TS-MD exhibits better performance than parallel cascade selection molecular dynamics, which is one of the state-of-the-art methods, for the folding of miniproteins, Chignolin and Trp-cage, in explicit water.

A Hybrid Hamiltonian for the Accelerated Sampling along Experimental Restraints

In this article, we present an enhanced sampling method based on a hybrid Hamiltonian which combines experimental distance restraints with a bias dependent from multiple path-dependent variables. This simulation method determines the bias-coordinates on the fly and does not require a priori knowledge about reaction coordinates. The hybrid Hamiltonian accelerates the sampling of proteins, and, combined with experimental distance information, the technique considers the restraints adaptively and in dependency of the system's intrinsic dynamics. We validate the methodology on the dipole relaxation of two water models and the conformational landscape of dialanine. Using experimental NMR-restraint data, we explore the folding landscape of the TrpCage mini-protein and in a second example apply distance restraints from chemical crosslinking/mass spectrometry experiments for the sampling of the conformation space of the Killer Cell Lectin-like Receptor Subfamily B Member 1A (NKR-P1A). The new methodology has the potential to adaptively introduce experimental restraints without affecting the conformational space of the system along an ergodic trajectory. Since only a limited number of input- and no-order parameters are required for the setup of the simulation, the method is broadly applicable and has the potential to be combined with coarse-graining methods.

Enriched Conformational Sampling of DNA and Proteins with a Hybrid Hamiltonian Derived from the Protein Data Bank

In this article, we present a method for the enhanced molecular dynamics simulation of protein and DNA systems called potential of mean force (PMF)-enriched sampling. The method uses partitions derived from the potentials of mean force, which we determined from DNA and protein structures in the Protein Data Bank (PDB). We define a partition function from a set of PDB-derived PMFs, which efficiently compensates for the error introduced by the assumption of a homogeneous partition function from the PDB datasets. The bias based on the PDB-derived partitions is added in the form of a hybrid Hamiltonian using a renormalization method, which adds the PMF-enriched gradient to the system depending on a linear weighting factor and the underlying force field. We validated the method using simulations of dialanine, the folding of TrpCage, and the conformational sampling of the Dickerson⁻Drew DNA dodecamer. Our results show the potential for the PMF-enriched simulation technique to enrich the conformational space of biomolecules along their order parameters, while we also observe a considerable speed increase in the sampling by factors ranging from 13.1 to 82. The novel method can effectively be combined with enhanced sampling or coarse-graining methods to enrich conformational sampling with a partition derived from the PDB.

Molecular enhanced sampling with autoencoders: On-the-fly collective variable discovery and accelerated free energy landscape exploration

Macromolecular and biomolecular folding landscapes typically contain high free energy barriers that impede efficient sampling of configurational space by standard molecular dynamics simulation. Biased sampling can artificially drive the simulation along prespecified collective variables (CVs), but success depends critically on the availability of good CVs associated with the important collective dynamical motions. Nonlinear machine learning techniques can identify such CVs but typically do not furnish an explicit relationship with the atomic coordinates necessary to perform biased sampling. In this work, we employ auto-associative artificial neural networks ("autoencoders") to learn nonlinear CVs that are explicit and differentiable functions of the atomic coordinates. Our approach offers substantial speedups in exploration of configurational space, and is distinguished from existing approaches by its capacity to simultaneously discover and directly accelerate along data-driven CVs. We demonstrate the approach in simulations of alanine dipeptide and Trp-cage, and have developed an open-source and freely available implementation within OpenMM. © 2018 Wiley Periodicals, Inc.

Hydrophobic Collapse in N-Methylacetamide-Water Mixtures

Aqueous N-methylacetamide solutions were investigated by polarization-resolved pump-probe and 2D infrared spectroscopy (2D IR), using the amide I mode as a reporter. The 2D IR results are compared with molecular dynamics simulations and spectral calculations to gain insight into the molecular structures in the mixture. N-Methylacetamide and water molecules tend to form clusters with "frozen" amide I dynamics. This is driven by a hydrophobic collapse as the methyl groups of the N-methylacetamide molecules cluster in the presence of water. Since the studied system can be considered as a simplified model for the backbone of proteins, the present study forms a convenient basis for understanding the structural and vibrational dynamics in proteins. It is particularly interesting to find out that a hydrophobic collapse as the one driving protein folding is observed in such a simple system.

Multiscale Simulation of Protein Hydration Using the SWINGER Dynamical Clustering Algorithm

To perform computationally efficient concurrent multiscale simulations of biological macromolecules in solution, where the all-atom (AT) models are coupled to supramolecular coarse-grained (SCG) solvent models, previous studies resorted to modified AT water models, such as the bundled-simple point charge (SPC) models, that use semiharmonic springs to restrict the relative movement of water molecules within a cluster. Those models can have a significant impact on the simulated biomolecules and can lead, for example, to a partial unfolding of a protein. In this work, we employ the recently developed alternative approach with a dynamical clustering algorithm, SWINGER, which enables a direct coupling of original unmodified AT and SCG water models. We perform an adaptive resolution molecular dynamics simulation of a Trp-Cage miniprotein in multiscale water, where the standard SPC water model is interfaced with the widely used MARTINI SCG model, and demonstrate that, compared to the corresponding full-blown AT simulations, the structural and dynamic properties of the solvated protein and surrounding solvent are well reproduced by our approach.

An adaptive bias - hybrid MD/kMC algorithm for protein folding and aggregation

In this paper, we present a novel hybrid Molecular Dynamics/kinetic Monte Carlo (MD/kMC) algorithm and apply it to protein folding and aggregation in explicit solvent. The new algorithm uses a dynamical definition of biases throughout the MD component of the simulation, normalized in relation to the unbiased forces. The algorithm guarantees sampling of the underlying ensemble in dependency of one average linear coupling factor 〈α〉_τ. We test the validity of the kinetics in simulations of dialanine and compare dihedral transition kinetics with long-time MD-simulations. We find that for low 〈α〉_τ values, kinetics are in good quantitative agreement. In folding simulations of TrpCage and TrpZip4 in explicit solvent, we also find good quantitative agreement with experimental results and prior MD/kMC simulations. Finally, we apply our algorithm to study growth of the Alzheimer Amyloid Aβ 16-22 fibril by monomer addition. We observe two possible binding modes, one at the extremity of the fibril (elongation) and one on the surface of the fibril (lateral growth), on timescales ranging from ns to 8 μs.

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.

Genetically encoded photocross-linkers determine the biological binding site of exendin-4 peptide in the N-terminal domain of the intact human glucagon-like peptide-1 receptor (GLP-1R)

The glucagon-like peptide-1 receptor (GLP-1R) is a key therapeutic target in the management of type II diabetes mellitus, with actions including regulation of insulin biosynthesis and secretion, promotion of satiety, and preservation of β-cell mass. Like most class B G protein-coupled receptors (GPCRs), there is limited knowledge linking biological activity of the GLP-1R with the molecular structure of an intact, full-length, and functional receptor·ligand complex. In this study, we have utilized genetic code expansion to site-specifically incorporate the photoactive amino acid p-azido-l-phenylalanine (azF) into N-terminal residues of a full-length functional human GLP-1R in mammalian cells. UV-mediated photolysis of azF was then carried out to induce targeted photocross-linking to determine the proximity of the azido group in the mutant receptor with the peptide exendin-4. Cross-linking data were compared directly with the crystal structure of the isolated N-terminal extracellular domain of the GLP-1R in complex with exendin(9-39), revealing both similarities as well as distinct differences in the mode of interaction. Generation of a molecular model to accommodate the photocross-linking constraints highlights the potential influence of environmental conditions on the conformation of the receptor·peptide complex, including folding dynamics of the peptide and formation of dimeric and higher order oligomeric receptor multimers. These data demonstrate that crystal structures of isolated receptor regions may not give a complete reflection of peptide/receptor interactions and should be combined with additional experimental constraints to reveal peptide/receptor interactions occurring in the dynamic, native, and full-length receptor state.

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.

Itinerary profiling to analyze a large number of protein-folding trajectories

Understanding how proteins fold through a vast number of unfolded states is a major subject in the study of protein folding. Herein, we present itinerary profiling as a simple method to analyze molecular dynamics trajectories, and apply this method to Trp-cage. In itinerary profiling, structural clusters included in a trajectory are represented by a bit sequence, and a number of trajectories, as well as the structural clusters, can be compared and classified. As a consequence, the structural clusters that characterize the foldability of trajectories were able to be identified. The connections between the clusters were then illustrated as a network and the structural features of the clusters were examined. We found that in the true folding funnel, Trp-cage formed a left-handed main-chain topology and the Trp6 side-chain was located at the front of the main-chain ring, even in the initial unfolded states. In contrast, in the false folding funnel of the pseudo-native states, in which the Trp6 side-chain is upside down in the protein core, Trp-cage had a right-handed main-chain topology and the Trp side-chain was at the back. The initial topological partition, as determined by the main-chain handedness and the location of the Trp residue, predetermines Trp-cage foldability and the destination of the trajectory to the native state or the pseudo-native states.

Folding of a Zinc-Finger ββα-Motif Investigated Using Two-Dimensional and Time-Resolved Vibrational Spectroscopy

Small proteins provide good model systems for studying the fundamental forces that control protein folding. Here, we investigate the folding dynamics of the 28-residue zinc-finger mutant FSD-1, which is designed to form a metal-independent folded ββα-motif, and which provides a testing ground for proteins containing a mixed α/β fold. Although the folding of FSD-1 has been actively studied, the folding mechanism remains largely unclear. In particular, it is unclear in what stage of folding the α-helix is formed. To address this issue we investigate the folding mechanism of FSD-1 using a combination of temperature-dependent UV circular dichroism (UV-CD), Fourier transform infrared (FTIR) spectroscopy, two-dimensional infrared (2D-IR) spectroscopy, and temperature-jump (T-jump) transient-IR spectroscopy. Our UV-CD and FTIR data show different thermal melting transitions, indicating multistate folding behavior. Temperature-dependent 2D-IR spectra indicate that the α-helix is the most stable structural element of FSD-1. To investigate the folding/unfolding re-equilibration dynamics of FSD-1, the conformational changes induced by a nanosecond T-jump are probed with transient-IR and transient dispersed-pump-probe (DPP) IR spectroscopy. We observe biexponential T-jump relaxation kinetics (with time constants of 80 ± 13 ns and 1300 ± 100 ns at 322 K), confirming that the folding involves an intermediate state. The IR and dispersed-pump-probe IR spectra associated with the two kinetic components suggest that the folding of FSD-1 involves early formation of the α-helix, followed by the formation of the β-hairpin and hydrophobic contacts.

Application of two-dimensional infrared spectroscopy to benchmark models for the amide I band of proteins

In this paper, we present a novel benchmarking method for validating the modelling of vibrational spectra for the amide I region of proteins. We use the linear absorption spectra and two-dimensional infrared spectra of four experimentally well-studied proteins as a reference and test nine combinations of molecular dynamics force fields, vibrational frequency mappings, and coupling models. We find that two-dimensional infrared spectra provide a much stronger test of the models than linear absorption does. The best modelling approach in the present study still leaves significant room for future improvement. The presented benchmarking scheme, thus, provides a way of validating future protocols for modelling the amide I band in proteins.

Coarse kMC-based replica exchange algorithms for the accelerated simulation of protein folding in explicit solvent

In this paper, we present a coarse replica exchange molecular dynamics (REMD) approach, based on kinetic Monte Carlo (kMC).

A hydrodynamic view of the first-passage folding of Trp-cage miniprotein

We study folding of Trp-cage miniprotein in the conditions when the native state of the protein is stable and unfolding events are improbable, which corresponds to physiological conditions. Using molecular dynamics simulations with an implicit solvent model, an ensemble of folding trajectories from unfolded (practically extended) states of the protein to the native state was generated. To get insight into the folding kinetics, the free energy surface and kinetic network projected on this surface were constructed. This, "conventional" analysis of the folding reaction was followed by a recently proposed hydrodynamic description of protein folding (Chekmarev et al. in Phys Rev Lett 100(1):018107, 2008), in which the process of the first-passage folding is viewed as a stationary flow of a folding "fluid" from the unfolded to native state. This approach is conceptually different from the previously used approaches and thus allows an alternative view of the folding dynamics and kinetics of Trp-cage, the conclusions about which are very diverse. In agreement with most previous studies, we observed two characteristic folding pathways: in one pathway (I), the collapse of the hydrophobic core precedes the formation of the [Formula: see text]-helix, and in the other pathway (II), these events occur in the reverse order. We found that although pathway II is complicated by a repeated partial protein unfolding, it contributes to the total folding flow as little as ≈10%, so that the folding kinetics remain essentially single-exponential.

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

To test whether our recently developed residue-specific force field RSFF2 can reproduce the mutational effect on the thermal stability of Trp-cage mini-protein and decipher its detailed folding mechanism, we carried out long-time replica-exchange molecular dynamics (REMD) simulations on five Trp-cage variants, including TC5b and TC10b. Initiated from their unfolded structures, the simulations not only well-reproduce their experimental structures but also their melting temperatures and folding enthalpies reasonably well. For each Trp-cage variant, the overall folding free energy landscape is apparently two-state, but some intermediate states can be observed when projected on more detailed coordinates. We also found different variants have the same major folding pathway, including the well formed PII-helix in the unfolded state, the formation of W6-P12/P18/P19 contacts and the α-helix before the transition state, the following formation of most native contacts, and the final native loop formation. The folding mechanism derived here is consistent with many previous simulations and experiments.

Influence of the ionic liquid [C4mpy][Tf2N] on the structure of the miniprotein Trp-cage

We examine the effect of the ionic liquid [C4mpy][Tf2N] on the structure of the miniprotein Trp-cage and contrast these results with the behavior of Trp-cage in water. We find the ionic liquid has a dramatic effect on Trp-cage, though many similarities with aqueous Trp-cage are observed. We assess Trp-cage folding by monitoring root mean square deviation from the crystallographic structure, radius of gyration, proline cis/trans isomerization state, protein secondary structure, amino acid contact formation and distance, and native and non-native contact formation. Starting from an unfolded configuration, Trp-cage folds in water at 298 K in less than 500 ns of simulation, but has very little mobility in the ionic liquid at the same temperature, which can be ascribed to the higher ionic liquid viscosity. At 365 K, the mobility of the ionic liquid is increased and initial stages of Trp-cage folding are observed, however Trp-cage does not reach the native folded state in 2 μs of simulation in the ionic liquid. Therefore, in addition to conventional molecular dynamics, we also employ scaled molecular dynamics to expedite sampling, and we demonstrate that Trp-cage in the ionic liquid does closely approach the aqueous folded state. Interestingly, while the reduced mobility of the ionic liquid is found to restrict Trp-cage motion, the ionic liquid does facilitate proline cis/trans isomerization events that are not seen in our aqueous simulations.

Parallel continuous simulated tempering and its applications in large-scale molecular simulations

In this paper, we introduce a parallel continuous simulated tempering (PCST) method for enhanced sampling in studying large complex systems. It mainly inherits the continuous simulated tempering (CST) method in our previous studies [C. Zhang and J. Ma, J. Chem. Phys. 130, 194112 (2009); C. Zhang and J. Ma, J. Chem. Phys. 132, 244101 (2010)], while adopts the spirit of parallel tempering (PT), or replica exchange method, by employing multiple copies with different temperature distributions. Differing from conventional PT methods, despite the large stride of total temperature range, the PCST method requires very few copies of simulations, typically 2-3 copies, yet it is still capable of maintaining a high rate of exchange between neighboring copies. Furthermore, in PCST method, the size of the system does not dramatically affect the number of copy needed because the exchange rate is independent of total potential energy, thus providing an enormous advantage over conventional PT methods in studying very large systems. The sampling efficiency of PCST was tested in two-dimensional Ising model, Lennard-Jones liquid and all-atom folding simulation of a small globular protein trp-cage in explicit solvent. The results demonstrate that the PCST method significantly improves sampling efficiency compared with other methods and it is particularly effective in simulating systems with long relaxation time or correlation time. We expect the PCST method to be a good alternative to parallel tempering methods in simulating large systems such as phase transition and dynamics of macromolecules in explicit solvent.

Correlation analysis for heat denaturation of Trp-cage miniprotein with explicit solvent

Energetics was analyzed for Trp-cage miniprotein in water to elucidate the solvation effect in heat denaturation. The solvation free energy was computed for a set of protein structures at room and high temperatures with all-atom molecular dynamics simulation combined with the solution theory in the energy representation, and its correlations were investigated against the intramolecular (structural) energy of the protein and the average interaction energy of the protein with the solvent water. It was observed both at room and high temperatures that the solvation free energy is anticorrelated to the structural energy and varies in parallel to the electrostatic component of the protein-water interaction energy without correlations to the van der Waals and excluded-volume components. When the set of folded structures sampled at room temperature was compared with the set of unfolded ones at high temperature, it was found that the preference order of the two sets is in correspondence to the van der Waals and excluded-volume components in the sum of the protein intramolecular and protein-water intermolecular interactions and is not distinguished by the electrostatic component.

Accelerated molecular dynamics simulations of protein folding

Folding of four fast-folding proteins, including chignolin, Trp-cage, villin headpiece and WW domain, was simulated via accelerated molecular dynamics (aMD). In comparison with hundred-of-microsecond timescale conventional molecular dynamics (cMD) simulations performed on the Anton supercomputer, aMD captured complete folding of the four proteins in significantly shorter simulation time. The folded protein conformations were found within 0.2-2.1 Å of the native NMR or X-ray crystal structures. Free energy profiles calculated through improved reweighting of the aMD simulations using cumulant expansion to the second-order are in good agreement with those obtained from cMD simulations. This allows us to identify distinct conformational states (e.g., unfolded and intermediate) other than the native structure and the protein folding energy barriers. Detailed analysis of protein secondary structures and local key residue interactions provided important insights into the protein folding pathways. Furthermore, the selections of force fields and aMD simulation parameters are discussed in detail. Our work shows usefulness and accuracy of aMD in studying protein folding, providing basic references in using aMD in future protein-folding studies.

Tuning the Attempt Frequency of Protein Folding Dynamics via Transition-State Rigidification: Application to Trp-Cage

The attempt frequency or prefactor (k0) of the transition-state rate equation of protein folding kinetics has been estimated to be on the order of 10(6) s(-1), which is many orders of magnitude smaller than that of chemical reactions. Herein we use the mini-protein Trp-cage to show that it is possible to significantly increase the value of k0 for a protein folding reaction by rigidifying the transition state. This is achieved by reducing the conformational flexibility of a key structural element (i.e., an α-helix) formed in the transition state via photoisomerization of an azobenzene cross-linker. We find that this strategy not only decreases the folding time of the Trp-cage peptide by more than an order of magnitude (to ∼100 ns at 25°C) but also exposes parallel folding pathways, allowing us to provide, to the best of our knowledge, the first quantitative assessment of the curvature of the transition-state free-energy surface of a protein.