Why Is Arginine the Only Amino Acid That Inhibits Polyglutamine Monomers from Taking on Toxic Conformations?

Polyglutamine (polyQ) diseases are devastating neurodegenerative disorders characterized by abnormal expansion of glutamine repeats within specific proteins. The aggregation of polyQ proteins is a critical pathological hallmark of these diseases. Arginine was identified as a promising inhibitory compound because it prevents polyQ-protein monomers from forming intra- and intermolecular β-sheet structures and hinders polyQ proteins from aggregating to form oligomers. Such an aggregation inhibitory effect was not observed in other amino acids. However, the underlying molecular mechanism of the aggregation inhibition and the factors that differentiate arginine from other amino acids, in terms of the inhibition of the polyQ-protein aggregation, remain poorly understood. Here, we performed replica-permutation molecular dynamics simulations to elucidate the molecular mechanism by which arginine inhibits the formation of the intramolecular β-sheet structure of a polyQ monomer. We found that the intramolecular β-sheet structure with more than four β-bridges of the polyQ monomer with arginine is more unstable than without any ligand and with lysine. We also found that arginine has 1.6-2.1 times more contact with polyQ than lysine. In addition, we revealed that arginine forms more hydrogen bonds with the main chain of the polyQ monomer than lysine. More hydrogen bonds formed between arginine and polyQ inhibit polyQ from forming the long intramolecular β-sheet structure. It is known that intramolecular β-sheet structure enhances intermolecular β-sheet structure between proteins. These effects are thought to be the reason for the inhibition of polyQ aggregation. This study provides insights into the molecular events underlying arginine's inhibition of polyQ-protein aggregation.

Vendi sampling for molecular simulations: Diversity as a force for faster convergence and better exploration

Molecular dynamics (MD) is the method of choice for understanding the structure, function, and interactions of molecules. However, MD simulations are limited by the strong metastability of many molecules, which traps them in a single conformation basin for an extended amount of time. Enhanced sampling techniques, such as metadynamics and replica exchange, have been developed to overcome this limitation and accelerate the exploration of complex free energy landscapes. In this paper, we propose Vendi Sampling, a replica-based algorithm for increasing the efficiency and efficacy of the exploration of molecular conformation spaces. In Vendi sampling, replicas are simulated in parallel and coupled via a global statistical measure, the Vendi Score, to enhance diversity. Vendi sampling allows for the recovery of unbiased sampling statistics and dramatically improves sampling efficiency. We demonstrate the effectiveness of Vendi sampling in improving molecular dynamics simulations by showing significant improvements in coverage and mixing between metastable states and convergence of free energy estimates for four common benchmarks, including Alanine Dipeptide and Chignolin.

Dissociation process of polyalanine aggregates by free electron laser irradiation

Polyalanine (polyA) disease-causative proteins with an expansion of alanine repeats can be aggregated. Although curative treatments for polyA diseases have not been explored, the dissociation of polyA aggregates likely reduces the cytotoxicity of polyA. Mid-infrared free electron laser (FEL) successfully dissociated multiple aggregates. However, whether the FEL dissociates polyA aggregates like other aggregates has not been tested. Here, we show that FEL at 6.1 μm experimentally weakened the extent of aggregation of a peptide with 13 alanine repeats (13A), and the irradiated 13A exerted lesser cytotoxicity to neuron-like cells than non-irradiated 13A. Then, we applied molecular dynamics (MD) simulation to follow the dissociation process by FEL. We successfully observed how the intermolecular β-sheet of polyA aggregates was dissociated and separated into monomers with helix structures upon FEL irradiation. After the dissociation by FEL, water molecules inhibited the reformation of polyA aggregates. We recently verified the same dissociation process using FEL-treated amyloid-β aggregates. Thus, a common mechanism underlies the dissociation of different protein aggregates that cause different diseases, polyA disease and Alzheimer's disease. However, MD simulation indicated that polyA aggregates are less easily dissociated than amyloid-β aggregates and require longer laser irradiation due to hydrophobic alanine repeats.

The Double-Layered Structure of Amyloid-β Assemblage on GM1-Containing Membranes Catalytically Promotes Fibrillization

Alzheimer's disease (AD) is associated with progressive accumulation of amyloid-β (Aβ) cross-β fibrils in the brain. Aβ species tightly associated with GM1 ganglioside, a glycosphingolipid abundant in neuronal membranes, promote amyloid fibril formation; therefore, they could be attractive clinical targets. However, the active conformational state of Aβ in GM1-containing lipid membranes is still unknown. The present solid-state nuclear magnetic resonance study revealed a nonfibrillar Aβ assemblage characterized by a double-layered antiparallel β-structure specifically formed on GM1 ganglioside clusters. Our data show that this unique assemblage was not transformed into fibrils on GM1-containing membranes but could promote conversion of monomeric Aβ into fibrils, suggesting that a solvent-exposed hydrophobic layer provides a catalytic surface evoking Aβ fibril formation. Our findings offer structural clues for designing drugs targeting catalytically active Aβ conformational species for the development of anti-AD therapeutics.

Ingenuity in performing replica permutation: How to order the state labels for improving sampling efficiency

In the replica-permutation method, an advanced version of the replica-exchange method, all combinations of replicas and parameters are considered for parameter permutation, and a list of all the combinations is prepared. Here, we report that the temperature transition probability depends on how the list is created, especially in replica permutation with solute tempering (RPST). We found that the transition probabilities decrease at large replica indices when the combinations are sequentially assigned to the state labels as in the originally proposed list. To solve this problem, we propose to modify the list by randomly assigning the combinations to the state labels. We performed molecular dynamics simulations of amyloid-β(16-22) peptides using RPST with the "randomly assigned" list (RPST-RA) and RPST with the "sequentially assigned" list (RPST-SA). The results show the decreases in the transition probabilities in RPST-SA are eliminated, and the sampling efficiency is improved in RPST-RA.

Sigmoid Accelerated Molecular Dynamics: An Efficient Enhanced Sampling Method for Biosystems

Gaussian accelerated molecular dynamics (GaMD) is recognized as a popular enhanced sampling method for tackling long-standing challenges in biomolecular simulations. Inspired by GaMD, Sigmoid accelerated molecular dynamics (SaMD) is proposed in this work by adding a Sigmoid boost potential to improve the balance between the highest acceleration and accurate reweighting. Compared with GaMD, SaMD extends the accessible time scale and improves the computational efficiency as tested in three tasks. In the alanine dipeptide task, SaMD can produce the free energy landscape with better accuracy and efficiency. In the chignolin folding task, the estimated Gibbs free energy difference can converge to the experimental value ∼30% faster. In the protein-ligand binding task, the bound conformations are closer to the crystal structure with a minimal ligand root-mean-square deviation of 1.7 Å. The binding of the ligand XK263 to the HIV protease is reproduced by SaMD in ∼60% less simulation time.

Probing protein misfolding and dissociation with an infrared free-electron laser

Misfolding is observed in the mutant proteins that are causative for neurodegenerative disorders such as polyglutamine diseases. These proteins are prone to aggregate in the cytoplasm and nucleus of cells. To reproduce cells with the aggregated proteins, gene expression system is usually applied, in which the expression construct having the mutated DNA sequence of the interest is transfected into cells. The transfected DNA is finally converted into the mutant protein, which is gradually aggregated in the cells. In addition, a simple method to prepare the cells having aggregates inside has been recently applied. Peptides were first aggregated by incubating them in water. The aggregates are spontaneously taken up by cells because aggregated proteins generally transfer between cells. Peptides with different degrees of aggregation can be made by changing the incubation times and temperatures, which enables to examine contribution of aggregation to the toxicity to the recipient cells. Moreover, such cells can be used for therapeutic researches of diseases in which aggregates are involved. In this chapter, we show methods to induce aggregation of peptides. The functional analyses of the cells with aggregates are also described. Then, experimental dissociation of the aggregates produced using this method by mid infrared free electron laser irradiation and its theoretical support by molecular dynamics simulation are introduced as the therapeutic research for neurodegenerative disorders.

LINES: Log-Probability Estimation via Invertible Neural Networks for Enhanced Sampling

It is very challenging to sample a molecular process with large activation energies using molecular dynamics simulations. Current enhanced sampling methodologies, such as umbrella sampling and metadynamics, rely on the identification of appropriate reaction coordinates for a system. In this paper, we developed a method for log-probability estimation via invertible neural networks for enhanced sampling (LINES). This iterative scheme utilizes a normalizing flow machine learning model to learn the underlying free energy surface (FES) of a system as a function of molecular coordinates and then applies a gradient-based optimization method to the learned normalizing flow to identify reaction coordinates. A biasing potential is then evaluated over a tabulated grid of the reaction coordinate values, which can be applied to the next round of simulations for enhanced sampling, resulting in more efficient sampling. We tested the accuracy and efficiency of the LINES method in sampling the FES using the alanine dipeptide system. We also demonstrated the effectiveness of identification of reaction coordinates through simulation of cyclobutanol unbinding from β-cyclodextrin and the folding/unfolding of CLN025─a variant of the peptide Chignolin. The LINES method can be extended to the study of large-scale protein systems with complex nonlinear reaction pathways.

State-of-the-Art Molecular Dynamics Simulation Studies of RNA-Dependent RNA Polymerase of SARS-CoV-2

Molecular dynamics (MD) simulations are powerful theoretical methods that can reveal biomolecular properties, such as structure, fluctuations, and ligand binding, at the level of atomic detail. In this review article, recent MD simulation studies on these biomolecular properties of the RNA-dependent RNA polymerase (RdRp), which is a multidomain protein, of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) are presented. Although the tertiary structures of RdRps in SARS-CoV-2 and SARS-CoV are almost identical, the RNA synthesis activity of RdRp of SARS-CoV is higher than SARS-CoV-2. Recent MD simulations observed a difference in the dynamic properties of the two RdRps, which may cause activity differences. RdRp is also a drug target for Coronavirus disease 2019 (COVID-19). Nucleotide analogs, such as remdesivir and favipiravir, are considered to be taken up by RdRp and inhibit RNA replication. Recent MD simulations revealed the recognition mechanism of RdRp for these drug molecules and adenosine triphosphate (ATP). The ligand-recognition ability of RdRp decreases in the order of remdesivir, favipiravir, and ATP. As a typical recognition process, it was found that several lysine residues of RdRp transfer these ligand molecules to the binding site such as a “bucket brigade.” This finding will contribute to understanding the mechanism of the efficient ligand recognition by RdRp. In addition, various simulation studies on the complexes of SARS-CoV-2 RdRp with several nucleotide analogs are reviewed, and the molecular mechanisms by which these compounds inhibit the function of RdRp are discussed. The simulation studies presented in this review will provide useful insights into how nucleotide analogs are recognized by RdRp and inhibit the RNA replication.

Investigation of the conformation of human prion protein in ethanol solution using molecular dynamics simulations

When the conformation of protein is changed from its natural state to a misfolded state, some diseases will happen like prion disease. Prion diseases are a set of deadly neurodegenerative diseases caused by prion protein misfolding and aggregation. Monohydric alcohols have a strong influence on the structure of protein. However, whether monohydric alcohols inhibit amyloid fibrosis remains uncertain. Here, to elucidate the effect of ethanol on the structural stability of human prion protein, molecular dynamics simulations were employed to analyze the conformational changes and dynamics characteristics of human prion proteins at different temperatures. The results show that the extension of β-sheet occurs more easily and the α-helix is more easily disrupted at high temperatures. We found that ethanol can destroy the hydrophobic interactions and make the hydrogen bonds stable, which protects the secondary structure of the protein, especially at 500 K.Communicated by Ramaswamy H. Sarma.

Molecular Dynamics Simulation Studies on the Aggregation of Amyloid-β Peptides and Their Disaggregation by Ultrasonic Wave and Infrared Laser Irradiation

Alzheimer’s disease is understood to be caused by amyloid fibrils and oligomers formed by aggregated amyloid-β (Aβ) peptides. This review article presents molecular dynamics (MD) simulation studies of Aβ peptides and Aβ fragments on their aggregation, aggregation inhibition, amyloid fibril conformations in equilibrium, and disruption of the amyloid fibril by ultrasonic wave and infrared laser irradiation. In the aggregation of Aβ, a β-hairpin structure promotes the formation of intermolecular β-sheet structures. Aβ peptides tend to exist at hydrophilic/hydrophobic interfaces and form more β-hairpin structures than in bulk water. These facts are the reasons why the aggregation is accelerated at the interface. We also explain how polyphenols, which are attracting attention as aggregation inhibitors of Aβ peptides, interact with Aβ. An MD simulation study of the Aβ amyloid fibrils in equilibrium is also presented: the Aβ amyloid fibril has a different structure at one end from that at the other end. The amyloid fibrils can be destroyed by ultrasonic wave and infrared laser irradiation. The molecular mechanisms of these amyloid fibril disruptions are also explained, particularly focusing on the function of water molecules. Finally, we discuss the prospects for developing treatments for Alzheimer’s disease using MD simulations.

Tardigrade Secretory-Abundant Heat-Soluble Protein Varies Entrance Propensity Depending on the Amino-Acid Sequence

Secretory-abundant heat-soluble (SAHS) proteins, which constitute a protein family unique to tardigrades, are thought to be essential for anhydrobiosis. Our previous study has revealed that one of the SAHS proteins of Ramazzottius varieornatus (RvSAHS1) has a more flexible entrance than a mammalian fatty-acid-binding protein, which has a crystal structure similar to that of RvSAHS1. Recently, SAHS paralogs that are expressed abundantly and specifically in the early embryos of this tardigrade and Hypsibius exemplaris have been identified. Comparing these amino-acid sequences with that of RvSAHS1, we have found characteristic differences as I113F and D146T. In this study, we investigate I113F and D146T mutants' properties of RvSAHS1 using molecular dynamics simulations and compare the structures and fluctuations of their entrances with those of the wild type. The two mutants exhibit different properties at the entrance of the β-barrel structure. The I113F mutant tends to close the entrance more than the wild type due to the enhanced hydrophobic network inside the cavity. The D146T mutant, in contrast to the I113F mutant, tends to open the entrance. The mechanism by which this mutation opens the entrance is also discussed. Even though only a single mutation located far from the entrance is added to the wild type, there is a clear difference in the tendency to open and close the β-barrel entrance. It indicates that the entrance properties of the SAHS protein are sensitive to the amino-acid sequence.

Replica permutation with solute tempering for molecular dynamics simulation and its application to the dimerization of amyloid-β fragments

We propose the replica permutation with solute tempering (RPST) by combining the replica-permutation method (RPM) and the replica exchange with solute tempering (REST). Temperature permutations are performed among more than two replicas in RPM, whereas temperature exchanges are performed between two replicas in the replica-exchange method (REM). The temperature transition in RPM occurs more efficiently than in REM. In REST, only the temperatures of the solute region, the solute temperatures, are exchanged to reduce the number of replicas compared to REM. Therefore, RPST is expected to be an improved method taking advantage of these methods. For comparison, we applied RPST, REST, RPM, and REM to two amyloid-β(16–22) peptides in explicit water. We calculated the transition ratio and the number of tunneling events in the temperature space and the number of dimerization events of amyloid-β(16–22) peptides. The results indicate that, in RPST, the number of replicas necessary for frequent random walks in the temperature and conformational spaces is reduced compared to the other three methods. In addition, we focused on the dimerization process of amyloid-β(16–22) peptides. The RPST simulation with a relatively small number of replicas shows that the two amyloid-β(16–22) peptides form the intermolecular antiparallel β-bridges due to the hydrophilic side-chain contact between Lys and Glu and hydrophobic side-chain contact between Leu, Val, and Phe, which stabilizes the dimer of the peptides.

All-Atom Molecular Dynamics Simulation Methods for the Aggregation of Protein and Peptides: Replica Exchange/Permutation and Nonequilibrium Simulations

Protein aggregates are associated with more than 40 serious human diseases. To understand the formation mechanism of protein aggregates at atomic level, all-atom molecular dynamics (MD) simulation is a powerful computational tool. In this chapter, we review the all-atom MD simulation methods that are useful for study on the protein aggregation. We first explain conventional MD simulation methods in physical statistical ensembles, such as the canonical and isothermal-isobaric ensembles. We then describe the generalized-ensemble algorithms such as replica-exchange and replica-permutation MD methods. These methods can overcome a difficulty, in which simulations tend to get trapped in local-minimum free-energy states. Finally we explain the nonequilibrium MD method. Some simulation results based on these methods are also presented.

Stress-dependent flexibility of a full-length human monoclonal antibody: Insights from molecular dynamics to support biopharmaceutical development

After several decades of advancements in drug discovery, product development of biopharmaceuticals remains a time- and resource-consuming endeavor. One of the main reasons is associated to the lack of fundamental understanding of conformational dynamics of such biologic entities, and how they respond to various stresses encountered during manufacturing. In this work, we have studied the conformational dynamics of human IgG₁κ b12 monoclonal antibody (mAb) using molecular dynamics simulations. The hundreds of nanoseconds long trajectories reveal that b12 mAb is highly flexible. Its variable domains show greater conformational fluctuations than the constant domains. Additionally, it collapses towards a more globular shape in response to thermal stress, leading to decrease in the total solvent exposed surface area and radius of gyration. This behavior is more pronounced for the deglycosylated b12 mAb, and it appears to correlate with increase in inter-domain contacts between specific regions of the antibody. Conformational fluctuations also cause temporary formation and disruption of hydrophobic and charged patches on the antibody surface, which is particularly important for the prediction of CMC properties during development phases of antibody-based biotherapeutics. The insights gained through these simulations may help the development of biologic drugs, especially with regards to manufacturing processes where antibodies may undergo significant thermal stress.

Deep learning the slow modes for rare events sampling

The development of enhanced sampling methods has greatly extended the scope of atomistic simulations, allowing long-time phenomena to be studied with accessible computational resources. Many such methods rely on the identification of an appropriate set of collective variables. These are meant to describe the system's modes that most slowly approach equilibrium under the action of the sampling algorithm. Once identified, the equilibration of these modes is accelerated by the enhanced sampling method of choice. An attractive way of determining the collective variables is to relate them to the eigenfunctions and eigenvalues of the transfer operator. Unfortunately, this requires knowing the long-term dynamics of the system beforehand, which is generally not available. However, we have recently shown that it is indeed possible to determine efficient collective variables starting from biased simulations. In this paper, we bring the power of machine learning and the efficiency of the recently developed on the fly probability-enhanced sampling method to bear on this approach. The result is a powerful and robust algorithm that, given an initial enhanced sampling simulation performed with trial collective variables or generalized ensembles, extracts transfer operator eigenfunctions using a neural network ansatz and then accelerates them to promote sampling of rare events. To illustrate the generality of this approach, we apply it to several systems, ranging from the conformational transition of a small molecule to the folding of a miniprotein and the study of materials crystallization.

"Bucket brigade" using lysine residues in RNA-dependent RNA polymerase of SARS-CoV-2

The RNA-dependent RNA polymerase (RdRp) of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is a promising drug target for coronavirus disease 2019 (COVID-19) because it plays the most important role in the replication of the RNA genome. Nucleotide analogs such as remdesivir and favipiravir are thought to interfere with the RNA replication by RdRp. More specifically, they are expected to compete with nucleoside triphosphates, such as ATP. However, the process in which these drug molecules and nucleoside triphosphates are taken up by RdRp remains unknown. In this study, we performed all-atom molecular dynamics simulations to clarify the recognition mechanism of RdRp for these drug molecules and ATP that were at a distance. The ligand recognition ability of RdRp decreased in the order of remdesivir, favipiravir, and ATP. We also identified six recognition paths. Three of them were commonly found in all ligands, and the remaining three paths were ligand-dependent ones. In the common two paths, it was observed that the multiple lysine residues of RdRp carried the ligands to the binding site like a "bucket brigade." In the remaining common path, the ligands directly reached the binding site. Our findings contribute to the understanding of the efficient ligand recognition by RdRp at the atomic level.

Dynamic properties of SARS-CoV and SARS-CoV-2 RNA-dependent RNA polymerases studied by molecular dynamics simulations

One of the promising drug targets against COVID-19 is an RNA-dependent RNA polymerase (RdRp) of SARS-CoV-2. The tertiary structures of the SARS-CoV-2 and SARS-CoV RdRps are almost the same. However, the RNA-synthesizing activity of the SARS-CoV RdRp is higher than that of the SARS-CoV-2 RdRp. We performed molecular dynamics simulations and found differences in their dynamic properties. In the SARS-CoV RdRp, motifs A-G, which form the active site, are up to 63% closer to each other. We also observed cooperative domain motion in the SARS-CoV RdRp. Such dynamic differences may cause the activity differences between the two RdRps.

Tardigrade Secretory-Abundant Heat-Soluble Protein Has a Flexible β-Barrel Structure in Solution and Keeps This Structure in Dehydration

Secretory-abundant heat-soluble (SAHS) proteins are unique heat-soluble proteins of Tardigrada and are believed to play an essential role in anhydrobiosis, a latent state of life induced by desiccation. To investigate the dynamic properties, molecular dynamics (MD) simulations of a SAHS protein, RvSAHS1, were performed in solution and under dehydrating conditions. For comparison purposes, MD simulations of a human liver-type fatty-acid binding protein (LFABP) were performed in solution. Furthermore, high-speed atomic force microscopy observations were conducted to ascertain the results of the MD simulations. Three properties of RvSAHS1 were found as follows. (1) The entrance region of RvSAHS1 is more flexible and can be more extensive in solutions compared with that of a human LFABP because there is no salt bridge between the βD and βE strands. (2) The intrinsically disordered domain in the N-terminal region significantly fluctuates and can form an amphiphilic α-helix. (3) The size of the entrance region gets smaller along with dehydration, keeping the β-barrel structure. Overall, the obtained results provide atomic-level dynamics of SAHS proteins.

Molecular dynamics simulations of amyloid-β(16-22) peptide aggregation at air-water interfaces

Oligomers of amyloid-β (Aβ) peptides are known to be related to Alzheimer's disease, and their formation is accelerated at hydrophilic-hydrophobic interfaces, such as the cell membrane surface and air-water interface. Here, we report molecular dynamics simulations of aggregation of Aβ(16-22) peptides at air-water interfaces. First, 100 randomly distributed Aβ(16-22) peptides moved to the interface. The high concentration of peptides then accelerated their aggregation and formation of antiparallel β-sheets. Two layers of oligomers were observed near the interface. In the first layer from the interface, the oligomer with less β-bridges exposed the hydrophobic residues to the air. The second layer consisted of oligomers with more β-bridges that protruded into water. They are more soluble in water because the hydrophobic residues are covered by N- and C-terminal hydrophilic residues that are aligned well along the oligomer edge. These results indicate that amyloid protofibril formation mainly occurs in the second layer.

Role of Water Molecules and Helix Structure Stabilization in the Laser-Induced Disruption of Amyloid Fibrils Observed by Nonequilibrium Molecular Dynamics Simulations

Water plays a crucial role in the formation and destruction of biomolecular structures. The mechanism for destroying biomolecular structures was thought to be an active breaking of hydrogen bonds by water molecules. However, using nonequilibrium molecular dynamics simulations, in which an amyloid-β amyloid fibril was destroyed via infrared free-electron laser (IR-FEL) irradiation, we discovered a new mechanism, in which water molecules disrupt protein aggregates. The intermolecular hydrogen bonds formed by C═O and N-H in the fibril are broken at each pulse of laser irradiation. These bonds spontaneously re-form after the irradiation in many cases. However, when a water molecule happens to enter the gap between C═O and N-H, it inhibits the re-formation of the hydrogen bonds. Such sites become defects in the regularly aligned hydrogen bonds, from which all hydrogen bonds in the intermolecular β-sheet are broken as the fraying spreads. This role of water molecules is entirely different from other known mechanisms. This new mechanism can explain the recent experiments showing that the amyloid fibrils are not destroyed by laser irradiation under dry conditions. Additionally, we found that helix structures form more after the amyloid disruption; this is because the resonance frequency is different in a helix structure. Our findings provide a theoretical basis for the application of IR-FEL to the future treatment of amyloidosis.

Universal Relation between Instantaneous Diffusivity and Radius of Gyration of Proteins in Aqueous Solution

Protein conformational fluctuations are highly complex and exhibit long-term correlations. Here, molecular dynamics simulations of small proteins demonstrate that these conformational fluctuations directly affect the protein's instantaneous diffusivity D_{I}. We find that the radius of gyration R_{g} of the proteins exhibits 1/f fluctuations that are synchronous with the fluctuations of D_{I}. Our analysis demonstrates the validity of the local Stokes-Einstein-type relation D_{I}∝1/(R_{g}+R_{0}), where R_{0}∼0.3  nm is assumed to be a hydration layer around the protein. From the analysis of different protein types with both strong and weak conformational fluctuations, the validity of the Stokes-Einstein-type relation appears to be a general property.

Structural basis for promiscuous action of monoterpenes on TRP channels

Monoterpenes are major constituents of plant-derived essential oils and have long been widely used for therapeutic and cosmetic applications. The monoterpenes menthol and camphor are agonists or antagonists for several TRP channels such as TRPM8, TRPV1, TRPV3 and TRPA1. However, which regions within TRPV1 and TRPV3 confer sensitivity to monoterpenes or other synthesized chemicals such as 2-APB are unclear. In this study we identified conserved arginine and glycine residues in the linker between S4 and S5 that are related to the action of these chemicals and validated these findings in molecular dynamics simulations. The involvement of these amino acids differed between TRPV3 and TRPV1 for chemical-induced and heat-evoked activation. These findings provide the basis for characterization of physiological function and biophysical properties of ion channels.

Dimerization of α-Synuclein Fragments Studied by Isothermal-Isobaric Replica-Permutation Molecular Dynamics Simulation

Aggregates and fibrils of intrinsically disordered α-synuclein are associated with Parkinson's disease. Within a non-amyloid β component (NAC) spanning from the 61st to the 95th residue of α-synuclein, an 11-residue segment called NACore (₆₈GAVVTGVTAVA₇₈) is an essential region for both fibril formation and cytotoxicity. Although NACore peptides alone are known to form aggregates and amyloid fibrils, the mechanisms of aggregation and fibrillation remain unknown. This study investigated the dimerization process of NACore peptides as the initial stage of the aggregation and fibrillation processes. We performed an isothermal-isobaric replica-permutation molecular dynamics simulation, which is one of the efficient sampling methods, for the two NACore peptides in explicit water over 96 μs. The simulation succeeded in sampling a variety of dimer structures. An analysis of secondary structure revealed that most of the NACore dimers form intermolecular β-bridges. In particular, more antiparallel β-bridges were observed than parallel β-bridges. We also found that intramolecular secondary structures such as α-helix and antiparallel β-bridge are stabilized in the pre-dimer state. However, we identified that the intermolecular β-bridges tend to form directly between residues with no specific structure rather than via the intramolecular β-bridges. This is because the NACore peptides still have a low propensity to form the intramolecular secondary structures even though they are stabilized in the pre-dimer state.

Promotion and Inhibition of Amyloid-β Peptide Aggregation: Molecular Dynamics Studies

Aggregates of amyloid-β (Aβ) peptides are known to be related to Alzheimer’s disease. Their aggregation is enhanced at hydrophilic–hydrophobic interfaces, such as a cell membrane surface and air-water interface, and is inhibited by polyphenols, such as myricetin and rosmarinic acid. We review molecular dynamics (MD) simulation approaches of a full-length Aβ peptide, Aβ40, and Aβ(16–22) fragments in these environments. Since these peptides have both hydrophilic and hydrophobic amino acid residues, they tend to exist at the interfaces. The high concentration of the peptides accelerates the aggregation there. In addition, Aβ40 forms a β-hairpin structure, and this structure accelerates the aggregation. We also describe the inhibition mechanism of the Aβ(16–22) aggregation by polyphenols. The aggregation of Aβ(16–22) fragments is caused mainly by the electrostatic attraction between charged amino acid residues known as Lys16 and Glu22. Since polyphenols form hydrogen bonds between their hydroxy and carboxyl groups and these charged amino acid residues, they inhibit the aggregation.

Involvement of pore helix in voltage-dependent inactivation of TRPM5 channel

The transient receptor potential melastatin 5 (TRPM5) channel is a monovalent-permeable cation channel that is activated by intracellular Ca²⁺. Expression of TRPM5 has been shown in taste cells, pancreas, brainstem and olfactory epithelium, and this channel is thought to be involved in controlling membrane potentials. In whole-cell patch-clamp recordings, TRPM5 exhibited voltage-dependent inactivation at negative membrane potentials and time constant of voltage-dependent inactivation of TRPM5 did not depend on the intracellular Ca²⁺ concentrations between 100 and 500 nM. Alanine substitution at Y913 and I916 in the pore helix of TRPM5 increased time constant of voltage-dependent inactivation. Meanwhile, voltage-dependent inactivation was reduced in TRPM5 mutants having glycine substitution at L901, Y913, Q915 and I916 in the pore helix. From these results, we conclude that the pore helix in the outer pore loop might play a role in voltage-dependent inactivation of TRPM5.

Energetics and kinetics of substrate analog-coupled staphylococcal nuclease folding revealed by a statistical mechanical approach

Protein conformational changes associated with ligand binding, especially those involving intrinsically disordered proteins, are mediated by tightly coupled intra- and intermolecular events. Such reactions are often discussed in terms of two limiting kinetic mechanisms, conformational selection (CS), where folding precedes binding, and induced fit (IF), where binding precedes folding. It has been shown that coupled folding/binding reactions can proceed along both CS and IF pathways with the flux ratio depending on conditions such as ligand concentration. However, the structural and energetic basis of such complex reactions remains poorly understood. Therefore, we used experimental, theoretical, and computational approaches to explore structural and energetic aspects of the coupled-folding/binding reaction of staphylococcal nuclease in the presence of the substrate analog adenosine-3',5'-diphosphate. Optically monitored equilibrium and kinetic data, combined with a statistical mechanical model, gave deeper insight into the relative importance of specific and Coulombic protein-ligand interactions in governing the reaction mechanism. We also investigated structural aspects of the reaction at the residue level using NMR and all-atom replica-permutation molecular dynamics simulations. Both approaches yielded clear evidence for accumulation of a transient protein-ligand encounter complex early in the reaction under IF-dominant conditions. Quantitative analysis of the equilibrium/kinetic folding revealed that the ligand-dependent CS-to-IF shift resulted from stabilization of the compact transition state primarily by weakly ligand-dependent Coulombic interactions with smaller contributions from specific binding energies. At a more macroscopic level, the CS-to-IF shift was represented as a displacement of the reaction "route" on the free energy surface, which was consistent with a flux analysis.

In situ and real-time insight into Rhizopus chinensis lipase under high pressure and temperature: Conformational traits and biobehavioural analysis

An in situ and real-time investigation was performed using an optical cell system and in-silico analysis to reveal the impacts of pressure and temperature on the conformational state and behaviours of Rhizopus chinensis lipase (RCL). The fluorescence intensity (FI) of RCL increased remarkably under high pressure, and part of this increase was recovered after depressurization. This result displayed the partially reversible conformational change of RCL, which may be associated with the local change of Trp224 near the catalytic centre. High temperature caused a significant loss of secondary structure, whereas the α-helical segments including the lid were preserved by high pressure even at temperatures over 60 °C. The parameters of enzymatic reaction monitored by UV showed that the hydrolysis rate was remarkably enhanced by the pressure of 200 MPa. In the pressure range of 0.1-200 MPa, the active volume measured by the in situ system decreased from -2.85 to -6.73 mL/mol with the temperature increasing from 20 °C to 40 °C. The high catalytic capacity of the lipase under high pressure and high temperature was primarily attributed to pressure protection on RCL.

Mesophilic Pyrophosphatase Function at High Temperature: A Molecular Dynamics Simulation Study

The mesophilic inorganic pyrophosphatase from Escherichia coli (EcPPase) retains function at 353 K, the physiological temperature of hyperthermophilic Thermococcus thioreducens, whereas the homolog protein (TtPPase) from this hyperthermophilic organism cannot function at room temperature. To explain this asymmetric behavior, we examined structural and dynamical properties of the two proteins using molecular dynamics simulations. The global flexibility of TtPPase is significantly higher than its mesophilic homolog at all tested temperature/pressure conditions. However, at 353 K, EcPPase reduces its solvent-exposed surface area and increases subunit compaction while maintaining flexibility in its catalytic pocket. In contrast, TtPPase lacks this adaptability and has increased rigidity and reduced protein/water interactions in its catalytic pocket at room temperature, providing a plausible explanation for its inactivity near room temperature.

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.

Molecular dynamics simulations for genetic interpretation in protein coding regions: where we are, where to go and when

The increasing ease with which massive genetic information can be obtained from patients or healthy individuals has stimulated the development of interpretive bioinformatics tools as aids in clinical practice. Most such tools analyze evolutionary information and simple physical-chemical properties to predict whether replacement of one amino acid residue with another will be tolerated or cause disease. Those approaches achieve up to 80-85% accuracy as binary classifiers (neutral/pathogenic). As such accuracy is insufficient for medical decision to be based on, and it does not appear to be increasing, more precise methods, such as full-atom molecular dynamics (MD) simulations in explicit solvent, are also discussed. Then, to describe the goal of interpreting human genetic variations at large scale through MD simulations, we restrictively refer to all possible protein variants carrying single-amino-acid substitutions arising from single-nucleotide variations as the human variome. We calculate its size and develop a simple model that allows calculating the simulation time needed to have a 0.99 probability of observing unfolding events of any unstable variant. The knowledge of that time enables performing a binary classification of the variants (stable-potentially neutral/unstable-pathogenic). Our model indicates that the human variome cannot be simulated with present computing capabilities. However, if they continue to increase as per Moore's law, it could be simulated (at 65°C) spending only 3 years in the task if we started in 2031. The simulation of individual protein variomes is achievable in short times starting at present. International coordination seems appropriate to embark upon massive MD simulations of protein variants.

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.

Insight into Improved Thermostability of Cold-Adapted Staphylococcal Lipase by Glycine to Cysteine Mutation

Thermostability remains one of the most desirable traits in many lipases. Numerous studies have revealed promising strategies to improve thermostability and random mutagenesis often leads to unexpected yet interesting findings in engineering stability. Previously, the thermostability of C-terminal truncated cold-adapted lipase from Staphylococcus epidermidis AT2 (rT-M386) was markedly enhanced by directed evolution. The newly evolved mutant, G210C, demonstrated an optimal temperature shift from 25 to 45 °C and stability up to 50 °C. Interestingly, a cysteine residue was randomly introduced on the loop connecting the two lids and accounted for the only cysteine found in the lipase. We further investigated the structural and mechanistic insights that could possibly cause the significant temperature shift. Both rT-M386 and G210C were modeled and simulated at 25 °C and 50 °C. The results clearly portrayed the effect of cysteine substitution primarily on the lid stability. Comparative molecular dynamics simulation analysis revealed that G210C exhibited greater stability than the wild-type at high temperature simulation. The compactness of the G210C lipase structure increased at 50 °C and resulted in enhanced rigidity hence stability. This observation is supported by the improved and stronger non-covalent interactions formed in the protein structure. Our findings suggest that the introduction of a single cysteine residue at the lid region of cold-adapted lipase may result in unexpected increased in thermostability, thus this approach could serve as one of the thermostabilization strategies in engineering lipase stability.

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 dynamics simulation of proteins under high pressure: Structure, function and thermodynamics

BACKGROUND

Molecular dynamics (MD) simulation is well-recognized as a powerful tool to investigate protein structure, function, and thermodynamics. MD simulation is also used to investigate high pressure effects on proteins. For conducting better MD simulation under high pressure, the main issues to be addressed are: (i) protein force fields and water models were originally developed to reproduce experimental properties obtained at ambient pressure; and (ii) the timescale to observe the pressure effect is often much longer than that of conventional MD simulations.

SCOPE OF REVIEW

First, we describe recent developments in MD simulation methodologies for studying the high-pressure structure and dynamics of protein molecules. These developments include force fields for proteins and water molecules, and enhanced simulation techniques. Then, we summarize recent studies of MD simulations of proteins in water under high pressure.

MAJOR CONCLUSIONS

Recent MD simulations of proteins in solution under pressure have reproduced various phenomena identified by experiments using high pressure, such as hydration, water penetration, conformational change, helix stabilization, and molecular stiffening.

GENERAL SIGNIFICANCE

MD simulations demonstrate differences in the properties of proteins and water molecules between ambient and high-pressure conditions. Comparing the results obtained by MD calculations with those obtained experimentally could reveal the mechanism by which biological molecular machines work well in collaboration with water molecules.

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.

Theoretical analysis on thermodynamic stability of chignolin

Understanding the dominant factor in thermodynamic stability of proteins remains an open challenge. Kauzmann's hydrophobic interaction hypothesis, which considers hydrophobic interactions between nonpolar groups as the dominant factor, has been widely accepted for about sixty years and attracted many scientists. The hypothesis, however, has not been verified or disproved because it is difficult, both theoretically and experimentally, to quantify the solvent effects on the free energy change in protein folding. Here, we developed a computational method for extracting the dominant factor behind thermodynamic stability of proteins and applied it to a small, designed protein, chignolin. The resulting free energy profile quantitatively agreed with the molecular dynamics simulations. Decomposition of the free energy profile indicated that intramolecular interactions predominantly stabilized collapsed conformations, whereas solvent-induced interactions, including hydrophobic ones, destabilized them. These results obtained for chignolin were consistent with the site-directed mutagenesis and calorimetry experiments for globular proteins with hydrophobic interior cores.

Development of isothermal-isobaric replica-permutation method for molecular dynamics and Monte Carlo simulations and its application to reveal temperature and pressure dependence of folded, misfolded, and unfolded states of chignolin

We developed a two-dimensional replica-permutation molecular dynamics method in the isothermal-isobaric ensemble. The replica-permutation method is a better alternative to the replica-exchange method. It was originally developed in the canonical ensemble. This method employs the Suwa-Todo algorithm, instead of the Metropolis algorithm, to perform permutations of temperatures and pressures among more than two replicas so that the rejection ratio can be minimized. We showed that the isothermal-isobaric replica-permutation method performs better sampling efficiency than the isothermal-isobaric replica-exchange method and infinite swapping method. We applied this method to a β-hairpin mini protein, chignolin. In this simulation, we observed not only the folded state but also the misfolded state. We calculated the temperature and pressure dependence of the fractions on the folded, misfolded, and unfolded states. Differences in partial molar enthalpy, internal energy, entropy, partial molar volume, and heat capacity were also determined and agreed well with experimental data. We observed a new phenomenon that misfolded chignolin becomes more stable under high-pressure conditions. We also revealed this mechanism of the stability as follows: TYR2 and TRP9 side chains cover the hydrogen bonds that form a β-hairpin structure. The hydrogen bonds are protected from the water molecules that approach the protein as the pressure increases.

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.

Theoretical approaches for dynamical ordering of biomolecular systems

BACKGROUND

Living systems are characterized by the dynamic assembly and disassembly of biomolecules. The dynamical ordering mechanism of these biomolecules has been investigated both experimentally and theoretically. The main theoretical approaches include quantum mechanical (QM) calculation, all-atom (AA) modeling, and coarse-grained (CG) modeling. The selected approach depends on the size of the target system (which differs among electrons, atoms, molecules, and molecular assemblies). These hierarchal approaches can be combined with molecular dynamics (MD) simulation and/or integral equation theories for liquids, which cover all size hierarchies.

SCOPE OF REVIEW

We review the framework of quantum mechanical/molecular mechanical (QM/MM) calculations, AA MD simulations, CG modeling, and integral equation theories. Applications of these methods to the dynamical ordering of biomolecular systems are also exemplified.

MAJOR CONCLUSIONS

The QM/MM calculation enables the study of chemical reactions. The AA MD simulation, which omits the QM calculation, can follow longer time-scale phenomena. By reducing the number of degrees of freedom and the computational cost, CG modeling can follow much longer time-scale phenomena than AA modeling. Integral equation theories for liquids elucidate the liquid structure, for example, whether the liquid follows a radial distribution function.

GENERAL SIGNIFICANCE

These theoretical approaches can analyze the dynamic behaviors of biomolecular systems. They also provide useful tools for exploring the dynamic ordering systems of biomolecules, such as self-assembly. This article is part of a Special Issue entitled "Biophysical Exploration of Dynamical Ordering of Biomolecular Systems" edited by Dr. Koichi Kato.

Structural and fluctuational difference between two ends of Aβ amyloid fibril: MD simulations predict only one end has open conformations

Aβ amyloid fibrils, which are related to Alzheimer’s disease, have a cross-β structure consisting of two β-sheets: β1 and β2. The Aβ peptides are thought to be serially arranged in the same molecular conformation along the fibril axis. However, to understand the amyloid extension mechanism, we must understand the amyloid fibril structure and fluctuation at the fibril end, which has not been revealed to date. Here, we reveal these features by all-atom molecular dynamics (MD) simulations of Aβ42 and Aβ40 fibrils in explicit water. The structure and fluctuation were observed to differ between the two ends. At the even end, the Aβ peptide always took a closed form wherein β1 and β2 were closely spaced. The Aβ peptide fluctuated more at the odd end and took an open form wherein the two β-sheets were well separated. The differences are attributed to the stronger β-sheet formation by the β1 exposed at the even end than the β2 exposed at the odd end. Along with the small fluctuations at the even end, these results explain why the fibril extends from one end only, as observed in experiments. Our MD results agree well with recent observations by high-speed atomic force microscopy.

Structural basis of TRPA1 inhibition by HC-030031 utilizing species-specific differences

Pain is a harmful sensation that arises from noxious stimuli. Transient receptor potential ankyrin 1 (TRPA1) is one target for studying pain mechanisms. TRPA1 is activated by various stimuli such as noxious cold, pungent natural products and environmental irritants. Since TRPA1 is an attractive target for pain therapy, a few TRPA1 antagonists have been developed and some function as analgesic agents. The responses of TRPA1 to agonists and antagonists vary among species and these species differences have been utilized to identify the structural basis of activation and inhibition mechanisms. The TRPA1 antagonist HC-030031 (HC) failed to inhibit frog TRPA1 (fTRPA1) and zebrafish TRPA1 activity induced by cinnamaldehyde (CA), but did inhibit human TRPA1 (hTRPA1) in a heterologous expression system. Chimeric studies between fTRPA1 and hTRPA1, as well as analyses using point mutants, revealed that a single amino acid residue (N855 in hTRPA1) significantly contributes to the inhibitory action of HC. Moreover, the N855 residue and the C-terminus region exhibited synergistic effects on the inhibition by HC. Molecular dynamics simulation suggested that HC stably binds to hTRPA1-N855. These findings provide novel insights into the structure-function relationship of TRPA1 and could lead to the development of more effective analgesics targeted to TRPA1.

Effect of high pressure processing on dispersive and aggregative properties of almond milk

BACKGROUND

A study was conducted to investigate the impact of high pressure (450 and 600 MPa at 30 °C) and thermal (72, 85 and 99 °C at 0.1 MPa) treatments on dispersive and aggregative characteristics of almond milk. Experiments were conducted using a kinetic pressure testing unit and water bath. Particle size distribution, microstructure, UV absorption spectra, pH and color changes of processed and unprocessed samples were analyzed.

RESULTS

Raw almond milk represented the mono model particle size distribution with average particle diameters of 2 to 3 µm. Thermal or pressure treatment of almond milk shifted the particle size distribution towards right and increased particle size by five- to six-fold. Micrographs confirmed that both the treatments increased particle size due to aggregation of macromolecules. Pressure treatment produced relatively more and larger aggregates than those produced by heat treated samples. The apparent aggregation rate constant for 450 MPa and 600 MPa processed samples were k450MPa,30°C  = 0.0058 s(-1) and k600MPa,30°C  = 0.0095 s(-1) respectively.

CONCLUSIONS

This study showed that dispersive and aggregative properties of high pressure and heat-treated almond milk were different due to differences in protein denaturation, particles coagulation and aggregates morphological characteristics. Knowledge gained from the study will help food processors to formulate novel plant-based beverages treated with high pressure. © 2015 Society of Chemical Industry.

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.