Folding Dynamics of 10-Residue β-Hairpin Peptide Chignolin

Folding Kinetics and Volume Variation of the β-Hairpin Peptide Chignolin upon Ultrafast pH-Jumps

In-depth characterization of fundamental folding steps of small model peptides is crucial for a better understanding of the folding mechanisms of more complex biomacromolecules. We have previously reported on the folding/unfolding kinetics of a model α-helix. Here, we study folding transitions in chignolin (GYDPETGTWG), a short β-hairpin peptide previously used as a model to study conformational changes in β-sheet proteins. Although previously suggested, until now, the role of the Tyr2-Trp9 interaction in the folding mechanism of chignolin was not clear. In the present work, pH-dependent conformational changes of chignolin were characterized by circular dichroism (CD), nuclear magnetic resonance (NMR), ultrafast pH-jump coupled with time-resolved photoacoustic calorimetry (TR-PAC), and molecular dynamics (MD) simulations. Taken together, our results present a comprehensive view of chignolin's folding kinetics upon local pH changes and the role of the Tyr2-Trp9 interaction in the folding process. CD data show that chignolin's β-hairpin formation displays a pH-dependent skew bell-shaped curve, with a maximum close to pH 6, and a large decrease in β-sheet content at alkaline pH. The β-hairpin structure is mainly stabilized by aromatic interactions between Tyr2 and Trp9 and CH-π interactions between Tyr2 and Pro4. Unfolding of chignolin at high pH demonstrates that protonation of Tyr2 is essential for the stability of the β-hairpin. Refolding studies were triggered by laser-induced pH-jumps and detected by TR-PAC. The refolding of chignolin from high pH, mainly due to the protonation of Tyr2, is characterized by a volume expansion (10.4 mL mol-1), independent of peptide concentration, in the microsecond time range (lifetime of 1.15 μs). At high pH, the presence of the deprotonated hydroxyl (tyrosinate) hinders the formation of the aromatic interaction between Tyr2 and Trp9 resulting in a more disorganized and dynamic tridimensional structure of the peptide. This was also confirmed by comparing MD simulations of chignolin under conditions mimicking neutral and high pH.

Flexible Gaussian Accelerated Molecular Dynamics to Enhance Biological Sampling

Molecular dynamics simulations often struggle to obtain sufficient sampling to study complex molecular events due to high energy barriers separating the minima of interest. Multiple enhanced sampling techniques have been developed and improved over the years to tackle this issue. Gaussian accelerated molecular dynamics (GaMD) is a recently developed enhanced sampling technique that works by adding a biasing potential, lifting the energy landscape up, and decreasing the height of its barriers. GaMD allows one to increase the sampling of events of interest without the need of a priori knowledge of the system or the relevant coordinates. All required acceleration parameters can be obtained from a previous search run. Upon its development, several improvements for the methodology have been proposed, among them selective GaMD in which the boosting potential is selectively applied to the region of interest. There are currently four selective GaMD methods that have shown promising results. However, all of these methods are constrained on the number, location, and scenarios in which this selective boosting potential can be applied to ligands, peptides, or protein-protein interactions. In this work, we showcase a GROMOS implementation of the GaMD methodology with a fully flexible selective GaMD approach that allows the user to define, in a straightforward way, multiple boosting potentials for as many regions as desired. We show and analyze the advantages of this flexible selective approach on two previously used test systems, the alanine dipeptide and the chignolin peptide, and extend these examples to study its applicability and potential to study conformational changes of glycans and glycosylated proteins.

A study of ab initio folding of chignolins using replica-exchange molecular dynamics simulations

More and more studies have confirmed the importance of polarization effects in hydrogen bonding interactions in protein folding simulations. In this paper, a recently developed charge update scheme termed polarized structure-specific backbone charge (PSBC) model was applied to the folding of 10-residue chignolin. A comparison between simulations performed using PSBC and a nonpolarizable (AMBER99SB) force field demonstrably showed the importance of the electrostatic polarization effect in the folding of the short β-hairpin peptide by a series of analyses such as DSSP, free-energy landscape, hydrogen bond occupancy, and melting curve. The PSBC model was further validated by folding two other chignolin variants.

Granger Causality Analysis of Chignolin Folding

Constantly advancing computer simulations of biomolecules provide huge amounts of data that are difficult to interpret. In particular, obtaining insights into functional aspects of macromolecular dynamics, often related to cascades of transient events, calls for methodologies that depart from the well-grounded framework of equilibrium statistical physics. One of the approaches toward the analysis of complex temporal data which has found applications in the fields of neuroscience and econometrics is Granger causality analysis. It allows determining which components of multidimensional time series are most influential for the evolution of the entire system, thus providing insights into causal relations within the dynamic structure of interest. In this work, we apply Granger analysis to a long molecular dynamics trajectory depicting repetitive folding and unfolding of a mini β-hairpin protein, CLN025. We find objective, quantitative evidence indicating that rearrangements within the hairpin turn region are determinant for protein folding and unfolding. On the contrary, interactions between hairpin arms score low on the causality scale. Taken together, these findings clearly favor the concept of zipperlike folding, which is one of two postulated β-hairpin folding mechanisms. More importantly, the results demonstrate the possibility of a conclusive application of Granger causality analysis to a biomolecular system.

Effects of Externally Applied Electric Fields on the Manipulation of Solvated-Chignolin Folding: Static- versus Alternating-Field Dichotomy at Play

The interaction between a protein and external electric field (EF) can alter its structure and dynamical behavior, which has a potential impact on the biological function of proteins and cause uncertain health consequences. Conversely, the application of EFs of judiciously selected intensity and frequency can help to treat disease, and optimization of this requires a greater understanding of EF-induced effects underpinning basic protein biophysics. In the present study, chignolin─an artificial protein sufficiently small to undergo fast-folding events and transitions─was selected as an ideal prototype to investigate how, and to what extent, externally applied electric fields may manipulate or influence protein-folding phenomena. Nonequilibrium molecular dynamics (NEMD) simulations have been performed of solvated chignolin to determine the distribution of folding states and their underlying transition dynamics, in the absence and presence of externally applied electric fields (both static and alternating); a key focus has been to ascertain how folding pathways are altered in an athermal sense by external fields. Compared to zero-field conditions, a dramatically different─indeed, bifurcated─behavior of chignolin-folding processes emerges between static- and alternating-field scenarios, especially vis-à-vis incipient stages of hydrophobic-core formation: in alternating fields, fold-state populations diversified, with an attendant acceleration of state-hopping folding kinetics, featuring the concomitant emergence of a new, quasi-stable structure compared to the native structure, in field-shifted energy landscapes.

Manipulating turn residues on de novo designed β-hairpin peptides for selectivity against drug-resistant bacteria

Synthetic β-hairpin antimicrobial peptides (AMPs) offer a useful source for the development of novel antimicrobial agents. β-hairpin peptides generally consist of two side strands bridged by a reverse turn. In literature, most studies focused on the modifications of the side strands to manipulate the stability and activity of β-hairpin peptides, and much less is known about the impact of the turn region. By designing a series of de novo β-hairpin peptides with identical side strands but varied turns, we demonstrated that mutations of only 2 to 4 amino acids at the turn region could impart a wide range of antimicrobial profiles among synthetic β-hairpin AMPs. BTT2-4 and BTT6 displayed selective potency against Gram-negative bacteria, with minimum inhibitory concentrations (MICs) of 4-8 µM. In contrast, BTT1 exhibited broad-spectrum activity, with MICs of 4-8 µM against both Gram-positive and Gram-negative strains. Additionally, BTT1 was potent against methicillin-resistant Staphylococcus aureus (MRSA) and colistin-resistant Enterobacterales. The antimicrobial potency of BTT1 persisted after 14 days of serial passage. Mechanistic studies revealed that interactions between lipopolysaccharide (LPS) and the peptides were critical to their membranolytic activity against the bacterial inner membrane. Aside from folding stability, we observed that a degree of conformational flexibility was required for disruptive membrane interactions. STATEMENT OF SIGNIFICANCE: By examining the significance of the turn region of β-hairpin peptides, we present valuable knowledge to the design toolkit of novel antimicrobial peptides as alternative therapeutics to overcome antibiotic resistance. Our de novo designed synthetic peptides displayed selective activity against Gram-negative bacteria and potent activity against clinically relevant antibiotic-resistant strains (e.g. colistin-resistant Enterobacterales and methicillin-resistant Staphylococcus aureus). The bactericidal activity of our peptides was shown to be robust in the presence of proteolytic trypsin and saline, conditions that could suppress peptide activity. Our peptides were also determined to be non-cytotoxic against a human cell line.

Machine learning implicit solvation for molecular dynamics

Accurate modeling of the solvent environment for biological molecules is crucial for computational biology and drug design. A popular approach to achieve long simulation time scales for large system sizes is to incorporate the effect of the solvent in a mean-field fashion with implicit solvent models. However, a challenge with existing implicit solvent models is that they often lack accuracy or certain physical properties compared to explicit solvent models as the many-body effects of the neglected solvent molecules are difficult to model as a mean field. Here, we leverage machine learning (ML) and multi-scale coarse graining (CG) in order to learn implicit solvent models that can approximate the energetic and thermodynamic properties of a given explicit solvent model with arbitrary accuracy, given enough training data. Following the previous ML-CG models CGnet and CGSchnet, we introduce ISSNet, a graph neural network, to model the implicit solvent potential of mean force. ISSNet can learn from explicit solvent simulation data and be readily applied to molecular dynamics simulations. We compare the solute conformational distributions under different solvation treatments for two peptide systems. The results indicate that ISSNet models can outperform widely used generalized Born and surface area models in reproducing the thermodynamics of small protein systems with respect to explicit solvent. The success of this novel method demonstrates the potential benefit of applying machine learning methods in accurate modeling of solvent effects for in silico research and biomedical applications.

Overcoming Free-Energy Barriers with a Seamless Combination of a Biasing Force and a Collective Variable-Independent Boost Potential

Amid collective-variable (CV)-based importance-sampling algorithms, a hybrid of the extended adaptive biasing force and the well-tempered metadynamics algorithms (WTM-eABF) has proven particularly cost-effective for exploring the rugged free-energy landscapes that underlie biological processes. However, as an inherently CV-based algorithm, this hybrid scheme does not explicitly accelerate sampling in the space orthogonal to the chosen CVs, thereby limiting its efficiency and accuracy, most notably in those cases where the slow degrees of freedom of the process at hand are not accounted for in the model transition coordinate. Here, inspired by Gaussian-accelerated molecular dynamics (GaMD), we introduce the same CV-independent harmonic boost potential into WTM-eABF, yielding a hybrid algorithm coined GaWTM-eABF. This algorithm leans on WTM-eABF to explore the transition coordinate with a GaMD-mollified potential and recovers the unbiased free-energy landscape through thermodynamic integration followed by proper reweighting. As illustrated in our numerical tests, GaWTM-eABF effectively overcomes the free-energy barriers in orthogonal space and correctly recovers the unbiased potential of mean force (PMF). Furthermore, applying both GaWTM-eABF and WTM-eABF to two biologically relevant processes, namely, the reversible folding of (i) deca-alanine and (ii) chignolin, our results indicate that GaWTM-eABF reduces the uncertainty in the PMF calculation and converges appreciably faster than WTM-eABF. Obviating the need of multiple-copy strategies, GaWTM-eABF is a robust, computationally efficient algorithm to surmount the free-energy barriers in orthogonal space and maps with utmost fidelity the free-energy landscape along selections of CVs. Moreover, our strategy that combines WTM-eABF with GaMD can be easily extended to other biasing-force algorithms.

Chirality Driven Twisting as a Driving Force of Primitive Folding in Binary Mixtures

The N-trifluoroacetylated α-aminoalcohols (TFAAAs) are able to form quasi-one-dimensional supramolecular fibers (strings) when chirally pure, and isometric precipitates in the racemate. The strings' formation leads to the reversible gelation of the solution. The fresh gels occupy all the available volume, however during the incubation, they contract and concentrate in the central region of the tube. The microscopic observations revealed the growth of the strings' diameter and their rotation in the course of the incubation at the hour time-scale. The rotation provides for the hairpins forming that serve as hooks on the rotating string, which provides for coiling of the strings, which was observed as gel contraction. The morphology of the twisted strings resembles the structures observed in modern proteins, which allows drawing an analogy between the folding of biopolymers and the formation of the clew of strings. In addition, the rotation found in the TFAAA gels is an example of a simple system converting the energy of intermolecular agglutination to the rotational movement, so they could be considered as molecular motors.

Mutation-induced change in chignolin stability from π-turn to α-turn

A mutation from threonine to proline at the eighth residue in chignolin changes π-turn to α-turn.

Analysis of molecular dynamics simulations of 10-residue peptide, chignolin, using statistical mechanics: Relaxation mode analysis and three-dimensional reference interaction site model theory

Molecular dynamics simulation is a fruitful tool for investigating the structural stability, dynamics, and functions of biopolymers at an atomic level. In recent years, simulations can be performed on time scales of the order of milliseconds using special purpose systems. Since the most stable structure, as well as meta-stable structures and intermediate structures, is included in trajectories in long simulations, it is necessary to develop analysis methods for extracting them from trajectories of simulations. For these structures, methods for evaluating the stabilities, including the solvent effect, are also needed. We have developed relaxation mode analysis to investigate dynamics and kinetics of simulations based on statistical mechanics. We have also applied the three-dimensional reference interaction site model theory to investigate stabilities with solvent effects. In this paper, we review the results for designing amino-acid substitution of the 10-residue peptide, chignolin, to stabilize the misfolded structure using these developed analysis methods.

Replica sub-permutation method for molecular dynamics and monte carlo simulations

We propose an improvement of the replica-exchange and replica-permutation methods, which we call the replica sub-permutation method (RSPM). Instead of considering all permutations, this method uses a new algorithm referred to as sub-permutation to perform parameter transition. The RSPM succeeds in reducing the number of combinations between replicas and parameters without the loss of sampling efficiency. For comparison, we applied the replica sub-permutation, replica-permutation, and replica-exchange methods to a β-hairpin mini protein, chignolin, in explicit water. We calculated the transition ratio and number of tunneling events in the parameter space, the number of folding-unfolding events, the autocorrelation function, and the autocorrelation time as measures of sampling efficiency. The results indicate that among the three methods, the proposed RSPM is the most efficient in both parameter and conformational spaces. © 2019 Wiley Periodicals, Inc.

Molecular simulations by generalized-ensemble algorithms in isothermal-isobaric ensemble

Generalized-ensemble algorithms are powerful techniques for investigating biomolecules such as protein, DNA, lipid membrane, and glycan. The generalized-ensemble algorithms were originally developed in the canonical ensemble. On the other hand, not only temperature but also pressure is controlled in experiments. Additionally, pressure is used as perturbation to study relationship between function and structure of biomolecules. For this reason, it is important to perform efficient conformation sampling based on the isothermal-isobaric ensemble. In this article, we review a series of the generalized-ensemble algorithms in the isothermal-isobaric ensemble: multibaric-multithermal, pressure- and temperature-simulated tempering, replica-exchange, and replica-permutation methods. These methods achieve more efficient simulation than the conventional isothermal-isobaric simulation. Furthermore, the isothermal-isobaric generalized-ensemble simulation samples conformations of biomolecules from wider range of temperature and pressure. Thus, we can estimate physical quantities more accurately at any temperature and pressure values. The applications to the biomolecular system are also presented.

Conformational change of a biomolecule studied by the weighted ensemble method: Use of the diffusion map method to extract reaction coordinates

We simulate the nonequilibrium ensemble dynamics of a biomolecule using the weighted ensemble method, which was introduced in molecular dynamics simulations by Huber and Kim and further developed by Zuckerman and co-workers. As the order parameters to characterize its conformational change, we here use the coordinates derived from the diffusion map (DM) method, one of the manifold learning techniques. As a concrete example, we study the kinetic properties of a small peptide, chignolin in explicit water, and calculate the conformational change between the folded and misfolded states in a nonequilibrium way. We find that the transition time scales thus obtained are comparable to those using previously employed hydrogen-bond distances as the order parameters. Since the DM method only uses the 3D Cartesian coordinates of a peptide, this shows that the DM method can extract the important distance information of the peptide without relying on chemical intuition. The time scales are compared well with the previous results using different techniques, non-Markovian analysis and core-set milestoning for a single long trajectory. We also find that the most significant DM coordinate turns out to extract a dihedral angle of glycine, and the previously studied relaxation modes are well correlated with the most significant DM coordinates.

Analysis of Structural Stability of Chignolin

We discuss the stability of an entire protein and the influence of main chains and side chains of individual amino acids to investigate the protein-folding mechanism. For this purpose, we calculated the solvation free-energy contribution of individual atoms using the three-dimensional reference interaction site model with the atomic decomposition method. We generated structures of chignolin miniprotein by a molecular dynamics simulation and classified them into six types: native 1, native 2, misfolded 1, misfolded 2, intermediate, and unfolded states. The total energies of the native (-171.1 kcal/mol) and misfolded (-171.2 kcal/mol) states were almost the same and lower than those of the intermediate (-158.5 kcal/mol) and unfolded (-148.1 kcal/mol) states; however, their components were different. In the native state, the side-chain interaction between Thr6 and Thr8 is important for the formation of π-turn. On the other hand, the hydrogen bonds between the atoms of the main chains in the misfolded state become stronger than those in the intermediate state.

Relaxation mode analysis for molecular dynamics simulations of proteins

Molecular dynamics simulation is a powerful method for investigating the structural stability, dynamics, and function of biopolymers at the atomic level. In recent years, it has become possible to perform simulations on time scales of the order of milliseconds using special hardware. However, it is necessary to derive the important factors contributing to structural change or function from the complicated movements of biopolymers obtained from long simulations. Although some analysis methods for protein systems have been developed using increasing simulation times, many of these methods are static in nature (i.e., no information on time). In recent years, dynamic analysis methods have been developed, such as the Markov state model and relaxation mode analysis (RMA), which was introduced based on spin and homopolymer systems. The RMA method approximately extracts slow relaxation modes and rates from trajectories and decomposes the structural fluctuations into slow relaxation modes, which characterize the slow relaxation dynamics of the system. Recently, this method has been applied to biomolecular systems. In this article, we review RMA and its improved versions for protein systems.

Modeling the mechanism of CLN025 beta-hairpin formation

Beta-hairpins are substructures found in proteins that can lend insight into more complex systems. Furthermore, the folding of beta-hairpins is a valuable test case for benchmarking experimental and theoretical methods. Here, we simulate the folding of CLN025, a miniprotein with a beta-hairpin structure, at its experimental melting temperature using a range of state-of-the-art protein force fields. We construct Markov state models in order to examine the thermodynamics, kinetics, mechanism, and rate-determining step of folding. Mechanistically, we find the folding process is rate-limited by the formation of the turn region hydrogen bonds, which occurs following the downhill hydrophobic collapse of the extended denatured protein. These results are presented in the context of established and contradictory theories of the beta-hairpin folding process. Furthermore, our analysis suggests that the AMBER-FB15 force field, at this temperature, best describes the characteristics of the full experimental CLN025 conformational ensemble, while the AMBER ff99SB-ILDN and CHARMM22* force fields display a tendency to overstabilize the native state.

FF12MC: A revised AMBER forcefield and new protein simulation protocol

Specialized to simulate proteins in molecular dynamics (MD) simulations with explicit solvation, FF12MC is a combination of a new protein simulation protocol employing uniformly reduced atomic masses by tenfold and a revised AMBER forcefield FF99 with (i) shortened CH bonds, (ii) removal of torsions involving a nonperipheral sp(3) atom, and (iii) reduced 1-4 interaction scaling factors of torsions ϕ and ψ. This article reports that in multiple, distinct, independent, unrestricted, unbiased, isobaric-isothermal, and classical MD simulations FF12MC can (i) simulate the experimentally observed flipping between left- and right-handed configurations for C14-C38 of BPTI in solution, (ii) autonomously fold chignolin, CLN025, and Trp-cage with folding times that agree with the experimental values, (iii) simulate subsequent unfolding and refolding of these miniproteins, and (iv) achieve a robust Z score of 1.33 for refining protein models TMR01, TMR04, and TMR07. By comparison, the latest general-purpose AMBER forcefield FF14SB locks the C14-C38 bond to the right-handed configuration in solution under the same protein simulation conditions. Statistical survival analysis shows that FF12MC folds chignolin and CLN025 in isobaric-isothermal MD simulations 2-4 times faster than FF14SB under the same protein simulation conditions. These results suggest that FF12MC may be used for protein simulations to study kinetics and thermodynamics of miniprotein folding as well as protein structure and dynamics. Proteins 2016; 84:1490-1516. © 2016 The Authors Proteins: Structure, Function, and Bioinformatics Published by Wiley Periodicals, Inc.

Relaxation mode analysis and Markov state relaxation mode analysis for chignolin in aqueous solution near a transition temperature

It is important to extract reaction coordinates or order parameters from protein simulations in order to investigate the local minimum-energy states and the transitions between them. The most popular method to obtain such data is principal component analysis, which extracts modes of large conformational fluctuations around an average structure. We recently applied relaxation mode analysis for protein systems, which approximately estimates the slow relaxation modes and times from a simulation and enables investigations of the dynamic properties underlying the structural fluctuations of proteins. In this study, we apply this relaxation mode analysis to extract reaction coordinates for a system in which there are large conformational changes such as those commonly observed in protein folding/unfolding. We performed a 750-ns simulation of chignolin protein near its folding transition temperature and observed many transitions between the most stable, misfolded, intermediate, and unfolded states. We then applied principal component analysis and relaxation mode analysis to the system. In the relaxation mode analysis, we could automatically extract good reaction coordinates. The free-energy surfaces provide a clearer understanding of the transitions not only between local minimum-energy states but also between the folded and unfolded states, even though the simulation involved large conformational changes. Moreover, we propose a new analysis method called Markov state relaxation mode analysis. We applied the new method to states with slow relaxation, which are defined by the free-energy surface obtained in the relaxation mode analysis. Finally, the relaxation times of the states obtained with a simple Markov state model and the proposed Markov state relaxation mode analysis are compared and discussed.

Nontargeted Parallel Cascade Selection Molecular Dynamics for Enhancing the Conformational Sampling of Proteins

Nontargeted parallel cascade selection molecular dynamics (nt-PaCS-MD) is proposed as an efficient conformational sampling method to enhance the conformational transitions of proteins, which is an extension of the original targeted PaCS-MD (t-PaCS-MD). The original PaCS-MD comprises cycles of (i) selection of initial structures for multiple independent MD simulations toward a predetermined target and (ii) conformational sampling by the independent MDs. In nt-PaCS-MD, structures that significantly deviate from an average are regarded as candidates that have high potential to address other metastable states and are chosen as the initial structures in the selection. To select significantly deviated structures, we examine the root-mean-square deviation (RMSD) of snapshots generated from the average structure based on Gram-Schmidt orthogonalization. nt-PaCS-MD was applied to the folding of the mini-protein chignolin in implicit solvent and to the open-closed conformational transitions of T4 lysozyme (T4L) and glutamine binding protein (QBP) in explicit solvent. We show that nt-PaCS-MD can reach chignolin's native state and can also cause the open-closed transition of T4L and QBP on a nanosecond time scale, which are very efficient in terms of conformational sampling and comparable to that with t-PaCS-MD.

Enhanced conformational sampling method for proteins based on the TaBoo SeArch algorithm: application to the folding of a mini-protein, chignolin

The conformational samplings are indispensible for obtaining reliable canonical ensembles, which provide statistical averages of physical quantities such as free energies. However, the samplings of vast conformational space of biomacromolecules by conventional molecular dynamics (MD) simulations might be insufficient, due to their inadequate accessible time-scales for investigating biological functions. Therefore, the development of methodologies for enhancing the conformational sampling of biomacromolecules still remains as a challenging issue in computational biology. To tackle this problem, we newly propose an efficient conformational search method, which is referred as TaBoo SeArch (TBSA) algorithm. In TBSA, an inverse energy histogram is used to select seeds for the conformational resampling so that states with high frequencies are inhibited, while states with low frequencies are efficiently sampled to explore the unvisited conformational space. As a demonstration, TBSA was applied to the folding of a mini-protein, chignolin, and automatically sampled the native structure (Cα root mean square deviation < 1.0 Å) with nanosecond order computational costs started from a completely extended structure, although a long-time 1-µs normal MD simulation failed to sample the native structure. Furthermore, a multiscale free energy landscape method based on the conformational sampling of TBSA were quantitatively evaluated through free energy calculations with both implicit and explicit solvent models, which enable us to find several metastable states on the folding landscape.

WExplore: hierarchical exploration of high-dimensional spaces using the weighted ensemble algorithm

As most relevant motions in biomolecular systems are inaccessible to conventional molecular dynamics simulations, algorithms that enhance sampling of rare events are indispensable. Increasing interest in intrinsically disordered systems and the desire to target ensembles of protein conformations (rather than single structures) in drug development motivate the need for enhanced sampling algorithms that are not limited to "two-basin" problems, and can efficiently determine structural ensembles. For systems that are not well-studied, this must often be done with little or no information about the dynamics of interest. Here we present a novel strategy to determine structural ensembles that uses dynamically defined sampling regions that are organized in a hierarchical framework. It is based on the weighted ensemble algorithm, where an ensemble of copies of the system ("replicas") is directed to new regions of configuration space through merging and cloning operations. The sampling hierarchy allows for a large number of regions to be defined, while using only a small number of replicas that can be balanced over multiple length scales. We demonstrate this algorithm on two model systems that are analytically solvable and examine the 10-residue peptide chignolin in explicit solvent. The latter system is analyzed using a configuration space network, and novel hydrogen bonds are found that facilitate folding.

Parallel Cascade Selection Molecular Dynamics (PaCS-MD) to generate conformational transition pathway

Parallel Cascade Selection Molecular Dynamics (PaCS-MD) is proposed as a molecular simulation method to generate conformational transition pathway under the condition that a set of "reactant" and "product" structures is known a priori. In PaCS-MD, the cycle of short multiple independent molecular dynamics simulations and selection of the structures close to the product structure for the next cycle are repeated until the simulated structures move sufficiently close to the product. Folding of 10-residue mini-protein chignolin from the extended to native structures and open-close conformational transition of T4 lysozyme were investigated by PaCS-MD. In both cases, tens of cycles of 100-ps MD were sufficient to reach the product structures, indicating the efficient generation of conformational transition pathway in PaCS-MD with a series of conventional MD without additional external biases. Using the snapshots along the pathway as the initial coordinates, free energy landscapes were calculated by the combination with multiple independent umbrella samplings to statistically elucidate the conformational transition pathways.

QM/MM investigations of organic chemistry oriented questions

About 35 years after its first suggestion, QM/MM became the standard theoretical approach to investigate enzymatic structures and processes. The success is due to the ability of QM/MM to provide an accurate atomistic picture of enzymes and related processes. This picture can even be turned into a movie if nuclei-dynamics is taken into account to describe enzymatic processes. In the field of organic chemistry, QM/MM methods are used to a much lesser extent although almost all relevant processes happen in condensed matter or are influenced by complicated interactions between substrate and catalyst. There is less importance for theoretical organic chemistry since the influence of nonpolar solvents is rather weak and the effect of polar solvents can often be accurately described by continuum approaches. Catalytic processes (homogeneous and heterogeneous) can often be reduced to truncated model systems, which are so small that pure quantum-mechanical approaches can be employed. However, since QM/MM becomes more and more efficient due to the success in software and hardware developments, it is more and more used in theoretical organic chemistry to study effects which result from the molecular nature of the environment. It is shown by many examples discussed in this review that the influence can be tremendous, even for nonpolar reactions. The importance of environmental effects in theoretical spectroscopy was already known. Due to its benefits, QM/MM can be expected to experience ongoing growth for the next decade.In the present chapter we give an overview of QM/MM developments and their importance in theoretical organic chemistry, and review applications which give impressions of the possibilities and the importance of the relevant effects. Since there is already a bunch of excellent reviews dealing with QM/MM, we will discuss fundamental ingredients and developments of QM/MM very briefly with a focus on very recent progress. For the applications we follow a similar strategy.

Population based reweighting of scaled molecular dynamics

Molecular dynamics simulation using enhanced sampling methods is one of the powerful computational tools used to explore protein conformations and free energy landscapes. Enhanced sampling methods often employ either an increase in temperature or a flattening of the potential energy surface to rapidly sample phase space, and a corresponding reweighting algorithm is used to recover the Boltzmann statistics. However, potential energies of complex biomolecules usually involve large fluctuations on a magnitude of hundreds of kcal/mol despite minimal structural changes during simulation. This leads to noisy reweighting statistics and complicates the obtainment of accurate final results. To overcome this common issue in enhanced conformational sampling, we propose a scaled molecular dynamics method, which modifies the biomolecular potential energy surface and employs a reweighting scheme based on configurational populations. Statistical mechanical theory is applied to derive the reweighting formula, and the canonical ensemble of simulated structures is recovered accordingly. Test simulations on alanine dipeptide and the fast folding polypeptide Chignolin exhibit sufficiently enhanced conformational sampling and accurate recovery of free energy surfaces and thermodynamic properties. The results are comparable to long conventional molecular dynamics simulations and exhibit better recovery of canonical statistics over methods which employ a potential energy term in reweighting.

Raising the speed limit for β-hairpin formation

Understanding the folding of the β-hairpin is a crucial step in studying how β-rich proteins fold. We have studied CLN025, an optimized ten residue synthetic peptide, which adopts a compact, well-structured β-hairpin conformation. Formation of the component β-sheet and β-turn structures of CLN025 was probed independently using a combination of equilibrium Fourier transform infrared spectroscopy and laser-induced temperature jump coupled with time-resolved infrared and fluorescence spectroscopies. We find that CLN025 is an ultrafast folder due to its small free energy barrier to folding and that it exceeds the predicted speed limit for β-hairpin formation by an order of magnitude. We also find that the folding mechanism cannot be described by a simple two-state model, but rather is a heterogeneous process involving two independent parallel processes. Formation of stabilizing cross-strand hydrophobic interactions and turn alignment occur competitively, with relaxation lifetimes of 82 ± 10 and 124 ± 10 ns, respectively, at the highest probed temperature. The ultrafast and heterogeneous folding kinetics observed for CLN025 provide evidence for folding on a nearly barrierless free energy landscape, and recalibrate the speed limit for the formation of a β-hairpin.

Turn-directed folding dynamics of β-hairpin-forming de novo decapeptide Chignolin

Realistic mechanistic pictures of β-hairpin formation, offering valuable insights into some of the key early events in protein folding, are accessible through short designed polypeptides as they allow atomic-level scrutiny through simulations. Here, we present a detailed picture of the dynamics and mechanism of β-hairpin formation of Chignolin, a de novo decapeptide, using extensive, unbiased molecular dynamics simulations. The results provide clear evidence for turn-directed broken-zipper folding and reveal details of turn nucleation and cooperative progression of turn growth, hydrogen-bond formations, and eventual packing of the hydrophobic core. Further, we show that, rather than driving folding through hydrophobic collapse, cross-strand side-chain packing could in fact be rate-limiting as packing frustrations can delay formation of the native hydrophobic core prior to or during folding and even cause relatively long-living misfolded or partially folded states that may nucleate aggregative events in more complex situations. The results support the increasing evidence for turn-centric folding mechanisms for β-hairpin formation suggested recently for GB1 and Peptide 1 based on experiments and simulations but also point to the need for similar examinations of polypeptides with larger numbers of cross-strand hydrophobic residues.

Effect of the thermostat in the molecular dynamics simulation on the folding of the model protein chignolin

Molecular dynamics simulations of the model protein chignolin with explicit solvent were carried out, in order to analyze the influence of the Berendsen thermostat on the evolution and folding of the peptide. The dependence of the peptide behavior on temperature was tested with the commonly employed thermostat scheme consisting of one thermostat for the protein and another for the solvent. The thermostat coupling time of the protein was increased to infinity, when the protein is not in direct contact with the thermal bath, a situation known as minimally invasive thermostat. In agreement with other works, it was observed that only in the last situation the instantaneous temperature of the model protein obeys a canonical distribution. As for the folding studies, it was shown that, in the applications of the commonly utilized thermostat schemes, the systems are trapped in local minima regions from which it has difficulty escaping. With the minimally invasive thermostat the time that the protein needs to fold was reduced by two to three times. These results show that the obstacles to the evolution of the extended peptide to the folded structure can be overcome when the temperature of the peptide is not directly controlled.

PROTEIN FOLDING AS A PHYSICAL STOCHASTIC PROCESS

We model protein folding as a physical stochastic process as follows. The unfolded protein chain is treated as a random coil described by SAW (self-avoiding walk). Folding is induced by hydrophobic forces and other interactions, such as hydrogen bonding, which can be taken into account by imposing conditions on SAW. The resulting model is termed CSAW (conditioned self-avoiding walk). Conceptually, the mathematical basis is a generalized Langevin equation. In practice, the model is implemented on a computer by combining SAW and Monte Carlo. To illustrate the flexibility and capabilities of the model, we consider a number of examples, including folding pathways, elastic properties, helix formation, and collective modes.

Perturbation analyses of intermolecular interactions

Conformational fluctuations of a protein molecule are important to its function, and it is known that environmental molecules, such as water molecules, ions, and ligand molecules, significantly affect the function by changing the conformational fluctuations. However, it is difficult to systematically understand the role of environmental molecules because intermolecular interactions related to the conformational fluctuations are complicated. To identify important intermolecular interactions with regard to the conformational fluctuations, we develop herein (i) distance-independent and (ii) distance-dependent perturbation analyses of the intermolecular interactions. We show that these perturbation analyses can be realized by performing (i) a principal component analysis using conditional expectations of truncated and shifted intermolecular potential energy terms and (ii) a functional principal component analysis using products of intermolecular forces and conditional cumulative densities. We refer to these analyses as intermolecular perturbation analysis (IPA) and distance-dependent intermolecular perturbation analysis (DIPA), respectively. For comparison of the IPA and the DIPA, we apply them to the alanine dipeptide isomerization in explicit water. Although the first IPA principal components discriminate two states (the α state and PPII (polyproline II) + β states) for larger cutoff length, the separation between the PPII state and the β state is unclear in the second IPA principal components. On the other hand, in the large cutoff value, DIPA eigenvalues converge faster than that for IPA and the top two DIPA principal components clearly identify the three states. By using the DIPA biplot, the contributions of the dipeptide-water interactions to each state are analyzed systematically. Since the DIPA improves the state identification and the convergence rate with retaining distance information, we conclude that the DIPA is a more practical method compared with the IPA. To test the feasibility of the DIPA for larger molecules, we apply the DIPA to the ten-residue chignolin folding in explicit water. The top three principal components identify the four states (native state, two misfolded states, and unfolded state) and their corresponding eigenfunctions identify important chignolin-water interactions to each state. Thus, the DIPA provides the practical method to identify conformational states and their corresponding important intermolecular interactions with distance information.

Exploring the folding free energy landscape of a β-hairpin miniprotein, chignolin, using multiscale free energy landscape calculation method

The folding process for a β-hairpin miniprotein, chignolin, was investigated by free energy landscape (FEL) calculations using the recently proposed multiscale free energy landscape calculation method (MSFEL). First, coarse-grained molecular dynamics simulations searched a broad conformational space, then multiple independent, all-atom molecular dynamics simulations with explicit solvent determined the detailed local FEL using massively distributed computing. The combination of the two models enabled efficient calculation of the free energy landscapes. The MSFEL analysis showed that chignolin has an intermediate state as well as a misfolded state. The folding process is initiated by the formation of a β-hairpin turn, followed by the formation of contacts in the hydrophobic core between Tyr2 and Trp9. Furthermore, mutation of Tyr2 shifts the population to the misfolded conformation. The results indicate that the hydrophobic core plays an important role in stabilizing the native state of chignolin.

The CLN025 decapeptide retains a β-hairpin conformation in urea and guanidinium chloride

The conformational stability of the β-hairpin miniprotein, CLN025, a variant of chignolin in which the N- and C-terminal glycines are replaced by tyrosines, in various concentrations of guanidinium chloride (GdmCl) and urea was examined by molecular dynamics simulations and electronic circular dichroism (ECD) spectropolarimetry. The peptide maintains its β-hairpin conformation in GdmCl and urea solutions. In GdmCl, Gly7 influences the turn to reduce the number of Asp3-Gly7 H-bonds and the Tyr1-Trp9 H-bond is lost. The structure of the peptide is less stable in 3 M GdmCl than in water or 6 M GdmCl, because the number of Asp3-Thr8 and Tyr1-Tyr10 H-bonds are reduced and the Tyr2 side chain moves away from the Pro4 and Trp9 side chains and toward the Tyr10 side chain. This reduces the number of Tyr2-Pro4 CH-π interactions and Tyr2-Trp9 and Tyr1-Tyr10 aromatic-aromatic (Ar-Ar) interactions and increases the number of Tyr2-Tyr10 Ar-Ar interactions. In 6 M GdmCl at 300 and 333 K, the number of Tyr1-Tyr10 and Asp3-Thr8 H-bonds increases, but fewer structures have Tyr2-Pro4 CH-π and Tyr1-Tyr10 and Tyr2-Trp9 Ar-Ar interactions. In urea, Gly7 is in a mixture of β-turn and random meander structures and the number of Asp3-Thr6 and Tyr1-Tyr10 H-bonds are reduced as are the number of Tyr2-Pro4 CH-π interactions and Tyr1-Tyr10 and Tyr2-Trp9 Ar-Ar interactions. In 4 M urea, a shorter turn places Gly7 into the β-sheet region and Tyr10 is pushed out into the solvent. In 8 M urea, the number of Asp3-Glu5 H-bonds is increased and the β-sheet is lost, but the electrostatic interaction between the charged termini is restored and a cation-π interaction between the indolyl ring of Trp9 and the positively charged N-terminus is formed. In 8 M urea at 333 K, the β-hairpin conformation is almost lost. The structure of CLN025 is stable, because the weakly polar interactions and H-bonds maintain the β-hairpin conformation in the various environments. CLN025 should not be considered a miniprotein, because it lacks a well-defined tertiary structure, it is resistant to denaturation, it does not have an increased heat capacity near its melting temperature, and the structures near and above the melting temperature retain a β-hairpin conformation.

Dual folding pathways of an alpha/beta protein from all-atom ab initio folding simulations

Successful ab initio folding of proteins with both alpha-helix and beta-sheet requires a delicate balance among a variety of forces in the simulation model, which may explain that the successful folding of any alpha/beta proteins to within experimental error has yet to be reported. Here we demonstrate that it is an achievable goal to fold alpha/beta proteins with a force field emphasizing the balance between the two major secondary structures. Using our newly developed force field, we conducted extensive ab initio folding simulations on an alpha/beta protein full sequence design (FSD) employing both conventional molecular dynamics and replica exchange molecular dynamics in combination with a generalized-Born solvation model. In these simulations, the folding of FSD to the native state with high population (>64.2%) and high fidelity (C(alpha)-Root Mean Square Deviation of 1.29 A for the most sampled conformation when compared to the experimental structure) was achieved. The folding of FSD was found to follow two pathways. In the major pathway, the folding started from the formation of the helix. In the minor pathway, however, folding of the beta-hairpin started first. Further examination revealed that the helix initiated from the C-terminus and propagated toward the N-terminus. The formation of the hydrophobic contacts coincided with the global folding. Therefore the hydrophobic force does not appear to be the driving force of the folding of this protein.