Energy landscape of a small peptide revealed by dihedral angle principal component analysis

Principal component analysis on a torus: Theory and application to protein dynamics

A dimensionality reduction method for high-dimensional circular data is developed, which is based on a principal component analysis (PCA) of data points on a torus. Adopting a geometrical view of PCA, various distance measures on a torus are introduced and the associated problem of projecting data onto the principal subspaces is discussed. The main idea is that the (periodicity-induced) projection error can be minimized by transforming the data such that the maximal gap of the sampling is shifted to the periodic boundary. In a second step, the covariance matrix and its eigendecomposition can be computed in a standard manner. Adopting molecular dynamics simulations of two well-established biomolecular systems (Aib₉ and villin headpiece), the potential of the method to analyze the dynamics of backbone dihedral angles is demonstrated. The new approach allows for a robust and well-defined construction of metastable states and provides low-dimensional reaction coordinates that accurately describe the free energy landscape. Moreover, it offers a direct interpretation of covariances and principal components in terms of the angular variables. Apart from its application to PCA, the method of maximal gap shifting is general and can be applied to any other dimensionality reduction method for circular data.

A Probabilistic Framework for Constructing Temporal Relations in Replica Exchange Molecular Trajectories

Knowledge of the structure and dynamics of biomolecules is essential for elucidating the underlying mechanisms of biological processes. Given the stochastic nature of many biological processes, like protein unfolding, it is almost impossible that two independent simulations will generate the exact same sequence of events, which makes direct analysis of simulations difficult. Statistical models like Markov chains, transition networks, etc. help in shedding some light on the mechanistic nature of such processes by predicting long-time dynamics of these systems from short simulations. However, such methods fall short in analyzing trajectories with partial or no temporal information, for example, replica exchange molecular dynamics or Monte Carlo simulations. In this work, we propose a probabilistic algorithm, borrowing concepts from graph theory and machine learning, to extract reactive pathways from molecular trajectories in the absence of temporal data. A suitable vector representation was chosen to represent each frame in the macromolecular trajectory (as a series of interaction and conformational energies), and dimensionality reduction was performed using principal component analysis (PCA). The trajectory was then clustered using a density-based clustering algorithm, where each cluster represents a metastable state on the potential energy surface (PES) of the biomolecule under study. A graph was created with these clusters as nodes with the edges learned using an iterative expectation maximization algorithm. The most reactive path is conceived as the widest path along this graph. We have tested our method on RNA hairpin unfolding trajectory in aqueous urea solution. Our method makes the understanding of the mechanism of unfolding in the RNA hairpin molecule more tractable. As this method does not rely on temporal data, it can be used to analyze trajectories from Monte Carlo sampling techniques and replica exchange molecular dynamics (REMD).

Protein Folding Free Energy Landscape along the Committor - the Optimal Folding Coordinate

Recent advances in simulation and experiment have led to dramatic increases in the quantity and complexity of produced data, which makes the development of automated analysis tools very important. A powerful approach to analyze dynamics contained in such data sets is to describe/approximate it by diffusion on a free energy landscape - free energy as a function of reaction coordinates (RC). For the description to be quantitatively accurate, RCs should be chosen in an optimal way. Recent theoretical results show that such an optimal RC exists; however, determining it for practical systems is a very difficult unsolved problem. Here we describe a solution to this problem. We describe an adaptive nonparametric approach to accurately determine the optimal RC (the committor) for an equilibrium trajectory of a realistic system. In contrast to alternative approaches, which require a functional form with many parameters to approximate an RC and thus extensive expertise with the system, the suggested approach is nonparametric and can approximate any RC with high accuracy without system specific information. To avoid overfitting for a realistically sampled system, the approach performs RC optimization in an adaptive manner by focusing optimization on less optimized spatiotemporal regions of the RC. The power of the approach is illustrated on a long equilibrium atomistic folding simulation of HP35 protein. We have determined the optimal folding RC - the committor, which was confirmed by passing a stringent committor validation test. It allowed us to determine a first quantitatively accurate protein folding free energy landscape. We have confirmed the recent theoretical results that diffusion on such a free energy profile can be used to compute exactly the equilibrium flux, the mean first passage times, and the mean transition path times between any two points on the profile. We have shown that the mean squared displacement along the optimal RC grows linear with time as for simple diffusion. The free energy profile allowed us to obtain a direct rigorous estimate of the pre-exponential factor for the folding dynamics.

From a Highly Disordered to a Metastable State: Uncovering Insights of α-Synuclein

α-Synuclein (αS) is a major constituent of Lewy bodies, the insoluble aggregates that are the hallmark of one of the most prevalent neurodegenerative disorders, Parkinson's disease (PD). The vast majority of experiments in vitro and in vivo provide extensive evidence that a disordered monomeric form is the predominant state of αS in water solution, and it undergoes a large-scale disorder-to-helix transition upon binding to vesicles of different types. Recently, another form, tetrameric, of αS with a stable helical structure was identified experimentally. It has been shown that a dynamic intracellular population of metastable αS tetramers and monomers coexists normally; and the tetramer plays an essential role in maintaining αS homeostasis. Therefore, it is of interest to know whether the tetramer can serve as a means of preventing or delaying the start of PD. Before answering this very important question, it is, first, necessary to find out, on an atomistic level, a correlation between tetramers and monomers; what mediates tetramer formation and what makes a tetramer stable. We address these questions here by investigating both monomeric and tetrameric forms of αS. In particular, by examining correlations between the motions of the side chains and the main chain, steric parameters along the amino-acid sequence, and one- and two-dimensional free-energy landscapes along the coarse-grained dihedral angles γ and δ and principal components, respectively, in monomeric and tetrameric αS, we were able to shed light on a fundamental relationship between monomers and tetramers, and the key residues involved in mediating formation of a tetramer. Also, the reasons for the stability of tetrameric αS and inability of monomeric αS to fold are elucidated here.

Machine Learning of Biomolecular Reaction Coordinates

We present a systematic approach to reduce the dimensionality of a complex molecular system. Starting with a data set of molecular coordinates (obtained from experiment or simulation) and an associated set of metastable conformational states (obtained from clustering the data), a supervised machine learning model is trained to assign unknown molecular structures to the set of metastable states. In this way, the model learns to determine the features of the molecular coordinates that are most important to discriminate the states. Using a new algorithm that exploits this feature importance via an iterative exclusion principle, we identify the essential internal coordinates (such as specific interatomic distances or dihedral angles) of the system, which are shown to represent versatile reaction coordinates that account for the dynamics of the slow degrees of freedom and explain the mechanism of the underlying processes. Moreover, these coordinates give rise to a free energy landscape that may reveal previously hidden intermediate states of the system.

Designing Well-Structured Cyclic Pentapeptides Based on Sequence-Structure Relationships

Cyclic peptides are a promising class of molecules for unique applications. Unfortunately, cyclic peptide design is severely limited by the difficulty in predicting the conformations they will adopt in solution. In this work, we use explicit-solvent molecular dynamics simulations to design well-structured cyclic peptides by studying their sequence-structure relationships. Critical to our approach is an enhanced sampling method that exploits the essential transitional motions of cyclic peptides to efficiently sample their conformational space. We simulated a range of cyclic pentapeptides from all-glycine to a library of cyclo-(X1X2AAA) peptides to map their conformational space and determine cooperative effects of neighboring residues. By combining the results from all cyclo-(X1X2AAA) peptides, we developed a scoring function to predict the structural preferences for X1-X2 residues within cyclic pentapeptides. Using this scoring function, we designed a cyclic pentapeptide, cyclo-(GNSRV), predicted to be well structured in aqueous solution. Subsequent circular dichroism and NMR spectroscopy revealed that this cyclic pentapeptide is indeed well structured in water, with a nuclear Overhauser effect and J-coupling values consistent with the predicted structure.

Molecular Mechanism, Dynamics, and Energetics of Protein-Mediated Dinucleotide Flipping in a Mismatched DNA: A Computational Study of the RAD4-DNA Complex

DNA damage alters genetic information and adversely affects gene expression pathways leading to various complex genetic disorders and cancers. DNA repair proteins recognize and rectify DNA damage and mismatches with high fidelity. A critical molecular event that occurs during most protein-mediated DNA repair processes is the extrusion of orphaned bases at the damaged site facilitated by specific repairing enzymes. The molecular-level understanding of the mechanism, dynamics, and energetics of base extrusion is necessary to elucidate the molecular basis of protein-mediated DNA damage repair. The present article investigates the molecular mechanism of dinucleotide extrusion in a mismatched DNA (containing a stretch of three contiguous thymidine-thymidine base pairs) facilitated by Radiation sensitive 4 (RAD4), a key DNA repair protein, on an atom-by-atom basis using molecular dynamics (MD) and umbrella-sampling (US) simulations. Using atomistic models of RAD4-free and RAD4-bound mismatched DNA, the free energy profiles associated with extrusion of mismatched partner bases are determined for both systems. The mismatched bases adopted the most stable intrahelical conformation, and their extrusion was unfavorable in RAD4-free mismatched DNA due to the presence of prohibitively high barriers (>12.0 kcal/mol) along the extrusion pathways. Upon binding of RAD4 to the DNA, the global free energy minimum is shifted to the extrahelical state indicating the key role of RAD4-DNA interactions in catalyzing the dinucleotide base extrusion in the DNA-RAD4 complex. The critical residues of RAD4 contributing to the conformational stability of the mismatched bases are identified, and the energetics of insertion of a β-hairpin of RAD4 into the DNA duplex is examined. The conformational energy landscape-based mechanistic insight into RAD4-mediated base extrusion provided here may serve as a useful baseline to understand the molecular basis of xeroderma pigmentosum C (XPC)-mediated DNA damage repair in humans.

Aggregation of PrP106-126 on surfaces of neutral and negatively charged membranes studied by molecular dynamics simulations

Prion diseases are neurodegenerative disorders characterized by the aggregation of an abnormal form of prion protein. The interaction of prion protein and cellular membrane is crucial to elucidate the occurrence and development of prion diseases. Its fragment, residues 106-126, has been proven to maintain the pathological properties of misfolded prion and was used as a model peptide. In this study, explicit solvent molecular dynamics (MD) simulations were carried out to investigate the adsorption, folding and aggregation of PrP106-126 with different sizes (2-peptides, 4-peptides and 6-peptides) on the surface of both pure neutral POPC (1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine) and negatively charged POPC/POPG (1-palmitoyl-2-oleoyl-sn-glycero-3-phosphoglycerol) (3:1) lipids. MD simulation results show that PrP106-126 display strong affinity with POPC/POPG but does not interact with pure POPC. The positively charged and polar residues participating hydrogen bonding with membrane promote the adsorption of PrP106-126. The presence of POPC and POPC/POPG exert limited influence on the secondary structures of PrP106-126 and random coil structures are predominant in all simulation systems. Upon the adsorption on the POPC/POPG surface, the aggregation states of PrP106-126 have been changed and more small oligomers were observed. This work provides insights into the interactions of PrP106-126 and membranes with different compositions in atomic level, which expand our understanding the role membrane plays in the development of prion diseases. This article is part of a Special Issue entitled: Protein Aggregation and Misfolding at the Cell Membrane Interface edited by Ayyalusamy Ramamoorthy.

CAMERRA: An analysis tool for the computation of conformational dynamics by evaluating residue-residue associations

A computational method which extracts the dominant motions from an ensemble of biomolecular conformations via a correlation analysis of residue-residue contacts is presented. The algorithm first renders the structural information into contact matrices, then constructs the collective modes based on the correlated dynamics of a selected set of dynamic contacts. Associated programs can bridge the results for further visualization using graphics software. The aim of this method is to provide an analysis of conformations of biopolymers from the contact viewpoint. It may assist a systematical uncovering of conformational switching mechanisms existing in proteins and biopolymer systems in general by statistical analysis of simulation snapshots. In contrast to conventional correlation analyses of Cartesian coordinates (such as distance covariance analysis and Cartesian principal component analysis), this program also provides an alternative way to locate essential collective motions in general. Herein, we detail the algorithm in a stepwise manner and comment on the importance of the method as applied to decoding allosteric mechanisms. © 2018 Wiley Periodicals, Inc.

Mapping the sequence-structure relationships of simple cyclic hexapeptides

Cyclic peptides are promising protein-protein interaction modulators with high binding affinities and specificities, as well as enhanced stabilities and oral availabilities over linear analogs. Despite their relatively small size and cyclic architecture, it is currently difficult to predict the favored conformation(s) of most classes of cyclic peptides. An improved understanding of the sequence-structure relationships for cyclic peptides will offer an avenue for the rational design of cyclic peptides as possible therapeutics. In this work, we systematically explored the sequence-structure relationships for two cyclic hexapeptide systems using molecular dynamics simulation techniques. Starting with an all-glycine cyclic hexapeptide, cyclo-G₆, we systematically replaced glycine residues with alanines and characterized the structural ensembles of different variants. The same process was repeated with valines to investigate the effects of larger side chains. An analysis of the origin of structure preferences was performed using thermodynamics decomposition and several general observations are reported.

Toward accurately modeling N-methylated cyclic peptides

Cyclic peptides have unique properties and can target protein surfaces specifically and potently. N-Methylation provides a promising way to further optimize the pharmacokinetic and structural profiles of cyclic peptides. The capability to accurately model structures adopted by N-methylated cyclic peptides would facilitate rational design of this interesting and useful class of molecules. We apply molecular dynamics simulations with advanced enhanced sampling methods to efficiently characterize the structural ensembles of N-methylated cyclic peptides, while simultaneously evaluating the overall performance of several simulation force fields. We find that one of the residue-specific force fields, RSFF2, is able to recapitulate experimental structures of the N-methylated cyclic peptide benchmarks tested here when the correct amide isomers are used as initial configurations and enforced during the simulations. Thus, using our simulation approach, it is possible to accurately and efficiently predict the structures of N-methylated cyclic peptides if sufficient information is available to determine the correct amide cis/trans configuration. Moreover, our results suggest that, upon further optimization of RSFF2 to more reliably predict cis/trans isomers, molecular dynamics simulations will be able to de novo predict N-methylated cyclic peptides in the near future, strongly motivating such continued optimization.

Computational insights into the destabilization of α-helical conformations formed by leucine zipper peptides in response to temperature

Recent experiments in our lab (Phys. Chem. Chem. Phys., 2016, 18, 10129-10137) suggested using leucine zipper peptides to enhance the thermosensitivity of liposomes. To understand the mechanisms of temperature-responsive control by the leucine zipper peptide in liposomes, we firstly performed quantum mechanics calculations and implicit-solvent replica exchange molecular dynamics simulations to study the thermo-stability of two leucine zipper peptides, CH₃(CH₂)₄-CO-[VAQLEVK-VAQLESK-VSKLESK-VSSLESK] (termed the capped peptide) and A-[VAQLEVK-VAQLESK-VSKLESK-VSSLESK] (termed the ALA peptide). The analysis of dihedral angle principal components and protein secondary structures was conducted to determine the temperature-dependence conformation transition of the two peptides. Simulation results revealed that our computed transition temperature of the capped peptide is 319.1 K that accords with experimental measurement, 321.1 K. Later, explicit-solvent conventional molecular dynamics simulations were carried out to examine the process of folding and unfolding of the ALA and capped peptides complexed with a lipid bilayer and water in the vicinity of their transition temperatures. A further analysis of conformation and energy of the folded peptides showed that the increase of temperature gives rise to a notable decrease in the number of intra-chain hydrogen bonds and a significant increase in the potential energy of the peptides, thereby reducing the folding stability of the two peptides. As compared to the ALA peptide, a lower transition temperature caused by less intra-chain hydrogen bonds was observed in the capped peptide, which is closer to the temperature of tumor cells. This fact suggests that the capped peptide is more suitable to produce highly sensitive liposomes for the delivery of cancer drugs.

Effects of force fields on the conformational and dynamic properties of amyloid β(1-40) dimer explored by replica exchange molecular dynamics simulations

The conformational space and structural ensembles of amyloid beta (Aβ) peptides and their oligomers in solution are inherently disordered and proven to be challenging to study. Optimum force field selection for molecular dynamics (MD) simulations and the biophysical relevance of results are still unknown. We compared the conformational space of the Aβ(1-40) dimers by 300 ns replica exchange MD simulations at physiological temperature (310 K) using: the AMBER-ff99sb-ILDN, AMBER-ff99sb*-ILDN, AMBER-ff99sb-NMR, and CHARMM22* force fields. Statistical comparisons of simulation results to experimental data and previously published simulations utilizing the CHARMM22* and CHARMM36 force fields were performed. All force fields yield sampled ensembles of conformations with collision cross sectional areas for the dimer that are statistically significantly larger than experimental results. All force fields, with the exception of AMBER-ff99sb-ILDN (8.8 ± 6.4%) and CHARMM36 (2.7 ± 4.2%), tend to overestimate the α-helical content compared to experimental CD (5.3 ± 5.2%). Using the AMBER-ff99sb-NMR force field resulted in the greatest degree of variance (41.3 ± 12.9%). Except for the AMBER-ff99sb-NMR force field, the others tended to under estimate the expected amount of β-sheet and over estimate the amount of turn/bend/random coil conformations. All force fields, with the exception AMBER-ff99sb-NMR, reproduce a theoretically expected β-sheet-turn-β-sheet conformational motif, however, only the CHARMM22* and CHARMM36 force fields yield results compatible with collapse of the central and C-terminal hydrophobic cores from residues 17-21 and 30-36. Although analyses of essential subspace sampling showed only minor variations between force fields, secondary structures of lowest energy conformers are different.

Conformational analysis of spiro-epoxides by principal component analysis of molecular dynamics trajectories

A new procedure for full conformational analyses comprising the statistical analysis of molecular dynamics trajectories was developed and applied. This method included a coordinate space for sampling using molecular dynamics simulations, reduction of dimensionality using tensor decomposition tools, determination of probability distributions in a reduced space, and finally the search for all of the strict extrema points of probability distributions. These extracted extrema points formed an initial guess for geometry optimization and clustering of conformers. A complete conformational space of 1-oxaspiro[2,5]octane and its cis- and trans-4-, 5- and 6-methyl substituted derivatives was also determined. In each case, eight conformers were found with two chair-like conformers predominant at room temperature. It was found that chair-like conformers with an epoxide ring oxygen atom in the pseudo-axial position had less strain, as well as all of their conformers with the methyl substituent in an equatorial position on a cyclohexane moiety.

Conformational heterogeneity in tails of DNA-binding proteins is augmented by proline containing repeats

A cationic terminal extension or tail is a common feature of many DNA-binding proteins. We show that a particular type of tail rich in proline, alanine and lysine belongs to the class of 'flexible disorder' and consists of characteristic pentapeptide repeats. Our designed peptides, (AAKKA)_1-4 and (PAKKA)_1-4, represent the tails of several bacterial DNA-binding proteins. Enhanced conformational sampling of these representative peptides using accelerated molecular dynamic simulations supported by circular dichroism spectroscopy and nuclear magnetic resonance studies demonstrates the role of frequent and interspersed prolines in augmenting conformational heterogeneity of the peptide backbone. Analysis of circular variance of backbone dihedral angles indicates alternating regions of relative rigidity and flexibility along the peptide sequence due to prolines. Preferred placement of lysines in the regions of higher backbone flexibility might improve DNA-binding by conformational selection. Our results could be relevant for rational de novo design of disordered peptides.

Molecular dynamics approach to probe the antigenicity of PagN - an outer membrane protein of Salmonella Typhi

PagN is a highly immunogenic 27-kDa outer membrane adhesin present in Salmonella Typhi. It plays a major role in the pathogenesis of typhoid fever and has emerged as a strong vaccine candidate. In this report, we predict the three-dimensional structure of PagN and describe the conformational dynamics associated with its four extracellular loops based on two 100-ns molecular dynamics simulations at 300 and 310 K. The formation and deformation of the secondary structures on these loops were also investigated during the simulations which revealed loops L1 and L2 to be highly flexible, whereas the relative flexibility of loops L3 and L4 was minimal. Essential dynamics and principal component analysis deciphered more realistic dynamic behaviours of the loops, particularly at 310 K. Moreover, our epitope predictions suggest that the antigenic peptides for B-cell recognition are located within the loops L1 and L2, while those for T-cell recognition are located within the loops L3 and L4. The binding specificities of the antigenic peptides towards specific human MHC-I and MHC-II HLA alleles closely resembled the stability of the loops L3 and L4 inferred from the simulations. Finally, we identified potential antigenic peptides in the flexible (L1 and L2) as well as stable (L3 and L4) regions of PagN for both B- and T-cell recognitions, which can help in developing effective sub-unit vaccines.

Structural review of PPARγ in complex with ligands: Cartesian- and dihedral angle principal component analyses of X-ray crystallographic data

Two decades of research into the ligand-dependent modulation of the activity of the peroxisome proliferator-activated receptor γ (PPARγ) have demonstrated the heterogeneous modes of action of PPARγ ligands, in terms of their interaction surfaces in the ligand-binding pocket, binding stoichiometry and ability to interact with functionally important parts of the receptor, through both direct and allosteric mechanisms. These findings signal the complex mechanistic bases of the distinct biological effects of different classes of PPARγ ligands. Today, the development of PPARγ ligands focuses on partial- and non-agonists as opposed to classical agonists, due to the severe side effects observed with PPARγ classical agonists as therapeutic agents. To aid this development, we performed principal component analyses of the atomic (Cartesian) coordinates (cPCA) and dihedral angles (dPCA) of the structures of human PPARγ from X-ray crystallography, available in the public domain, seeking to reveal ligand-induced trends. In the cPCA, projections of the structures along the principal components (PCs) demonstrated a moderate correlation between cPC1 and structural parameters related to the stabilization of helix 12, which is central to the transcriptional activation by PPARγ classical agonists. Consequently, the presented cPCA mapping of the PPARγ-ligand complexes may guide in silico drug discovery programs seeking to avoid stabilization of helix 12 in their development of partial- and non-agonistic PPARγ ligands. Notably, while the dPCA could identify key regions of dihedral fluctuation in the structural ensemble, the distributions along dPC1 - 2 could not be classified according to the same parameters as the distribution along cPC1. Proteins 2017; 85:1684-1698. © 2017 Wiley Periodicals, Inc.

EGCG inhibits the oligomerization of amyloid beta (16-22) hexamer: Theoretical studies

An extensive replica exchange molecular dynamics (REMD) simulation was performed to investigate the progress patterns of the inhibition of (-)-epigallocatechin-3-gallate (EGCG) on the Aβ_16-22 hexamer. Structural variations of the oligomers without and with EGCG were monitored and analyzed in detail. It has been found that EGCG prevents the formation of Aβ oligomer through two different ways by either accelerating the Aβ oligomerization or reducing the β-content of the hexamer. It also decreases the potential "highly toxic" conformations of Aβ oligomer, which is related to the conformations having high order β-sheet sizes. Both electrostatic and van der Waals interaction energies are found to be involved to the binding process. Computed results using quantum chemical methods show that the π-π stacking is a critical factor of the interaction between EGCG and the peptides. As a result, the binding free energy of the EGCG to the Aβ peptides is slightly larger than that of the curcumin.

Nano-assembly of amyloid β peptide: role of the hairpin fold

Structural investigations have revealed that β hairpin structures are common features in amyloid fibrils, suggesting that these motifs play an important role in amyloid assembly. To test this hypothesis, we characterized the effect of the hairpin fold on the aggregation process using a model β hairpin structure, consisting of two Aβ(14–23) monomers connected by a turn forming YNGK peptide. AFM studies of the assembled aggregates revealed that the hairpin forms spherical structures whereas linear Aβ(14–23) monomers form fibrils. Additionally, an equimolar mixture of the monomer and the hairpin assembles into non-fibrillar aggregates, demonstrating that the hairpin fold dramatically changes the morphology of assembled amyloid aggregates. To understand the molecular mechanism underlying the role of the hairpin fold on amyloid assembly, we performed single-molecule probing experiments to measure interactions between hairpin and monomer and two hairpin complexes. The studies reveal that the stability of hairpin-monomer complexes is much higher than hairpin-hairpin complexes. Molecular dynamics simulations revealed a novel intercalated complex for the hairpin and monomer and Monte Carlo modeling further demonstrated that such nano-assemblies have elevated stability compared with stability of the dimer formed by Aβ(14–23) hairpin. The role of such folding on the amyloid assembly is also discussed.

Systematic Dissociation Pathway Searches Guided by Principal Component Modes

We introduce a novel method, Pathway Search guided by Internal Motions (PSIM), that efficiently finds molecular dissociation pathways of a ligand-receptor system with guidance from principal component (PC) modes obtained from molecular dynamics (MD) simulations. Modeling ligand-receptor dissociation pathways can provide insights into molecular recognition and has practical applications, including understanding kinetic mechanisms and barriers to binding/unbinding as well as design of drugs with desired kinetic properties. PSIM uses PC modes in multilayer internal coordinates to identify natural molecular motions that guide the search for conformational switches and unbinding pathways. The new multilayer internal coordinates overcome problems with Cartesian and classical internal coordinates that fail to smoothly present dihedral rotation or generate nonphysical distortions. We used HIV-1 protease, which has large-scale flap motions, as an example protein to demonstrate use of the multilayer internal coordinates. We provide examples of algorithms and implementation of PSIM with alanine dipeptide and chemical host-guest systems, 2-naphthyl ethanol-β-cyclodextrin and tetramethylammonium-cryptophane complexes. Tetramethylammonium-cryptophane has slow binding/unbinding kinetics. Its residence time, the length to dissociate tetramethylammonium from the host, is ∼14 s from experiments, and PSIM revealed 4 dissociation pathways in approximately 150 CPU h. We also searched the releasing pathways for the product glyceraldehyde-3-phosphate from tryptophan synthase, and one complete dissociation pathway was constructed after running multiple search iterations in approximately 300 CPU h. With guidance by internal PC modes from MD simulations, the PSIM method has advantages over simulation-based methods to search for dissociation pathways of molecular systems with slow noncovalent kinetic behavior.

Imaging Metastable States and Transitions in Proteins by Trajectory Map

It has been a long-standing and intriguing issue to develop robust methods to identify metastable states and interstate transitions from simulations or experimental data to understand the functional conformational changes of proteins. It is usually hard to define the complicated boundaries of the states in the conformational space using most of the existing methods, and they often lead to parameter-sensitive results. Here, we present a new approach, visualized Trajectory Map (vTM), to identify the metastable states and the rare interstate transitions, by considering both the conformational similarity and the temporal successiveness of conformations. The vTM is able to give a nonambiguous description of slow dynamics. The case study of a β-hairpin peptide shows that the vTM can reveal the states and transitions from all-atom MD trajectory data even when a single observable (i.e, one-dimensional reaction coordinate) is used. We also use the vTM to refine the folding/unfolding mechanism of HP35 in explicit water by analyzing a 125 μs all-atom MD trajectory and obtain folding/unfolding rates of about 1/μs, which are in good agreement with the experimental values.

Membrane localization and dynamics of geranylgeranylated Rab5 hypervariable region

The small GTPase Rab5 is a key regulator of endosomal trafficking processes and a marker for the early endosome. The C-terminal hypervariable region (HVR) of Rab5 is post-translationally modified at residues Cys²¹² and Cys²¹³ to accommodate two geranylgeranyl anchors (C20 carbon chain length) in order to associate Rab5 with the membrane. The structural role of the HVR regarding protein-early endosome membrane recruitment is not resolved due to its high degree of flexibility and lack of crystallographic information. Here, full-atomistic and coarse-grained molecular dynamics simulations of the truncated Rab5 HVR^206-215 in three model membranes of increasing complexity (pure phospholipid bilayer, ternary membrane with cholesterol, six-component early endosome) were performed. Specific electrostatic interactions between the HVR^206-215 Arg²⁰⁹ residue and the phosphate group of the inositol ring of PI(3)P were detected. This shows that PI(3)P acts as a first contact site of protein recruitment to the early endosome. The free energy change of HVR^206-215 extraction from the bilayer was largest for the physiological negatively charged membrane. 5μs coarse-grained simulations revealed an active recruitment of PI(3)P to the HVR^206-215 supporting the formation of Rab5- and PI(3)P enriched signaling platforms.

Structural properties of amyloid β(1-40) dimer explored by replica exchange molecular dynamics simulations

Replica exchange molecular dynamics simulations (300 ns) were used to study the dimerization of amyloid β(1-40) (Aβ(1-40)) polypeptide. Configurational entropy calculations revealed that at physiological temperature (310 K, 37°C) dynamic dimers are formed by randomly docked monomers. Free energy of binding of the two chains to each other was -93.56 ± 6.341 kJ mol^-1 . Prevalence of random coil conformations was found for both chains with the exceptions of increased β-sheet content from residues 16-21 and 29-32 of chain A and residues 15-21 and 30-33 of chain B with β-turn/β-bend conformations in both chains from residues 1-16, 21-29 of chain A, 1-16, and 21-29 of chain B. There is a mixed β-turn/β-sheet region from residues 33-38 of both chains. Analysis of intra- and interchain residue distances shows that, although the individual chains are highly flexible, the dimer system stays in a loosely packed antiparallel β-sheet configuration with contacts between residues 17-21 of chain A with residues 17-21 and 31-36 of chain B as well as residues 31-36 of chain A with residues 17-21 and 31-36 of chain B. Based on dihedral principal component analysis, the antiparallel β-sheet-loop-β-sheet conformational motif is favored for many low energy sampled conformations. Our results show that Aβ(1-40) can form dynamic dimers in aqueous solution that have significant conformational flexibility and are stabilized by collapse of the central and C-terminal hydrophobic cores with the expected β-sheet-loop-β-sheet conformational motif. Proteins 2017; 85:1024-1045. © 2017 Wiley Periodicals, Inc.

Optimized parameter selection reveals trends in Markov state models for protein folding

As molecular dynamics simulations access increasingly longer time scales, complementary advances in the analysis of biomolecular time-series data are necessary. Markov state models offer a powerful framework for this analysis by describing a system's states and the transitions between them. A recently established variational theorem for Markov state models now enables modelers to systematically determine the best way to describe a system's dynamics. In the context of the variational theorem, we analyze ultra-long folding simulations for a canonical set of twelve proteins [K. Lindorff-Larsen et al., Science 334, 517 (2011)] by creating and evaluating many types of Markov state models. We present a set of guidelines for constructing Markov state models of protein folding; namely, we recommend the use of cross-validation and a kinetically motivated dimensionality reduction step for improved descriptions of folding dynamics. We also warn that precise kinetics predictions rely on the features chosen to describe the system and pose the description of kinetic uncertainty across ensembles of models as an open issue.

Protein Structure Classification and Loop Modeling Using Multiple Ramachandran Distributions

Recently, the study of protein structures using angular representations has attracted much attention among structural biologists. The main challenge is how to efficiently model the continuous conformational space of the protein structures based on the differences and similarities between different Ramachandran plots. Despite the presence of statistical methods for modeling angular data of proteins, there is still a substantial need for more sophisticated and faster statistical tools to model the large-scale circular datasets. To address this need, we have developed a nonparametric method for collective estimation of multiple bivariate density functions for a collection of populations of protein backbone angles. The proposed method takes into account the circular nature of the angular data using trigonometric spline which is more efficient compared to existing methods. This collective density estimation approach is widely applicable when there is a need to estimate multiple density functions from different populations with common features. Moreover, the coefficients of adaptive basis expansion for the fitted densities provide a low-dimensional representation that is useful for visualization, clustering, and classification of the densities. The proposed method provides a novel and unique perspective to two important and challenging problems in protein structure research: structure-based protein classification and angular-sampling-based protein loop structure prediction.

Elucidating Important Sites and the Mechanism for Amyloid Fibril Formation by Coarse-Grained Molecular Dynamics

Fibrils formed by the β-amyloid (Aβ) peptide play a central role in the development of Alzheimer's disease. In this study, the principles governing their growth and stability are investigated by analyzing canonical and replica-exchange molecular dynamics trajectories of Aβ_(9-40) fibrils. In particular, an unstructured monomer was allowed to interact freely with an Aβ fibril template. Trajectories were generated with the coarse-grained united-residue force field, and one- and two-dimensional free-energy landscapes (FELs) along the backbone virtual-bond angle θ and backbone virtual-bond-dihedral angle γ of each residue and principal components, respectively, were analyzed. Also, thermal unbinding (unfolding) of an Aβ peptide from the fibril template was investigated. These analyses enable us to illustrate the entire process of Aβ fibril elongation and to elucidate the key residues involved in it. Several different pathways were identified during the search for the fibril conformation by the monomer, which finally follows a dock-lock mechanism with two distinct locking stages. However, it was found that the correct binding, with native hydrogen bonds, of the free monomer to the fibril template at both stages is crucial for fibril elongation. In other words, if the monomer is incorrectly bound (with nonnative hydrogen bonds) to the fibril template during the first "docking" stage, it can remain attached to it for a long time before it dissociates and either attempts a different binding or allows another monomer to bind. This finding is consistent with an experimentally observed "stop-and-go" mechanism of fibril growth.

Lid domain plasticity and lipid flexibility modulate enzyme specificity in human monoacylglycerol lipase

Human monoacylglycerol lipase (MAGL) is a membrane-interacting enzyme that generates pro-inflammatory signaling molecules. For this reason, MAGL inhibition is a promising strategy to treat pain, cancer, and neuroinflammatory diseases. MAGL can hydrolyze monoacylglycerols bearing an acyl chain of different lengths and degrees of unsaturation, cleaving primarily the endocannabinoid 2-arachidonoylglycerol. Importantly, the enzymatic binding site of MAGL is confined by a 75-amino-acid-long, flexible cap domain, named 'lid domain', which is structurally similar to that found in several other lipases. However, it is unclear how lid domain plasticity affects catalysis in MAGL. By integrating extensive molecular dynamics simulations and free-energy calculations with mutagenesis and kinetic experiments, we here define a lid-domain-mediated mechanism for substrate selection and binding in MAGL catalysis. In particular, we clarify the key role of Phe159 and Ile179, two conserved residues within the lid domain, in regulating substrate specificity in MAGL. We conclude by proposing that other structurally related lipases may share this lid-domain-mediated mechanism for substrate specificity.

Replica exchange molecular dynamics study of the amyloid beta (11–40) trimer penetrating a membrane

The transmembrane Aβ_11–40 trimer is investigated for the first time using REMD and FEP.

The sensitivity of folding free energy landscapes of trpzips to mutations in the hydrophobic core

The sensitivity of the stability of folded states and free energy landscapes to the differences in the hydrophobic content of the core residues has been studied for the set of 16-residue trpzips, namely, Trpzip4, Trpzip5 and Trpzip6.

Molecular dynamics approach to understand the denaturing effect of a millimolar concentration of dodine on a λ-repressor and counteraction by trehalose

Here, we use a molecular dynamics approach to calculate the spatial distribution function of the ternary water–dodine–trehalose (1.0 M) system.

Replica exchange molecular dynamics study of the truncated amyloid beta (11–40) trimer in solution

The structure of the 3Aβ_11–40 oligomer is determined for the first time using T-REMD simulations.

Self-assembly of the full-length amyloid Aβ42 protein in dimers

The self-assembly of amyloid (Aβ) proteins into nano-aggregates is a hallmark of Alzheimer's disease (AD) development, yet the mechanism of how disordered monomers assemble into aggregates remains elusive. Here, we applied long-time molecular dynamics simulations to fully characterize the assembly of Aβ42 monomers into dimers. Monomers undergo conformational changes during their interaction, but the resulting dimer structures do not resemble those found in fibril structures. To identify natural conformations of dimers among a set of simulated ones, validation approaches were developed and applied, and a subset of dimer conformations were characterized. These dimers do not contain long β-strands that are usually found in fibrils. The dimers are stabilized primarily by interactions within the central hydrophobic regions and the C-terminal regions, with a contribution from local hydrogen bonding. The dimers are dynamic, as evidenced by the existence of a set of conformations and by the quantitative analyses of the dimer dissociation process.

Exploring Flexibility of Progesterone Receptor Ligand Binding Domain Using Molecular Dynamics

Progesterone receptor (PR), a member of nuclear receptor (NR) superfamily, plays a vital role for female reproductive tissue development, differentiation and maintenance. PR ligand, such as progesterone, induces conformation changes in PR ligand binding domain (LBD), thus mediates subsequent gene regulation cascades. PR LBD may adopt different conformations upon an agonist or an antagonist binding. These different conformations would trigger distinct transcription events. Therefore, the dynamics of PR LBD would be of general interest to biologists for a deep understanding of its structure-function relationship. However, no apo-form (non-ligand bound) of PR LBD model has been proposed either by experiments or computational methods so far. In this study, we explored the structural dynamics of PR LBD using molecular dynamics simulations and advanced sampling tools in both ligand-bound and the apo-forms. Resolved by the simulation study, helix 11, helix 12 and loop 895–908 (the loop between these two helices) are quite flexible in antagonistic conformation. Several residues, such as Arg899 and Glu723, could form salt-bridging interaction between helix 11 and helix 3, and are important for the PR LBD dynamics. And we also propose that helix 12 in apo-form PR LBD, not like other NR LBDs, such as human estrogen receptor α (ERα) LBD, may not adopt a totally extended conformation. With the aid of umbrella sampling and metadynamics simulations, several stable conformations of apo-form PR LBD have been sampled, which may work as critical structural models for further large scale virtual screening study to discover novel PR ligands for therapeutic application.

Distributions of experimental protein structures on coarse-grained free energy landscapes

Predicting conformational changes of proteins is needed in order to fully comprehend functional mechanisms. With the large number of available structures in sets of related proteins, it is now possible to directly visualize the clusters of conformations and their conformational transitions through the use of principal component analysis. The most striking observation about the distributions of the structures along the principal components is their highly non-uniform distributions. In this work, we use principal component analysis of experimental structures of 50 diverse proteins to extract the most important directions of their motions, sample structures along these directions, and estimate their free energy landscapes by combining knowledge-based potentials and entropy computed from elastic network models. When these resulting motions are visualized upon their coarse-grained free energy landscapes, the basis for conformational pathways becomes readily apparent. Using three well-studied proteins, T4 lysozyme, serum albumin, and sarco-endoplasmic reticular Ca(2+) adenosine triphosphatase (SERCA), as examples, we show that such free energy landscapes of conformational changes provide meaningful insights into the functional dynamics and suggest transition pathways between different conformational states. As a further example, we also show that Monte Carlo simulations on the coarse-grained landscape of HIV-1 protease can directly yield pathways for force-driven conformational changes.

Succession of Alkane Conformational Motifs Bound within Hydrophobic Supramolecular Capsular Assemblies

n-Alkane encapsulation experiments within dimeric octa-acid cavitand capsules in water reveal a succession of packing motifs from extended, to helical, to hairpin, to spinning top structures with increasing chain length. Here, we report a molecular simulation study of alkane conformational preferences within these host-guest assemblies to uncover the factors stabilizing distinct conformers. The simulated alkane conformers follow the trends inferred from ¹H NMR experiments, while guest proton chemical shifts evaluated from Gauge Invariant Atomic Orbital calculations provide further evidence our simulations capture guest packing within these assemblies. Analysis of chain length and dihedral distributions indicates that packing under confinement to minimize nonpolar guest and host interior contact with water largely drives the transitions. Mean intramolecular distance maps and transfer free energy differences suggest the extended and helical motifs are members of a larger family of linear guest structures, for which the guest gauche population increases with increasing chain length to accommodate the chains within the complex. Breaks observed between the helical/hairpin and hairpin/spinning top motifs, on the other hand, indicate the hairpin and spinning top conformations are distinct from the linear family. Our results represent the first bridging of empirical and simulation data for flexible guests encapsulated within confined nanospaces, and constitute an effective strategy by which guest packing motifs within artificial or natural compartments can be rationalized and/or predicted a priori.