Dihedral angle principal component analysis of molecular dynamics simulations

Tripleurin XIIc: Peptide Folding Dynamics in Aqueous and Hydrophobic Environment Mimic Using Accelerated Molecular Dynamics

Peptaibols are a special class of fungal peptides with an acetylated N-terminus and a C-terminal 1,2-amino alcohol along with non-standard amino acid residues. New peptaibols named tripleurins were recently identified from a strain of the filamentous fungal species Trichoderma pleuroti, which is known to cause green mould disease on cultivated oyster mushrooms. To understand the mode of action of these peptaibols, the three-dimensional structure of tripleurin (TPN) XIIc, an 18-mer peptide, was elucidated using an enhanced sampling method, accelerated MD, in water and chloroform solvents. Non-standard residues were parameterized by the Restrained Electrostatic Potential (RESP) charge fitting method. The dihedral distribution indicated towards a right-handed helical formation for TPN XIIc in both solvents. Dihedral angle based principal component analysis revealed a propensity for a slightly bent, helical folded conformation in water solvent, while two distinct conformations were revealed in chloroform: One that folds into highly bent helical structure that resembles a beta-hairpin and another with an almost straight peptide backbone appearing as a rare energy barrier crossing event. The hinge-like movement of the terminals was also observed and is speculated to be functionally relevant. The convergence and efficient sampling is addressed using Cartesian PCA and Kullback-Leibler divergence methods.

Folding Molecular Dynamics Simulation of a gp41-Derived Peptide Reconcile Divergent Structure Determinations

T-20 peptide is the first FDA-approved fusion inhibitor against AIDS/HIV-1 gp41 protein. Part of it, the gp41[659-671] peptide, that contains the complete epitope for the neutralizing 2F5 monoclonal antibody, has been found experimentally in a number of divergent structures. Herein, we attempt to reconcile them by using unbiased large-scale all-atom molecular dynamics folding simulations. We show that our approach can successfully capture the peptide's heterogeneity and reach each and every experimentally determined conformation in sub-angstrom accuracy, whilst preserving the peptide's disordered nature. Our analysis also unveils that the minor refinements within the AMBER99SB family of force fields can lead to appreciable differences in the predicted conformational stability arising from subtle differences in the helical- and β-region of the Ramachandran plot. Our work underlines the contribution of molecular dynamics simulation in structurally characterizing pharmacologically important peptides of ambiguous structure.

Perspective: Identification of collective variables and metastable states of protein dynamics

The statistical analysis of molecular dynamics simulations requires dimensionality reduction techniques, which yield a low-dimensional set of collective variables (CVs) {x _i } = x that in some sense describe the essential dynamics of the system. Considering the distribution P( x ) of the CVs, the primal goal of a statistical analysis is to detect the characteristic features of P( x ), in particular, its maxima and their connection paths. This is because these features characterize the low-energy regions and the energy barriers of the corresponding free energy landscape ΔG( x ) = -k _B T ln P( x ), and therefore amount to the metastable states and transition regions of the system. In this perspective, we outline a systematic strategy to identify CVs and metastable states, which subsequently can be employed to construct a Langevin or a Markov state model of the dynamics. In particular, we account for the still limited sampling typically achieved by molecular dynamics simulations, which in practice seriously limits the applicability of theories (e.g., assuming ergodicity) and black-box software tools (e.g., using redundant input coordinates). We show that it is essential to use internal (rather than Cartesian) input coordinates, employ dimensionality reduction methods that avoid rescaling errors (such as principal component analysis), and perform density based (rather than k-means-type) clustering. Finally, we briefly discuss a machine learning approach to dimensionality reduction, which highlights the essential internal coordinates of a system and may reveal hidden reaction mechanisms.

Collective variable discovery and enhanced sampling using autoencoders: Innovations in network architecture and error function design

Auto-associative neural networks ("autoencoders") present a powerful nonlinear dimensionality reduction technique to mine data-driven collective variables from molecular simulation trajectories. This technique furnishes explicit and differentiable expressions for the nonlinear collective variables, making it ideally suited for integration with enhanced sampling techniques for accelerated exploration of configurational space. In this work, we describe a number of sophistications of the neural network architectures to improve and generalize the process of interleaved collective variable discovery and enhanced sampling. We employ circular network nodes to accommodate periodicities in the collective variables, hierarchical network architectures to rank-order the collective variables, and generalized encoder-decoder architectures to support bespoke error functions for network training to incorporate prior knowledge. We demonstrate our approach in blind collective variable discovery and enhanced sampling of the configurational free energy landscapes of alanine dipeptide and Trp-cage using an open-source plugin developed for the OpenMM molecular simulation package.

Molecular modeling of conformational dynamics and its role in enzyme evolution

With increasing computational power, biomolecular simulations have become an invaluable tool for understanding enzyme mechanisms and the origins of enzyme catalysis. More recently, computational studies have started to focus on understanding how enzyme activity itself evolves, both in terms of enhancing the native or new activities on existing enzyme scaffolds, or completely de novo on previously non-catalytic scaffolds. In this context, both experiment and molecular modeling provided strong evidence for an important role of conformational dynamics in the evolution of enzyme functions. This contribution will present a brief overview of the current state of the art for computationally exploring enzyme conformational dynamics in enzyme evolution, and, using several showcase studies, illustrate the ways molecular modeling can be used to shed light on how enzyme function evolves, at the most fundamental molecular level.

Replica-Based Protein Structure Sampling Methods: Compromising between Explicit and Implicit Solvents

The structure of a protein is often not completely accessible by experiments. In silico, replica exchange molecular dynamics (REMD) is the standard sampling method for predicting the secondary and tertiary structures from the amino acid sequence, but it is computationally very expensive. Two recent adaptations from REMD, temperature intervals with global exchange of replicas (TIGER2) and TIGER2A, have been tested here in implicit and explicit solvents. Additionally, explicit, implicit, and hybrid solvent REMD are compared. On the basis of the hybrid REMD (REMDh) method, we present a new hybrid TIGER2h algorithm for faster structural sampling, while retaining good accuracy. The implementations of REMDh, TIGER2, TIGER2A, and TIGER2h are provided for nanoscale molecular dynamics (NAMD). All the methods were tested with two model peptides of known structure, (AAQAA)₃ and HP7, with helix and sheet motifs, respectively. The TIGER2 methods and REMDh were also applied to the unknown structure of the collagen type I telopeptides, which represent bigger proteins with some degree of disorder. We present simulations covering more than 180 μs and analyze the performance and convergence of the distributions of states between the particular methods by dihedral principal component and secondary structure analysis.

Time-lagged autoencoders: Deep learning of slow collective variables for molecular kinetics

Inspired by the success of deep learning techniques in the physical and chemical sciences, we apply a modification of an autoencoder type deep neural network to the task of dimension reduction of molecular dynamics data. We can show that our time-lagged autoencoder reliably finds low-dimensional embeddings for high-dimensional feature spaces which capture the slow dynamics of the underlying stochastic processes-beyond the capabilities of linear dimension reduction techniques.

Modulation of phase transition of thermosensitive liposomes with leucine zipper-structured lipopeptides

Targeted therapy for cancer requires thermosensitive components in drug carriers for controlled drug release against viral cells. The conformational transition characteristic of leucine zipper-structured lipopeptides is utilized in our lab to modulate the phase transition temperature of liposomes, thus achieving temperature-responsive control. In this study, we computationally examined the conformational transition behaviors of leucine zipper-structured lipopeptides that were modified at the N-terminus by distinct functional groups. The conformational transition temperatures of these lipopeptides were determined by structural analysis of the implicit-solvent replica exchange molecular dynamics simulation trajectories using the dihedral angle principal component analysis and the dictionary of protein secondary structure method. Our calculations revealed that the computed transition temperatures of the lipopeptides are in good agreement with the experimental measurements. The effect of hydrogen bonds on the conformational stability of the lipopeptide dimers was examined in conventional explicit-solvent molecular dynamics simulations. A quantitative correlation of the degree of structural dissociation of the dimers and their binding strength is well described by an exponential fit of the binding free energies to the conformation transition temperatures of the lipopeptides.

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.

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.

Modeling of alcohol oxidase enzyme of Candida boidinii and in silico analysis of competitive binding of proton ionophores and FAD with enzyme

Alcohol oxidase (AOX) is an important flavin adenine dinucleotide (FAD) dependent oxidoreductase, which is responsible for converting methanol into formaldehyde and hydrogen peroxide for the growth of methylotrophic yeast Candida boidinii. Although AOX plays a crucial role in methanol catabolism, the experimental structure of AOX from Candida boidinii has not been elucidated. This study reports the first complete in silico model of AOX from C. boidinii. This paper also reports the AOX structure modeled using the threading approach, followed by structure analysis and molecular dynamics simulation. The modeled structure was compared with the aryl alcohol oxidase structure (a glucose-methanol-choline family member, pdbID: 3fim). A docking study was performed to analyze the interaction between AOX and its cofactor FAD. The AOX modeled structure also exhibited high similarity with respect to the FAD binding sites, which are the substrate binding sites as seen with 3fim. It was observed that the adenosine part of FAD was deeply buried inside AOX while the isoalloxazine ring stuck to the surface. This paper reports the interaction of selective proton ionophores (CCCP and DNP) with AOX and also reports their binding sites. These proton ionophores showed competitive binding with FAD. The occupancy of the FAD binding sites by the proton ionophore may lead to blocking of the entry of FAD and thereby disruption of AOX import into peroxisomes.

Atomistic molecular dynamics simulations of bioactive engrailed 1 interference peptides (EN1-iPeps)

The neural-specific transcription factor Engrailed 1 - is overexpressed in basal-like breast tumours. Synthetic interference peptides - comprising a cell-penetrating peptide/nuclear localisation sequence and the Engrailed 1-specific sequence from the N-terminus have been engineered to produce a strong apoptotic response in tumour cells overexpressing EN1, with no toxicity to normal or non Engrailed 1-expressing cells. Here scaled molecular dynamics simulations were used to study the conformational dynamics of these interference peptides in aqueous solution to characterise their structure and dynamics. Transitions from disordered to α-helical conformation, stabilised by hydrogen bonds and proline-aromatic interactions, were observed throughout the simulations. The backbone of the wild-type peptide folds to a similar conformation as that found in ternary complexes of anterior Hox proteins with conserved hexapeptide motifs important for recognition of pre-B-cell leukemia Homeobox 1, indicating that the motif may possess an intrinsic preference for helical structure. The predicted NMR chemical shifts of these peptides are consistent with the Hox hexapeptides in solution and Engrailed 2 NMR data. These findings highlight the importance of aromatic residues in determining the structure of Engrailed 1 interference peptides, shedding light on the rational design strategy of molecules that could be adopted to inhibit other transcription factors overexpressed in other cancer types, potentially including other transcription factor families that require highly conserved and cooperative protein-protein partnerships for biological activity.

Structure of APP-C991-99 and implications for role of extra-membrane domains in function and oligomerization

The 99 amino acid C-terminal fragment of Amyloid Precursor Protein APP-C99 (C99) is cleaved by γ-secretase to form Aβ peptide, which plays a critical role in the etiology of Alzheimer's Disease (AD). The structure of C99 consists of a single transmembrane domain flanked by intra and intercellular domains. While the structure of the transmembrane domain has been well characterized, little is known about the structure of the flanking domains and their role in C99 processing by γ-secretase. To gain insight into the structure of full-length C99, REMD simulations were performed for monomeric C99 in model membranes of varying thickness. We find equilibrium ensembles of C99 from simulation agree with experimentally-inferred residue insertion depths and protein backbone chemical shifts. In thin membranes, the transmembrane domain structure is correlated with extra-membrane structural states and the extra-membrane domain structural states become less correlated to each other. Mean and variance of the transmembrane and G₃₇G₃₈ hinge angles are found to increase with thinning membrane. The N-terminus of C99 forms β-strands that may seed aggregation of Aβ on the membrane surface, promoting amyloid formation. In thicker membranes the N-terminus forms α-helices that interact with the nicastrin domain of γ-secretase. The C-terminus of C99 becomes more α-helical as the membrane thickens, forming structures that may be suitable for binding by cytoplasmic proteins, while C-terminal residues essential to cytotoxic function become α-helical as the membrane thins. The heterogeneous but discrete extra-membrane domain states analyzed here open the path to new investigations of the role of C99 structure and membrane in amyloidogenesis. This article is part of a Special Issue entitled: Protein Aggregation and Misfolding at the Cell Membrane Interface edited by Ayyalusamy Ramamoorthy.

Molecular dynamics simulation of metal free structure of Lmb, a laminin-binding adhesin of Streptococcus agalactiae: metal removal and its structural implications

Metal-binding receptors are one of the extracellular components of ATP-binding cassette transporters that are essential for regulation of metal homeostasis in bacteria. Laminin-binding adhesin (Lmb) of Streptococcus agalactiae falls under this class of solute binding proteins. It binds to zinc with a high affinity. Crystal structure of Lmb solved previously by our group reveals that the zinc is tetrahedrally coordinated by three histidines and a glutamate at the interdomain cleft. Lmb contains a long disordered loop close to the metal-binding site whose precise function is unknown. Several experimental attempts to produce apo-Lmb failed and this prompted us to carry out in silico studies to analyse the structural importance of the metal in Lmb. Here, we present the results of the molecular dynamics (MD) simulation studies of native, apo-(metal removed) and the long loop truncated Lmb models along with a homologous protein, TroA from Treponema pallidum that was taken up for validating the MD results of Lmb. Absence of a metal results in significant structural changes in Lmb, particularly at the metal-binding pocket and with the long loop, although the overall fold is retained. This study thus revealed that the Lmb can exist in different conformational states with subtle differences in the overall fold based on the presence or absence of the metal. This could be functionally important for a putative metal uptake and release and also for the adhesive function of Lmb in recognizing laminin, which contains a high number of zinc finger motifs.

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.

Challenges and Opportunities in Drug Discovery of Biased Ligands

The observation of biased agonism in G protein-coupled receptors (GPCRs) has provided new approaches for the development of more efficacious and safer drugs. However, in order to rationally design biased drugs, one must understand the molecular basis of this phenomenon. Computational approaches can help in exploring the conformational universe of GPCRs and detecting conformational states with relevance for distinct functional outcomes. This information is extremely valuable for the development of new therapeutic agents that promote desired conformational receptor states and responses while avoiding the ones leading to undesired side-effects.This book chapter intends to introduce the reader to powerful computational approaches for sampling the conformational space of these receptors, focusing first on molecular dynamics and the analysis of the produced data through methods such as dimensionality reduction, Markov State Models and adaptive sampling. Then, we show how to seek for compounds that target distinct conformational states via docking and virtual screening. In addition, we describe how to detect receptor-ligand interactions that drive signaling bias and comment current challenges and opportunities of presented methods.

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.

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.

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.

Dynamical footprint of cross-reactivity in a human autoimmune T-cell receptor

The present work focuses on the dynamical aspects of cross-reactivity between myelin based protein (MBP) self-peptide and two microbial peptides (UL15, PMM) for Hy.1B11 T-cell receptor (TCR). This same TCR was isolated from a patient suffering from multiple sclerosis (MS). The study aims at highlighting the chemical interactions underlying recognition mechanisms between TCR and the peptides presented by Major Histocompatibility Complex (MHC) proteins, which form a crucial component in adaptive immune response against foreign antigens. Since the ability of a TCR to recognize different peptide antigens presented by MHC depends on its cross-reactivity, we used molecular dynamics methods to obtain atomistic detail on TCR-peptide-MHC complexes. Our results show how the dynamical basis of Hy.1B11 TCR’s cross-reactivity is rooted in a similar bridging interaction pattern across the TCR-peptide-MHC interface. Our simulations confirm the importance of TCR CDR3α E98 residue interaction with MHC and a predominant role of P6 peptide residue in MHC binding affinity. Altogether, our study provides energetic and dynamical insights into factors governing peptide recognition by the cross-reactive Hy.1B11 TCR, found in MS patient.

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.

Insights into the Structural Dynamics of Nucleocytoplasmic Transport of tRNA by Exportin-t

Exportin-t (Xpot) transports mature 5'- and 3'-end processed tRNA from the nucleus to the cytoplasm by associating with a small G-protein Ran (RAs-related nuclear protein), in the nucleus. The release of tRNA in cytoplasm involves RanGTP hydrolysis. Despite the availability of crystal structures of nuclear and cytosolic forms of Xpot, the molecular details regarding the sequential events leading to tRNA release and subsequent conformational changes occurring in Xpot remain unknown. We have performed a combination of classical all-atom and accelerated molecular dynamics simulations on a set of complexes involving Xpot to study a range of features including conformational flexibility of free and cargo-bound Xpot and functionally critical contacts between Xpot and its cargo. The systems investigated include free Xpot and its different complexes, bound either to Ran (GTP/GDP) or tRNA or both. This approach provided a statistically reliable estimate of structural dynamics of Xpot after cargo release. The mechanistic basis for Xpot opening after cargo release has been explained in terms of dynamic structural hinges, about which neighboring region could be displaced to facilitate the nuclear to cytosolic state transition. Post-RanGTP hydrolysis, a cascade of events including local conformational change in RanGTP and loss of critical contacts at Xpot/tRNA interface suggest factors responsible for eventual release of tRNA. The level of flexibility in different Xpot complexes varied depending on the arrangement of individual HEAT repeats. Current study provides one of the most comprehensive and robust analysis carried out on this protein using molecular dynamics schemes.

Hierarchical Time-Lagged Independent Component Analysis: Computing Slow Modes and Reaction Coordinates for Large Molecular Systems

Analysis of molecular dynamics, for example using Markov models, often requires the identification of order parameters that are good indicators of the rare events, i.e. good reaction coordinates. Recently, it has been shown that the time-lagged independent component analysis (TICA) finds the linear combinations of input coordinates that optimally represent the slow kinetic modes and may serve in order to define reaction coordinates between the metastable states of the molecular system. A limitation of the method is that both computing time and memory requirements scale with the square of the number of input features. For large protein systems, this exacerbates the use of extensive feature sets such as the distances between all pairs of residues or even heavy atoms. Here we derive a hierarchical TICA (hTICA) method that approximates the full TICA solution by a hierarchical, divide-and-conquer calculation. By using hTICA on distances between heavy atoms we identify previously unknown relaxation processes in the bovine pancreatic trypsin inhibitor.

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.

Mechanisms for allosteric activation of protease DegS by ligand binding and oligomerization as revealed from molecular dynamics simulations

A local perturbation of a protein may lead to functional changes at some distal site, a phenomenon denoted as allostery. Here, we study the allosteric control of a protease using molecular dynamics simulations. The system considered is the bacterial protein DegS which includes a protease domain activated on ligand binding to an adjacent PDZ domain. Starting from crystallographic structures of DegS homo-trimers, we perform simulations of the ligand-free and -bound state of DegS at equilibrium. Considering a single protomer only, the trimeric state was mimicked by applying restraints on the residues in contact with other protomers in the DegS trimer. In addition, the bound state was also simulated without any restraints to mimic the monomer. Our results suggest that not only ligand release but also disassembly of a DegS trimer inhibits proteolytic activity. Considering various observables for structural changes, we infer allosteric pathways from the interface with other protomers to the active site. Moreover, we study how ligand release leads to (i) catalytically relevant changes involving residues 199-201 and (ii) a transition from a stretched to a bent conformation for residues 217-219 (which prohibits proper substrate binding). Finally, based on ligand-induced C_α shifts we identify residues in contact with other protomers in the DegS trimer that likely transduce the perturbation from ligand release from a given protomer to adjacent protomers. These residues likely play a key role in the experimentally known effect of ligand release from a protomer on the proteolytic activity of the other protomers. Proteins 2016; 84:1690-1705. © 2016 Wiley Periodicals, Inc.

In silico Exploration of the Conformational Universe of GPCRs

The structural plasticity of G protein coupled receptors (GPCRs) leads to a conformational universe going from inactive to active receptor states with several intermediate states. Many of them have not been captured yet and their role for GPCR activation is not well understood. The study of this conformational space and the transition dynamics between different receptor populations is a major challenge in molecular biophysics. The rational design of effector molecules that target such receptor populations allows fine-tuning receptor signalling with higher specificity to produce drugs with safer therapeutic profiles. In this minireview, we outline highly conserved receptor regions which are considered determinant for the establishment of distinct receptor states. We then discuss in-silico approaches such as dimensionality reduction methods and Markov State Models to explore the GPCR conformational universe and exploit the obtained conformations through structure-based drug design.

Free-energy landscape of a hyperstable RNA tetraloop

We report the characterization of the energy landscape and the folding/unfolding thermodynamics of a hyperstable RNA tetraloop obtained through high-performance molecular dynamics simulations at microsecond timescales. Sampling of the configurational landscape is conducted using temperature replica exchange molecular dynamics over three isochores at high, ambient, and negative pressures to determine the thermodynamic stability and the free-energy landscape of the tetraloop. The simulations reveal reversible folding/unfolding transitions of the tetraloop into the canonical A-RNA conformation and the presence of two alternative configurations, including a left-handed Z-RNA conformation and a compact purine Triplet. Increasing hydrostatic pressure shows a stabilizing effect on the A-RNA conformation and a destabilization of the left-handed Z-RNA. Our results provide a comprehensive description of the folded free-energy landscape of a hyperstable RNA tetraloop and highlight the significant advances of all-atom molecular dynamics in describing the unbiased folding of a simple RNA secondary structure motif.

Cations Stiffen Actin Filaments by Adhering a Key Structural Element to Adjacent Subunits

Ions regulate the assembly and mechanical properties of actin filaments. Recent work using structural bioinformatics and site-specific mutagenesis favors the existence of two discrete and specific divalent cation binding sites on actin filaments, positioned in the long axis between actin subunits. Cation binding at one site drives polymerization, while the other modulates filament stiffness and plays a role in filament severing by the regulatory protein, cofilin. Existing structural methods have not been able to resolve filament-associated cations, and so in this work we turn to molecular dynamics simulations to suggest a candidate binding pocket geometry for each site and to elucidate the mechanism by which occupancy of the “stiffness site” affects filament mechanical properties. Incorporating a magnesium ion in the “polymerization site” does not seem to require any large-scale change to an actin subunit’s conformation. Binding of a magnesium ion in the “stiffness site” adheres the actin DNase-binding loop (D-loop) to its long-axis neighbor, which increases the filament torsional stiffness and bending persistence length. Our analysis shows that bound D-loops occupy a smaller region of accessible conformational space. Cation occupancy buries key conserved residues of the D-loop, restricting accessibility to regulatory proteins and enzymes that target these amino acids.

Robust Density-Based Clustering To Identify Metastable Conformational States of Proteins

A density-based clustering method is proposed that is deterministic, computationally efficient, and self-consistent in its parameter choice. By calculating a geometric coordinate space density for every point of a given data set, a local free energy is defined. On the basis of these free energy estimates, the frames are lumped into local free energy minima, ultimately forming microstates separated by local free energy barriers. The algorithm is embedded into a complete workflow to robustly generate Markov state models from molecular dynamics trajectories. It consists of (i) preprocessing of the data via principal component analysis in order to reduce the dimensionality of the problem, (ii) proposed density-based clustering to generate microstates, and (iii) dynamical clustering via the most probable path algorithm to construct metastable states. To characterize the resulting state-resolved conformational distribution, dihedral angle content color plots are introduced which identify structural differences of protein states in a concise way. To illustrate the performance of the method, three well-established model problems are adopted: conformational transitions of hepta-alanine, folding of villin headpiece, and functional dynamics of bovine pancreatic trypsin inhibitor.

Molecular Dynamics Analysis of Binding of Kinase Inhibitors to WT EGFR and the T790M Mutant

Epidermal growth factor receptor (EGFR) inhibitors interrupt EGFR-dependent cellular signaling pathways that lead to accelerated tumor growth and proliferation. Mutation of a threonine in the inhibitor binding pocket, known as the "gatekeeper", to methionine (T790M) confers acquired resistance to several EGFR-selective inhibitors. We studied interactions between EGFR inhibitors and the gatekeeper residues of the target protein. Thermodynamic integration (TI) with Amber14 indicates that the binding energies of gefitinib and AEE788 to the active state of the T790M mutant EGFR is 3 kcal/mol higher than to the wild type (WT), whereas ATP binding energy to the mutant is similar to the WT. Using metadynamics MD simulations with NAMD v2.9, the conformational equilibrium between the inactive resting state and the catalytically competent activate state was determined for the WT EGFR. When combined with the results obtained by Sutto and Gervasio, our simulations showed that the T790M point mutation lowers the free energy of the active state by 5 kcal/mol relative to the inactive state of the enzyme. Relative to the WT, the T790M mutant binds gefitinib more strongly. The T790M mutation is nevertheless resistant due to its increased binding of ATP. By contrast, the binding of AEE788 to the active state causes a conformational change in the αC-helix adjacent to the inhibitor binding pocket, that results in a 2 kcal/mol energy penalty. The energy penalty explains why the binding of AEE788 to the T790M mutant is unfavorable relative to binding to WT EGFR. These results establish the role of the gatekeeper mutation on inhibitor selectivity. Additional molecular dynamics (MD) simulations, TI, and metadynamics MD simulations reveal the origins of the changes in binding energy of WT and mutants.

The Effect of Molecular Conformation on the Accuracy of Theoretical (1)H and (13)C Chemical Shifts Calculated by Ab Initio Methods for Metabolic Mixture Analysis

NMR spectroscopy is a powerful method for analyzing metabolic mixtures. The information obtained from an NMR spectrum is in the form of physical parameters, such as chemical shifts, and construction of databases for many metabolites will be useful for data interpretation. To increase the accuracy of theoretical chemical shifts for development of a database for a variety of metabolites, the effects of sets of conformations (structural ensembles) and the levels of theory on computations of theoretical chemical shifts were systematically investigated for a set of 29 small molecules in the present study. For each of the 29 compounds, 101 structures were generated by classical molecular dynamics at 298.15 K, and then theoretical chemical shifts for 164 (1)H and 123 (13)C atoms were calculated by ab initio quantum chemical methods. Six levels of theory were used by pairing Hartree-Fock, B3LYP (density functional theory), or second order Møller-Plesset perturbation with 6-31G or aug-cc-pVDZ basis set. The six average fluctuations in the (1)H chemical shift were ±0.63, ± 0.59, ± 0.70, ± 0.62, ± 0.75, and ±0.66 ppm for the structural ensembles, and the six average errors were ±0.34, ± 0.27, ± 0.32, ± 0.25, ± 0.32, and ±0.25 ppm. The results showed that chemical shift fluctuations with changes in the conformation because of molecular motion were larger than the differences between computed and experimental chemical shifts for all six levels of theory. In conclusion, selection of an appropriate structural ensemble should be performed before theoretical chemical shift calculations for development of an accurate database for a variety of metabolites.

Long-Range Conformational Response of a PDZ Domain to Ligand Binding and Release: A Molecular Dynamics Study

The binding of a ligand to a protein may induce long-range structural or dynamical changes in the biomacromolecule even at sites physically well separated from the binding pocket. A system for which such behavior has been widely discussed is the PDZ2 domain of human tyrosine phosphatase 1E. Here, we present results from equilibrium trajectories of the PDZ2 domain in the free and ligand-bound state, as well as nonequilibrium simulations of the relaxation of PDZ2 after removal of its peptide ligand. The study reveals changes in inter-residue contacts, backbone dihedral angles, and C(α) positions upon ligand release. Our findings show a long-range conformational response of the PDZ2 domain to ligand release in the form of a collective shift of the secondary structure elements α2, β2, β3, α1-β4, and the C terminal loop relative to the rest of the protein away from the N-terminus, and a shift of the loops β2-β3 and β1-β2 in the opposite direction. The shifts lead to conformational changes in the backbone, especially in the β2-β3 loop but also in the β5-α2 and the α2-β6 loop, and are accompanied by changes of inter-residue contacts mainly within the β2-β3 loop as well as between the α2 helix and other segments. The residues showing substantial changes of inter-residue contacts, backbone conformations, or C(α) positions are considered "key residues" for the long-range conformational response of PDZ2. By comparing these residues with various sets of residues highlighted by previous studies of PDZ2, we investigate the statistical correlation of the various approaches. Interestingly, we find a considerable correlation of our findings with several works considering structural changes but no significant correlations with approaches considering energy flow or networks based on inter-residue energies.

Allosteric inhibition of c-Met kinase in sub-microsecond molecular dynamics simulations induced by its inhibitor, tivantinib

Molecular dynamics simulations showed that conformation transition of c-Met from DFG-in to DFG-out may accomplish rapidly in the presence of tivantinib. A unique binding mode of tivantinib was found to be critical for this “DFG-flip”.

Characterizing a partially ordered miniprotein through folding molecular dynamics simulations: Comparison with the experimental data

The villin headpiece helical subdomain (HP36) is one of the best known model systems for computational studies of fast-folding all-α miniproteins. HP21 is a peptide fragment-derived from HP36-comprising only the first and second helices of the full domain. Experimental studies showed that although HP21 is mostly unfolded in solution, it does maintain some persistent native-like structure as indicated by the analysis of NMR-derived chemical shifts. Here we compare the experimental data for HP21 with the results obtained from a 15-μs long folding molecular dynamics simulation performed in explicit water and with full electrostatics. We find that the simulation is in good agreement with the experiment and faithfully reproduces the major experimental findings, namely that (a) HP21 is disordered in solution with <10% of the trajectory corresponding to transiently stable structures, (b) the most highly populated conformer is a native-like structure with an RMSD from the corresponding portion of the HP36 crystal structure of <1 Å, (c) the simulation-derived chemical shifts-over the whole length of the trajectory-are in reasonable agreement with the experiment giving reduced χ(2) values of 1.6, 1.4, and 0.8 for the Δδ(13) C(α) , Δδ(13) CO, and Δδ(13) C(β) secondary shifts, respectively (becoming 0.8, 0.7, and 0.3 when only the major peptide conformer is considered), and finally, (d) the secondary structure propensity scores are in very good agreement with the experiment and clearly indicate the higher stability of the first helix. We conclude that folding molecular dynamics simulations can be a useful tool for the structural characterization of even marginally stable peptides.