Deconstructing the native state: energy landscapes, function, and dynamics of globular proteins

NumSimEX: A method using EXX hydrogen exchange mass spectrometry to map the energetics of protein folding landscapes

Hydrogen exchange mass spectrometry (HXMS) is a powerful tool to understand protein folding pathways and energetics. However, HXMS experiments to date have used exchange conditions termed EX1 or EX2 which limit the information that can be gained compared to the more general EXX exchange regime. If EXX behavior could be understood and analyzed, a single HXMS timecourse on an intact protein could fully map its folding landscape without requiring denaturation. To address this challenge, we developed a numerical simulation method called NumSimEX that models EXX exchange for arbitrarily complex folding pathways. NumSimEx fits protein folding dynamics to experimental HXMS data by iteratively comparing the simulated and experimental timecourses, allowing for determination of both kinetic and thermodynamic protein folding parameters. After analytically verifying NumSimEX's accuracy, we demonstrated its power on HXMS data from beta-2 microglobulin (β2M), a protein involved in dialysis-related amyloidosis. In particular, using NumSimEX, we identified three-state kinetics that near-perfectly matched experimental observation. This proof-of-principle application of NumSimEX sets the stage for harnessing HXMS to expand our understanding of proteins currently excluded from traditional protein folding methods. NumSimEX is freely available at https://github.com/JaswalLab/NumSimEX_Public.

Effect of Pressure on the Conformational Landscape of Human γD-Crystallin from Replica Exchange Molecular Dynamics Simulations

Human γD-crystallin belongs to a crucial family of proteins known as crystallins located in the fiber cells of the human lens. Since crystallins do not undergo any turnover after birth, they need to possess remarkable thermodynamic stability. However, their sporadic misfolding and aggregation, triggered by environmental perturbations or genetic mutations, constitute the molecular basis of cataracts, which is the primary cause of blindness in the globe according to the World Health Organization. Here, we investigate the impact of high pressure on the conformational landscape of wild-type HγD-crystallin using replica exchange molecular dynamics simulations augmented with principal component analysis. We find pressure to have a modest impact on global measures of protein stability, such as root-mean-square displacement and radius of gyration. Upon projecting our trajectories along the first two principal components from principal component analysis, however, we observe the emergence of distinct free energy basins at high pressures. By screening local order parameters previously shown or hypothesized as markers of HγD-crystallin stability, we establish correlations between a tyrosine-tyrosine aromatic contact within the N-terminal domain and the protein's end-to-end distance with projections along the first and second principal components, respectively. Furthermore, we observe the simultaneous contraction of the hydrophobic core and its intrusion by water molecules. This exploration sheds light on the intricate responses of HγD-crystallin to elevated pressures, offering insights into potential mechanisms underlying its stability and susceptibility to environmental perturbations, crucial for understanding cataract formation.

Permutationally Invariant Networks for Enhanced Sampling (PINES): Discovery of Multimolecular and Solvent-Inclusive Collective Variables

The typically rugged nature of molecular free-energy landscapes can frustrate efficient sampling of the thermodynamically relevant phase space due to the presence of high free-energy barriers. Enhanced sampling techniques can improve phase space exploration by accelerating sampling along particular collective variables (CVs). A number of techniques exist for the data-driven discovery of CVs parametrizing the important large-scale motions of the system. A challenge to CV discovery is learning CVs invariant to the symmetries of the molecular system, frequently rigid translation, rigid rotation, and permutational relabeling of identical particles. Of these, permutational invariance has proved a persistent challenge in frustrating the data-driven discovery of multimolecular CVs in systems of self-assembling particles and solvent-inclusive CVs for solvated systems. In this work, we integrate permutation invariant vector (PIV) featurizations with autoencoding neural networks to learn nonlinear CVs invariant to translation, rotation, and permutation and perform interleaved rounds of CV discovery and enhanced sampling to iteratively expand the sampling of configurational phase space and obtain converged CVs and free-energy landscapes. We demonstrate the permutationally invariant network for enhanced sampling (PINES) approach in applications to the self-assembly of a 13-atom argon cluster, association/dissociation of a NaCl ion pair in water, and hydrophobic collapse of a C45H92 n-pentatetracontane polymer chain. We make the approach freely available as a new module within the PLUMED2 enhanced sampling libraries.

Girsanov Reweighting Enhanced Sampling Technique (GREST): On-the-Fly Data-Driven Discovery of and Enhanced Sampling in Slow Collective Variables

Molecular dynamics simulations of microscopic phenomena are limited by the short integration time steps which are required for numerical stability but which limit the practically achievable simulation time scales. Collective variable (CV) enhanced sampling techniques apply biases to predefined collective coordinates to promote barrier crossing, phase space exploration, and sampling of rare events. The efficacy of these techniques is contingent on the selection of good CVs correlated with the molecular motions governing the long-time dynamical evolution of the system. In this work, we introduce Girsanov Reweighting Enhanced Sampling Technique (GREST) as an adaptive sampling scheme that interleaves rounds of data-driven slow CV discovery and enhanced sampling along these coordinates. Since slow CVs are inherently dynamical quantities, a key ingredient in our approach is the use of both thermodynamic and dynamical Girsanov reweighting corrections for rigorous estimation of slow CVs from biased simulation data. We demonstrate our approach on a toy 1D 4-well potential, a simple biomolecular system alanine dipeptide, and the Trp-Leu-Ala-Leu-Leu (WLALL) pentapeptide. In each case GREST learns appropriate slow CVs and drives sampling of all thermally accessible metastable states starting from zero prior knowledge of the system. We make GREST accessible to the community via a publicly available open source Python package.

Nonlinear Reaction Coordinate of an Enzyme Catalyzed Proton Transfer Reaction

We present an in-depth study on the theoretical calculation of an optimum reaction coordinate as a linear or nonlinear combination of important collective variables (CVs) sampled from an ensemble of reactive transition paths for an intramolecular proton transfer reaction catalyzed by the enzyme human carbonic anhydrase (HCA) II. The linear models are optimized by likelihood maximization for a given number of CVs. The nonlinear models are based on an artificial neural network with the same number of CVs and optimized by minimizing the root-mean-square error in comparison to a training set of committor estimators generated for the given transition. The nonlinear reaction coordinate thus obtained yields the free energy of activation and rate constant as 9.46 kcal mol^-1 and 1.25 × 10⁶ s^-1, respectively. These estimates are found to be in quantitative agreement with the known experimental results. We have also used an extended autoencoder model to show that a similar analysis can be carried out using a single CV only. The resultant free energies and kinetics of the reaction slightly overestimate the experimental data. The implications of these results are discussed using a detailed microkinetic scheme of the proton transfer reaction catalyzed by HCA II.

Bioinformatic and computational analysis for predominant mutations of the Nrf2/Keap1 complex in pediatric leukemia

The levels of reactive oxygen species (ROS) are tightly controlled and regulated by Nuclear Factor Erythroid-2-Like 2 (Nrf2) transcription factor, which is the main regulator of antioxidant responses and its suppressor protein Kelch-like ECH-associated protein 1 (Keap1). Our previous study has identified six novel changes in Nrf2/Keap1 pathway in pediatric ALL, which were described for the first time. These changes in the pathway are likely to alter the evolutionary process of amino acids and cause structural changes in the final products of genes. In this study, we aimed to compare the pathogenicity of eight determined mutations reported in our previous study by utilizing different programs with different algorithms and molecular dynamics simulation. Since it is too difficult to handle each existing mutation in a wet laboratory, in silico methods may give suggestion to choose the important mutations for further analysis and to establish the appropriate patient population and conduct wet laboratory studies. For this purpose, four different algorithms were used to evaluate the effects of single amino acid mutation. In addition, root-mean-square deviation, root-mean-square fluctuation and free-energy landscape analyses were performed to observe stability, flexibility and energetically favorable conformations, respectively, for each amino acid mutation. As a result, our study emphasizes the importance of Keap1 mutations in pediatric ALL Nrf2/Keap1 pathway, a total of eight mutations, two of which were shown for the first time in our study. Especially the mutations in the Keap1 Broad-Complex, Tramtrack and Bric-à-brac domain are worthy of attention.Communicated by Ramaswamy H. Sarma.

Analytical Method Using a Scaled Hypersphere Search for High-Dimensional Metadynamics Simulations

Metadynamics (MTD) is one of the most effective methods for calculating the free energy surface and finding rare events. Nevertheless, numerous studies using MTD have been carried out using 3D or lower dimensional collective variables (CVs), as higher dimensional CVs require costly computational resources and the obtained results are too complex to understand the important events. The latter issue can be conveniently solved by utilizing the free energy reaction network (FERN), which is a graph structure consisting of edges of minimum free energy paths (MFEPs), nodes of equation (EQ) points, and transition state (TS) points. In the present article, a new method for exploring FERNs on high-dimensional CVs using MTD and the scaled hypersphere search (SHS) method is described. A test calculation based on the MTD-SHS simulation of met-enkephalin in explicit water with 7 CVs was conducted. As a result, 889 EQ points and 1805 TS points were found. The MTD-SHS approach can find MFEPs exhaustively; therefore, the FERNs can be estimated without any a priori knowledge of the EQ and TS points.

Learning dynamical information from static protein and sequencing data

Many complex processes, from protein folding to neuronal network dynamics, can be described as stochastic exploration of a high-dimensional energy landscape. Although efficient algorithms for cluster detection in high-dimensional spaces have been developed over the last two decades, considerably less is known about the reliable inference of state transition dynamics in such settings. Here we introduce a flexible and robust numerical framework to infer Markovian transition networks directly from time-independent data sampled from stationary equilibrium distributions. We demonstrate the practical potential of the inference scheme by reconstructing the network dynamics for several protein-folding transitions, gene-regulatory network motifs, and HIV evolution pathways. The predicted network topologies and relative transition time scales agree well with direct estimates from time-dependent molecular dynamics data, stochastic simulations, and phylogenetic trees, respectively. Owing to its generic structure, the framework introduced here will be applicable to high-throughput RNA and protein-sequencing datasets, and future cryo-electron microscopy (cryo-EM) data.

Using Dimensionality Reduction to Analyze Protein Trajectories

In recent years the analysis of molecular dynamics trajectories using dimensionality reduction algorithms has become commonplace. These algorithms seek to find a low-dimensional representation of a trajectory that is, according to a well-defined criterion, optimal. A number of different strategies for generating projections of trajectories have been proposed but little has been done to systematically compare how these various approaches fare when it comes to analysing trajectories for biomolecules in explicit solvent. In the following paper, we have thus analyzed a molecular dynamics trajectory of the C-terminal fragment of the immunoglobulin binding domain B1 of protein G of Streptococcus modeled in explicit solvent using a range of different dimensionality reduction algorithms. We have then tried to systematically compare the projections generated using each of these algorithms by using a clustering algorithm to find the positions and extents of the basins in the high-dimensional energy landscape. We find that no algorithm outshines all the other in terms of the quality of the projection it generates. Instead, all the algorithms do a reasonable job when it comes to building a projection that separates some of the configurations that lie in different basins. Having said that, however, all the algorithms struggle to project the basins because they all have a large intrinsic dimensionality.

Simulating and analysing configurational landscapes of protein-protein contact formation

Interacting proteins can form aggregates and protein-protein interfaces with multiple patterns and different stabilities. Using molecular simulation one would like to understand the formation of these aggregates and which of the observed states are relevant for protein function and recognition. To characterize the complex configurational ensemble of protein aggregates, one needs a quantitative measure for the similarity of structures. We present well-suited descriptors that capture the essential features of non-covalent protein contact formation and domain motion. This set of collective variables is used with a nonlinear multi-dimensional scaling-based dimensionality reduction technique to obtain a low-dimensional representation of the configurational landscape of two ubiquitin proteins from coarse-grained simulations. We show that this two-dimensional representation is a powerful basis to identify meaningful states in the ensemble of aggregated structures and to calculate distributions and free energy landscapes for different sets of simulations. By using a measure to quantitatively compare free energy landscapes we can show how the introduction of a covalent bond between two ubiquitin proteins at different positions alters the configurational states of these dimers.

Interspecies Bombolitins Exhibit Structural Diversity upon Membrane Binding, Leading to Cell Specificity

Bombolitins, a class of peptides produced by bees of the genus Bombus, target and disrupt cellular membranes, leading to lysis. Antimicrobial peptides exhibit various mechanisms of action resulting from the interplay between peptide structure, lipid composition, and cellular target membrane selectivity. Herein, two bombolitins displaying significant amino-acid-sequence similarity, BII and BL6, were assessed for antimicrobial activity as well as correlated dodecylphosphocholine (DPC) micelle binding and membrane-induced peptide conformational changes. Infrared and circular dichroism spectroscopies were used to assess the structure-function relationship of each bombolitin, and the results indicate that BII forms a rigid and helically ordered secondary structure upon binding to DPC micelles, whereas BL6 largely lacks secondary structural order. Moreover, the binding affinity of each peptide to DPC micelles was determined, revealing that BL6 displayed a difference in binding affinity by over two orders of magnitude. Further investigations into the growth-inhibitory activity of the two bombolitins were performed against Escherichia coli and Saccharomyces cerevisiae. Interestingly, BII specifically targeted S. cerevisiae, whereas BL6 more effectively inhibited E. coli growth. Overall, the antimicrobial selectivity and specificity of BII and BL6 are largely dependent on the primary as well as secondary structural content of the peptides and the membrane composition.

Using Data-Reduction Techniques to Analyze Biomolecular Trajectories

This chapter discusses the way in which dimensionality reduction algorithms such as diffusion maps and sketch-map can be used to analyze molecular dynamics trajectories. The first part discusses how these various algorithms function as well as practical issues such as landmark selection and how these algorithms can be used when the data to be analyzed comes from enhanced sampling trajectories. In the later part a comparison between the results obtained by applying various algorithms to two sets of sample data is performed and discussed. This section is then followed by a summary of how one algorithm in particular, sketch-map, has been applied to a range of problems. The chapter concludes with a discussion on the directions that we believe this field is currently moving.

Recovery of Protein Folding Funnels from Single-Molecule Time Series by Delay Embeddings and Manifold Learning

The stability and folding of proteins is governed by the underlying single-molecule free energy surface (smFES) mapping the free energy of the molecule as a function of configurational state. Ascertaining the smFES is of great value in understanding and engineering protein structure and function. By integrating tools from dynamical systems theory and nonlinear manifold learning, we describe an approach to reconstruct the multidimensional smFES for a protein from a time series in a single experimentally measurable observable. We employ Takens' delay embeddings to project the time series into a high-dimensional space in which the projected dynamics are C¹-equivalent to the true system dynamics and employ diffusion maps to recover a low-dimensional reconstruction of the smFES that is equivalent to the true smFES up to a smooth and invertible transformation. We validate the approach in molecular dynamics simulations of Trp-cage, Villin, and BBA to demonstrate that landscapes recovered from univariate time series in the head-to-tail distance are topologically identical-they precisely preserve the metastable states and folding pathways-and topographically approximate-the free energy barrier heights and well depths are approximately preserved-to the true landscapes determined from complete knowledge of all atomic coordinates. We go on to show that the reconstructed landscapes reliably predict temperature denaturation and identify point mutations and groups of mutations critical to folding. These results demonstrate that protein folding funnels can be reconstructed from experimentally measurable time series and used to understand and engineer folding.

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

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

Building maps in collective variable space

Enhanced sampling techniques such as umbrella sampling and metadynamics are now routinely used to provide information on how the thermodynamic potential, or free energy, depends on a small number of collective variables (CVs). The free energy surfaces that one extracts by using these techniques provide a simplified or coarse-grained representation of the configurational ensemble. In this work, we discuss how auxiliary variables can be mapped in CV space. We show that maps of auxiliary variables allow one to analyze both the physics of the molecular system under investigation and the quality of the reduced representation of the system that is encoded in a set of CVs. We apply this approach to analyze the degeneracy of CVs and to compute entropy and enthalpy surfaces in CV space both for conformational transitions in alanine dipeptide and for phase transitions in carbon dioxide molecular crystals under pressure.

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.

Machine learning and data science in soft materials engineering

In many branches of materials science it is now routine to generate data sets of such large size and dimensionality that conventional methods of analysis fail. Paradigms and tools from data science and machine learning can provide scalable approaches to identify and extract trends and patterns within voluminous data sets, perform guided traversals of high-dimensional phase spaces, and furnish data-driven strategies for inverse materials design. This topical review provides an accessible introduction to machine learning tools in the context of soft and biological materials by 'de-jargonizing' data science terminology, presenting a taxonomy of machine learning techniques, and surveying the mathematical underpinnings and software implementations of popular tools, including principal component analysis, independent component analysis, diffusion maps, support vector machines, and relative entropy. We present illustrative examples of machine learning applications in soft matter, including inverse design of self-assembling materials, nonlinear learning of protein folding landscapes, high-throughput antimicrobial peptide design, and data-driven materials design engines. We close with an outlook on the challenges and opportunities for the field.

Histidine substitution in the most flexible fragments of firefly luciferase modifies its thermal stability

Molecular dynamics (MD) at two temperatures of 300 and 340 K identified two histidine residues, His461 and His489, in the most flexible regions of firefly luciferase, a light emitting enzyme. We therefore designed four protein mutants H461D, H489K, H489D and H489M to investigate their enzyme kinetic and thermodynamic stability changes. Substitution of His461 by aspartate (H461D) decreased ATP binding affinity, reduced the melting temperature of protein by around 25 °C and shifted its optimum temperature of activity to 10 °C. In line with the common feature of psychrophilic enzymes, the MD data showed that the overall flexibility of H461D was relatively high at low temperature, probably due to a decrease in the number of salt bridges around the mutation site. On the other hand, substitution of His489 by aspartate (H489D) introduced a new salt bridge between the C-terminal and N-terminal domains and increased protein rigidity but only slightly improved its thermal stability. Similar changes were observed for H489K and, to a lesser degree, H489M mutations. Based on our results we conclude that the MD simulation-based rational substitution of histidines by salt-bridge forming residues can modulate conformational dynamics in luciferase and shift its optimal temperature activity.

BayesWHAM: A Bayesian approach for free energy estimation, reweighting, and uncertainty quantification in the weighted histogram analysis method

The weighted histogram analysis method (WHAM) is a powerful approach to estimate molecular free energy surfaces (FES) from biased simulation data. Bayesian reformulations of WHAM are valuable in proving statistically optimal use of the data and providing a transparent means to incorporate regularizing priors and estimate statistical uncertainties. In this work, we develop a fully Bayesian treatment of WHAM to generate statistically optimal FES estimates in any number of biasing dimensions under arbitrary choices of the Bayes prior. Rigorous uncertainty estimates are generated by Metropolis-Hastings sampling from the Bayes posterior. We also report a means to project the FES and its uncertainties into arbitrary auxiliary order parameters beyond those in which biased sampling was conducted. We demonstrate the approaches in applications of alanine dipeptide and the unthreading of a synthetic mimic of the astexin-3 lasso peptide. Open-source MATLAB and Python implementations of our codes are available for free public download. © 2017 Wiley Periodicals, Inc.

Role of single-point mutations and deletions on transition temperatures in ideal proteinogenic heteropolymer chains in the gas phase

A coarse-grained statistical mechanics-based model for ideal heteropolymer proteinogenic chains of non-interacting residues is presented in terms of the size K of the chain and the set of helical propensities [Formula: see text] associated with each residue j along the chain. For this model, we provide an algorithm to compute the degeneracy tensor [Formula: see text] associated with energy level [Formula: see text] where [Formula: see text] is the number of residues with a native contact in a given conformation. From these results, we calculate the equilibrium partition function [Formula: see text] and characteristic temperature [Formula: see text] at which a transition from a low to a high entropy states is observed. The formalism is applied to analyze the effect on characteristic temperatures [Formula: see text] of single-point mutations and deletions of specific amino acids [Formula: see text] along the chain. Two probe systems are considered. First, we address the case of a random heteropolymer of size K and given helical propensities [Formula: see text] on a conformational phase space. Second, we focus our attention to a particular set of neuropentapeptides, [Met-5] and [Leu-5] enkephalins whose thermodynamic stability is a key feature on their coupling to [Formula: see text] and [Formula: see text] receptors and the triggering of biochemical responses.

Nonlinear reconstruction of single-molecule free-energy surfaces from univariate time series

The stable conformations and dynamical fluctuations of polymers and macromolecules are governed by the underlying single-molecule free energy surface. By integrating ideas from dynamical systems theory with nonlinear manifold learning, we have recovered single-molecule free energy surfaces from univariate time series in a single coarse-grained system observable. Using Takens' Delay Embedding Theorem, we expand the univariate time series into a high dimensional space in which the dynamics are equivalent to those of the molecular motions in real space. We then apply the diffusion map nonlinear manifold learning algorithm to extract a low-dimensional representation of the free energy surface that is diffeomorphic to that computed from a complete knowledge of all system degrees of freedom. We validate our approach in molecular dynamics simulations of a C(24)H(50) n-alkane chain to demonstrate that the two-dimensional free energy surface extracted from the atomistic simulation trajectory is - subject to spatial and temporal symmetries - geometrically and topologically equivalent to that recovered from a knowledge of only the head-to-tail distance of the chain. Our approach lays the foundations to extract empirical single-molecule free energy surfaces directly from experimental measurements.

Machine learning assembly landscapes from particle tracking data

Bottom-up self-assembly offers a powerful route for the fabrication of novel structural and functional materials. Rational engineering of self-assembling systems requires understanding of the accessible aggregation states and the structural assembly pathways. In this work, we apply nonlinear machine learning to experimental particle tracking data to infer low-dimensional assembly landscapes mapping the morphology, stability, and assembly pathways of accessible aggregates as a function of experimental conditions. To the best of our knowledge, this represents the first time that collective order parameters and assembly landscapes have been inferred directly from experimental data. We apply this technique to the nonequilibrium self-assembly of metallodielectric Janus colloids in an oscillating electric field, and quantify the impact of field strength, oscillation frequency, and salt concentration on the dominant assembly pathways and terminal aggregates. This combined computational and experimental framework furnishes new understanding of self-assembling systems, and quantitatively informs rational engineering of experimental conditions to drive assembly along desired aggregation pathways.

Thermodynamic and functional characteristics of deep-sea enzymes revealed by pressure effects

Hydrostatic pressure analysis is an ideal approach for studying protein dynamics and hydration. The development of full ocean depth submersibles and high pressure biological techniques allows us to investigate enzymes from deep-sea organisms at the molecular level. The aim of this review was to overview the thermodynamic and functional characteristics of deep-sea enzymes as revealed by pressure axis analysis after giving a brief introduction to the thermodynamic principles underlying the effects of pressure on the structural stability and function of enzymes.

Mapping polymerization and allostery of hemoglobin S using point mutations

Hemoglobin is a complex system that undergoes conformational changes in response to oxygen, allosteric effectors, mutations, and environmental changes. Here, we study allostery and polymerization of hemoglobin and its variants by application of two previously described methods: (i) AllosMod for simulating allostery dynamics given two allosterically related input structures and (ii) a machine-learning method for dynamics- and structure-based prediction of the mutation impact on allostery (Weinkam et al. J. Mol. Biol. 2013, 425, 647-661), now applicable to systems with multiple coupled binding sites, such as hemoglobin. First, we predict the relative stabilities of substates and microstates of hemoglobin, which are determined primarily by entropy within our model. Next, we predict the impact of 866 annotated mutations on hemoglobin's oxygen binding equilibrium. We then discuss a subset of 30 mutations that occur in the presence of the sickle cell mutation and whose effects on polymerization have been measured. Seven of these HbS mutations occur in three predicted druggable binding pockets that might be exploited to directly inhibit polymerization; one of these binding pockets is not apparent in the crystal structure, but only in structures generated by AllosMod. For the 30 mutations, we predict that mutation-induced conformational changes within a single tetramer tend not to significantly impact polymerization; instead, these mutations more likely impact polymerization by directly perturbing a polymerization interface. Finally, our analysis of allostery allows us to hypothesize why hemoglobin evolved to have multiple subunits and a persistent low frequency sickle cell mutation.

Population based reweighting of scaled molecular dynamics

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

Correlation between local structural dynamics of proteins inferred from NMR ensembles and evolutionary dynamics of homologues of known structure

Conformational changes in proteins are extremely important for their biochemical functions. Correlation between inherent conformational variations in a protein and conformational differences in its homologues of known structure is still unclear. In this study, we have used a structural alphabet called Protein Blocks (PBs). PBs are used to perform abstraction of protein 3-D structures into a 1-D strings of 16 alphabets (a-p) based on dihedral angles of overlapping pentapeptides. We have analyzed the variations in local conformations in terms of PBs represented in the ensembles of 801 protein structures determined using NMR spectroscopy. In the analysis of concatenated data over all the residues in all the NMR ensembles, we observe that the overall nature of inherent local structural variations in NMR ensembles is similar to the nature of local structural differences in homologous proteins with a high correlation coefficient of .94. High correlation at the alignment positions corresponding to helical and β-sheet regions is only expected. However, the correlation coefficient by considering only the loop regions is also quite high (.91). Surprisingly, segregated position-wise analysis shows that this high correlation does not hold true to loop regions at the structurally equivalent positions in NMR ensembles and their homologues of known structure. This suggests that the general nature of local structural changes is unique; however most of the local structural variations in loop regions of NMR ensembles do not correlate to their local structural differences at structurally equivalent positions in homologues.

Thermodynamics of ideal proteinogenic homopolymer chains as a function of the energy spectrum E, helical propensity ω and enthalpic energy barrier

A reformulation and generalization of the Zwanzig model (ZW model) for ideal homopolymer chains poly-X, where X represents any of the twenty naturally occurring proteinogenic amino acid residues is presented. This reformulation and generalization provides a direct connection between coarse-grained parameters originally proposed in the ZW model with variables from the Lifson-Roig (LR) theory, such as the helical propensity per residue ω, and new variables introduced here, such as the energy gap Δ between unfolded and folded structures, as well as the ratio f of the energy scales involved. This enables us to discover the relevance of the energy spectrum E to the onset of configurational phase transitions. From the configurational partition function Q, thermodynamic properties such as the configurational entropy S, specific heat v and average energy <E> are calculated in terms of the number of residues K, temperature T, helical propensity ω and energy barrier ΔH for different poly-X chains in vacuo. Results obtained here provide substantial evidence that configurational phase transitions for ideal poly-X chains correspond to first-order phase transitions. An anomalous behavior of the thermodynamic functions <E>, Cv, S with respect to the number K of residues is also highlighted. On-going methods of solution are outlined.

Intramolecular interactions stabilizing compact conformations of the intrinsically disordered kinase-inhibitor domain of Sic1: a molecular dynamics investigation

Cyclin-dependent kinase inhibitors (CKIs) are key regulatory proteins of the eukaryotic cell cycle, which modulate cyclin-dependent kinase (Cdk) activity. CKIs perform their inhibitory effect by the formation of ternary complexes with a target kinase and its cognate cyclin. These regulators generally belong to the class of intrinsically disordered proteins (IDPs), which lack a well-defined and organized three-dimensional (3D) structure in their free state, undergoing folding upon binding to specific partners. Unbound IDPs are not merely random-coil structures, but can present intrinsically folded structural units (IFSUs) and collapsed conformations. These structural features can be relevant to protein function in vivo. The yeast CKI Sic1 is a 284-amino acid IDP that binds to Cdk1 in complex with the Clb5,6 cyclins, preventing phosphorylation of G1 substrates and, therefore, entrance to the S phase. Sic1 degradation, triggered by multiple phosphorylation events, promotes cell-cycle progression. Previous experimental studies pointed out a propensity of Sic1 and its isolated domains to populate both extended and compact conformations. The present contribution provides models for compact conformations of the Sic1 kinase-inhibitory domain (KID) by all-atom molecular dynamics (MD) simulations in explicit solvent and in the absence of interactors. The results are integrated by spectroscopic and spectrometric data. Helical IFSUs are identified, along with networks of intramolecular interactions. The results identify a group of putative hub residues and networks of electrostatic interactions, which are likely to be involved in the stabilization of the globular states.

Can the energy gap in the protein-ligand binding energy landscape be used as a descriptor in virtual ligand screening?

The ranking of scores of individual chemicals within a large screening library is a crucial step in virtual screening (VS) for drug discovery. Previous studies showed that the quality of protein-ligand recognition can be improved using spectrum properties and the shape of the binding energy landscape. Here, we investigate whether the energy gap, defined as the difference between the lowest energy pose generated by a docking experiment and the average energy of all other generated poses and inferred to be a measure of the binding energy landscape sharpness, can improve the separation power between true binders and decoys with respect to the use of the best docking score. We performed retrospective single- and multiple-receptor conformation VS experiments in a diverse benchmark of 40 domains from 38 therapeutically relevant protein targets. Also, we tested the performance of the energy gap on 36 protein targets from the Directory of Useful Decoys (DUD). The results indicate that the energy gap outperforms the best docking score in its ability to discriminate between true binders and decoys, and true binders tend to have larger energy gaps than decoys. Furthermore, we used the energy gap as a descriptor to measure the height of the native binding phase and obtained a significant increase in the success rate of near native binding pose identification when the ligand binding conformations within the boundaries of the native binding phase were considered. The performance of the energy gap was also evaluated on an independent test case of VS-identified PKR-like ER-localized eIF2α kinase (PERK) inhibitors. We found that the energy gap was superior to the best docking score in its ability to more highly rank active compounds from inactive ones. These results suggest that the energy gap of the protein-ligand binding energy landscape is a valuable descriptor for use in VS.

Regulation of the H4 tail binding and folding landscapes via Lys-16 acetylation

Intrinsically disordered proteins (IDP) are a broad class of proteins with relatively flat energy landscapes showing a high level of functional promiscuity, which are frequently regulated through posttranslational covalent modifications. Histone tails, which are the terminal segments of the histone proteins, are prominent IDPs that are implicated in a variety of signaling processes, which control chromatin organization and dynamics. Although a large body of work has been done on elucidating the roles of posttranslational modifications in functional regulation of IDPs, molecular mechanisms behind the observed behaviors are not fully understood. Using extensive atomistic molecular dynamics simulations, we found in this work that H4 tail mono-acetylation at LYS-16, which is a key covalent modification, induces a significant reorganization of the tail's conformational landscape, inducing partial ordering and enhancing the propensity for alpha-helical segments. Furthermore, our calculations of the potentials of mean force between the H4 tail and a DNA fragment indicate that contrary to the expectations based on simple electrostatic reasoning, the Lys-16 mono-acetylated H4 tail binds to DNA stronger than the unacetylated protein. Based on these results, we propose a molecular mechanism for the way Lys-16 acetylation might lead to experimentally observed disruption of compact chromatin fibers.

Principal component and clustering analysis on molecular dynamics data of the ribosomal L11·23S subdomain

With improvements in computer speed and algorithm efficiency, MD simulations are sampling larger amounts of molecular and biomolecular conformations. Being able to qualitatively and quantitatively sift these conformations into meaningful groups is a difficult and important task, especially when considering the structure-activity paradigm. Here we present a study that combines two popular techniques, principal component (PC) analysis and clustering, for revealing major conformational changes that occur in molecular dynamics (MD) simulations. Specifically, we explored how clustering different PC subspaces effects the resulting clusters versus clustering the complete trajectory data. As a case example, we used the trajectory data from an explicitly solvated simulation of a bacteria’s L11·23S ribosomal subdomain, which is a target of thiopeptide antibiotics. Clustering was performed, using K-means and average-linkage algorithms, on data involving the first two to the first five PC subspace dimensions. For the average-linkage algorithm we found that data-point membership, cluster shape, and cluster size depended on the selected PC subspace data. In contrast, K-means provided very consistent results regardless of the selected subspace. Since we present results on a single model system, generalization concerning the clustering of different PC subspaces of other molecular systems is currently premature. However, our hope is that this study illustrates a) the complexities in selecting the appropriate clustering algorithm, b) the complexities in interpreting and validating their results, and c) by combining PC analysis with subsequent clustering valuable dynamic and conformational information can be obtained.

Coupled motions during dynamics reveal a tunnel toward the active site regulated by the N-terminal α-helix in an acylaminoacyl peptidase

Acylaminoacyl peptidase (AAP) subfamily belongs to the prolyl oligopeptidase (POP) family of serine-proteases. There is a great interest in the definition of molecular mechanisms related to the activity and substrate recognition of these complex multi-domain enzymes. The active site relies at the interface between the C-terminal catalytic domain and the β-propeller domain, whose N-terminal region acts as a bridge to the hydrolase domain. In AAP, the N-terminal extension is characterized by a structurally conserved α1-helix, which is known to affect thermal stability and thermal dependence of the catalytic activity. In the present contribution, results from hundreds nanosecond all-atom molecular dynamics simulations, along with analyses of the networks of cross-correlated motions of a member of the AAP subfamily are discussed. The MD investigation identifies a tunnel that from the surrounding of the N-terminal α1-helix bring to the catalytic site. This cavity seems to be regulated by conformational changes of the α1-helix itself during the dynamics. The evidence here provided can be a useful guide for a better understanding of the mechanistic aspects related to AAP activity, but also for drug design purposes.

Protein fluxes along the filopodium as a framework for understanding the growth-retraction dynamics: the interplay between diffusion and active transport

We present a picture of filopodial growth and retraction from physics perspective, where we emphasize the significance of the role played by protein fluxes due to spatially extended nature of the filopodium. We review a series of works, which used stochastic simulations and mean field analytical modeling to find the concentration profile of G-actin inside a filopodium, which, in turn, determines the stationary filopodial length. In addition to extensively reviewing the prior works, we also report some new results on the role of active transport in regulating the length of filopodia. We model a filopodium where delivery of actin monomers towards the tip can occur both through passive diffusion and active transport by myosin motors. We found that the concentration profile of G-actin along the filopodium is rather non-trivial, containing a narrow minimum near the base followed by a broad maximum. For efficient enough actin transport, this non-monotonous shape is expected to occur under a broad set of conditions. We also raise the issue of slow approach to the stationary length and the possibility of multiple steady state solutions.

Integrating diffusion maps with umbrella sampling: application to alanine dipeptide

Nonlinear dimensionality reduction techniques can be applied to molecular simulation trajectories to systematically extract a small number of variables with which to parametrize the important dynamical motions of the system. For molecular systems exhibiting free energy barriers exceeding a few k(B)T, inadequate sampling of the barrier regions between stable or metastable basins can lead to a poor global characterization of the free energy landscape. We present an adaptation of a nonlinear dimensionality reduction technique known as the diffusion map that extends its applicability to biased umbrella sampling simulation trajectories in which restraining potentials are employed to drive the system into high free energy regions and improve sampling of phase space. We then propose a bootstrapped approach to iteratively discover good low-dimensional parametrizations by interleaving successive rounds of umbrella sampling and diffusion mapping, and we illustrate the technique through a study of alanine dipeptide in explicit solvent.

Protonation/deprotonation effects on the stability of the Trp-cage miniprotein

The Trp-cage miniprotein is a 20 amino acid peptide that exhibits many of the properties of globular proteins. In this protein, the hydrophobic core is formed by a buried Trp side chain. The folded state is stabilized by an ion pair between aspartic acid and an arginine side chain. The effect of protonating the aspartic acid on the Trp-cage miniprotein folding/unfolding equilibrium is studied by explicit solvent molecular dynamics simulations of the protein in the charged and protonated Asp9 states. Unbiased Replica Exchange Molecular Dynamics (REMD) simulations, spanning a wide temperature range, are carried out to the microsecond time scale, using the AMBER99SB forcefield in explicit TIP3P water. The protein structural ensembles are studied in terms of various order parameters that differentiate the folded and unfolded states. We observe that in the folded state the root mean square distance (rmsd) from the backbone of the NMR structure shows two highly populated basins close to the native state with peaks at 0.06 nm and 0.16 nm, which are consistent with previous simulations using the same forcefield. The fraction of folded replicas shows a drastic decrease because of the absence of the salt bridge. However, significant populations of conformations with the arginine side chain exposed to the solvent, but within the folded basin, are found. This shows the possibility to reach the folded state without formation of the ion pair. We also characterize changes in the unfolded state. The equilibrium populations of the folded and unfolded states are used to characterize the thermodynamics of the system. We find that the change in free energy difference due to the protonation of the Asp amino acid is 3 kJ mol(-1) at 297 K, favoring the charged state, and resulting in ΔpK(1) = 0.5 units for Asp9. We also study the differences in the unfolded state ensembles for the two charge states and find significant changes at low temperature, where the protonated Asp side chain makes multiple hydrogen bonds to the protein backbone.

From the Cover: Simplifying the representation of complex free-energy landscapes using sketch-map

A new scheme, sketch-map, for obtaining a low-dimensional representation of the region of phase space explored during an enhanced dynamics simulation is proposed. We show evidence, from an examination of the distribution of pairwise distances between frames, that some features of the free-energy surface are inherently high-dimensional. This makes dimensionality reduction problematic because the data does not satisfy the assumptions made in conventional manifold learning algorithms We therefore propose that when dimensionality reduction is performed on trajectory data one should think of the resultant embedding as a quickly sketched set of directions rather than a road map. In other words, the embedding tells one about the connectivity between states but does not provide the vectors that correspond to the slow degrees of freedom. This realization informs the development of sketch-map, which endeavors to reproduce the proximity information from the high-dimensionality description in a space of lower dimensionality even when a faithful embedding is not possible.