Bridging between NMA and Elastic Network Models: Preserving All-Atom Accuracy in Coarse-Grained Models

Coarse-Graining Waters: Unveiling The Effective Hydrophilicity/Hydrophobicity of Individual Protein Atoms and The Roles of Waters' Hydrogens

There have been many coarse-graining methods developed that aim to reduce the sizes of simulated systems and their computational costs. In this work, we develop a new coarse-graining method, called coarse-graining-delta (or δ-CG in short), that reduces the degrees of freedom of the potential energy surface by coarse-graining relative locations of atoms from their unit centers. Our method extends and generalizes the methods used in the coarse-grained normal mode analysis and enables us to study the roles of the individual removed atoms in a system, which have been difficult to study in molecular dynamics simulations. By applying δ-CG to coarse-grain three-point water molecules into single-point solvent particles, we successfully identify the effective hydrophilicity and hydrophobicity of all the individual protein atom types, which collectively correlate well with the known hydrophilic, hydrophobic, and amphipathic characteristics of amino acids. Moreover, our investigation shows that water's hydrogens have two roles in interacting with protein atoms. First, water molecules adjust their poses around different amino acids and their atoms, and the statistical preferences of the hydrogen poses near the atoms determine the effective hydrophilicity and hydrophobicity of the atoms, which have not been successfully addressed before. Second, the collective dynamics of the hydrogens assist the water molecules in escaping from the potential energy wells of the hydrophilic atoms. Our method also shows that coarse-graining a system mathematically leads to breaking antisymmetry of the nonbonded interactions; as a result, two interacting coarse-grained units exert different forces on each other. Our study indicates that the accuracy of coarse-grained force fields, such as the MARTINI force field and the UNRES force field, can be improved in two ways: (i) refining their potential energy functions and coefficients by analyzing the coarse-grained potential energy surface using δ-CG, and (ii) introducing non-antisymmetric interactions.

Normal mode analysis and comparative study of intrinsic dynamics of alcohol oxidase enzymes from GMC protein family

Glucose-Methanol-Choline (GMC) family enzymes are very important in catalyzing the oxidation of a wide range of structurally diverse substrates. Enzymes that constitute the GMC family, share a common tertiary fold but < 25% sequence identity. Cofactor FAD, FAD binding signature motif, and similar structural scaffold of the active site are common features of oxidoreductase enzymes of the GMC family. Protein functionality mainly depends on protein three-dimensional structures and dynamics. In this study, we used the normal mode analysis method to search the intrinsic dynamics of GMC family enzymes. We have explored the dynamical behavior of enzymes with unique substrate catabolism and active site characteristics from different classes of the GMC family. Analysis of individual enzymes and comparative ensemble analysis of enzymes from different classes has shown conserved dynamic motion at FAD binding sites. The present study revealed that GMC enzymes share a strong dynamic similarity (Bhattacharyya coefficient >90% and root mean squared inner product >52%) despite low sequence identity across the GMC family enzymes. The study predicts that local deformation energy between atoms of the enzyme may be responsible for the catalysis of different substrates. This study may help that intrinsic dynamics can be used to make meaningful classifications of proteins or enzymes from different organisms.Communicated by Ramaswamy H. Sarma.

Coarse-Grained Models for Vault Normal Model Analysis

Recent experiments have shown that the molecular complex of vault has large conformational changes at its shoulder and cap regions in solution. From the comparison of two configuration structures, it has been found that the shoulder region can twist and move outward, while the cap region will rotate and push upward correspondingly. To further understand these experimental results, in this paper, we study the vault dynamics for the first time. Since vault has an extremely large-sized structure with around 63,336 Cα atoms, traditional normal mode method with the Cα coarse-grained representation will fall short. We employ a newly invented multiscale virtual particle-based anisotropic network model (MVP-ANM). To reduce the complexity, the 39-folder vault structure is coarse-grained to about 6000 virtual particles, which significantly reduces the computational cost while still maintaining the basic structure information. Among the 14 low frequency eigenmodes from Mode 7 to Mode 20, two eigenmodes, i.e., Mode 9 and Mode 20, are found to be directly associated with the experimental observations. In Mode 9, shoulder region undergoes a significant expansion while the cap part is lifted upward. In Mode 20, a clear rotation of both shoulder and cap regions is well observed. Our results are consistent with the experimental observations. More importantly, these low frequency eigenmodes indicate that the vault waist, shoulder and lower cap regions are the most likely regions for the opening of the vault particle. And the opening mechanism is highly likely to be rotation and expansion at these regions. As far as we know, this is the first work to provide the normal mode analysis for the vault complex.

Optimized Elastic Network Models With Direct Characterization of Inter-Residue Cooperativity for Protein Dynamics

The elastic network models (ENMs)are known as representative coarse-grained models to capture essential dynamics of proteins. Due to simple designs of the force constants as a decay with spatial distances of residue pairs in many previous studies, there is still much room for the improvement of ENMs. In this article, we directly computed the force constants with the inverse covariance estimation using a ridge-type operater for the precision matrix estimation (ROPE)on a large-scale set of NMR ensembles. Distance-dependent statistical analyses on the force constants were further comprehensively performed in terms of several paired types of sequence and structural information, including secondary structure, relative solvent accessibility, sequence distance and terminal. Various distinguished distributions of the mean force constants highlight the structural and sequential characteristics coupled with the inter-residue cooperativity beyond the spatial distances. We finally integrated these structural and sequential characteristics to build novel ENM variations using the particle swarm optimization for the parameter estimation. The considerable improvements on the correlation coefficient of the mean-square fluctuation and the mode overlap were achieved by the proposed variations when compared with traditional ENMs. This study opens a novel way to develop more accurate elastic network models for protein dynamics.

A Comparative Evaluation of the Structural and Dynamic Properties of Insect Odorant Binding Proteins

Insects devote a major part of their metabolic resources to the production of odorant binding proteins (OBPs). Although initially, these proteins were implicated in the solubilisation, binding and transport of semiochemicals to olfactory receptors, it is now recognised that they may play diverse, as yet uncharacterised, roles in insect physiology. The structures of these OBPs, the majority of which are known as “classical” OBPs, have shed some light on their potential functional roles. However, the dynamic properties of these proteins have received little attention despite their functional importance. Structural dynamics are encoded in the native protein fold and enable the adaptation of proteins to substrate binding. This paper provides a comparative review of the structural and dynamic properties of OBPs, making use of sequence/structure analysis, statistical and theoretical physics-based methods. It provides a new layer of information and additional methodological tools useful in unravelling the relationship between structure, dynamics and function of insect OBPs. The dynamic properties of OBPs, studied by means of elastic network models, reflect the similarities/dissimilarities observed in their respective structures and provides insights regarding protein motions that may have important implications for ligand recognition and binding. Furthermore, it was shown that the OBPs studied in this paper share conserved structural ‘core’ that may be of evolutionary and functional importance.

Modal Analysis of the Lysozyme Protein Considering All-Atom and Coarse-Grained Finite Element Models

Proteins are the fundamental entities of several organic activities. They are essential for a broad range of tasks in a way that their shapes and folding processes are crucial to achieving proper biological functions. Low-frequency modes, generally associated with collective movements at terahertz (THz) and sub-terahertz frequencies, have been appointed as critical for the conformational processes of many proteins. Dynamic simulations, such as molecular dynamics, are vastly applied by biochemical researchers in this field. However, in the last years, proposals that define the protein as a simplified elastic macrostructure have shown appealing results when dealing with this type of problem. In this context, modal analysis based on different modelization techniques, i.e., considering both an all-atom (AA) and coarse-grained (CG) representation, is proposed to analyze the hen egg-white lysozyme. This work presents new considerations and conclusions compared to previous analyses. Experimental values for the B-factor, considering all the heavy atoms or only one representative point per amino acid, are used to evaluate the validity of the numerical solutions. In general terms, this comparison allows the assessment of the regional flexibility of the protein. Besides, the low computational requirements make this approach a quick method to extract the protein’s dynamic properties under scrutiny.

Identification of key sites controlling protein functional motions by using elastic network model combined with internal coordinates

The elastic network model (ENM) is an effective method to extract the intrinsic dynamical properties encoded in protein tertiary structures. We have proposed a new ENM-based analysis method to reveal the motion modes directly responsible for a specific protein function, in which an internal coordinate related to the specific function was introduced to construct the internal/Cartesian hybrid coordinate space. In the present work, the function-related internal coordinates combined with a linear perturbation method were applied to identify the key sites controlling specific protein functional motions. The change in the fluctuations of the internal coordinate in response to residue perturbation was calculated in the hybrid coordinate space by using the linear response theory. The residues with the large fluctuation changes were identified to be the key sites that allosterically control the specific protein function. Two proteins, i.e., human DNA polymerase β and the chaperonin from Methanococcus maripaludis, were investigated as case studies, in which several collective and local internal coordinates were applied to identify the functionally key residues of these two studied proteins. The calculation results are consistent with the experimental observations. It is found that different collective internal coordinates lead to similar results, where the predicted functionally key sites are located at similar positions in the protein structure. While for the local internal coordinates, the predicted key sites tend to be situated at the region near to the coordinate-involving residues. Our studies provide a starting point for further exploring other function-related internal coordinates for other interesting proteins.

A complex multiscale virtual particle model based elastic network model (CMVP-ENM) for the normal mode analysis of biomolecular complexes

The CMVP-ENM for virus normal mode analysis. With a special ratio parameter, CMVP-ENM can characterize the multi-material properties of biomolecular complexes and systematically enhance or suppress the modes for different components.

A time and memory efficient recipe for fast normal mode computations of complexes with icosahedral symmetry

With the recent breakthroughs in experimental technologies, structure determination of extremely large assemblies, many with icosahedral symmetry, has been rapidly accelerating. Computational studies of their dynamics are important to deciphering their functions as well as to structural refinement but are challenged by their extremely large size, which ranges from hundreds of thousands to even millions of atoms. Group theory can be used to significantly speed up the normal mode computations of these symmetric complexes, but the derivation is often obscured by the complexity of group theory and consequently is not widely accessible. To address this problem, this work presents an easy recipe for normal mode computations of complexes with icosahedral symmetry. The recipe details how the Hessian matrix in symmetry coordinates can be constructed in a few easy steps of matrix multiplications, without going through the complexity of group theory. All the "ingredient" matrices required in the recipe are fully provided in the Supplemental Information for easy reproduction. The work is timely considering the expected large in-flux of many more icosahedral assemblies in the near future. The recipe uses a minimum amount of memory and solves the normal modes in a significantly reduced amount of time, making it feasible to perform normal mode computations of these assemblies on most computer systems.

Kinematic Flexibility Analysis: Hydrogen Bonding Patterns Impart a Spatial Hierarchy of Protein Motion

Elastic network models (ENMs) and constraint-based, topological rigidity analysis are two distinct, coarse-grained approaches to study conformational flexibility of macromolecules. In the two decades since their introduction, both have contributed significantly to insights into protein molecular mechanisms and function. However, despite a shared purpose of these approaches, the topological nature of rigidity analysis, and thereby the absence of motion modes, has impeded a direct comparison. Here, we present an alternative, kinematic approach to rigidity analysis, which circumvents these drawbacks. We introduce a novel protein hydrogen bond network spectral decomposition, which provides an orthonormal basis for collective motions modulated by noncovalent interactions, analogous to the eigenspectrum of normal modes. The zero modes decompose proteins into rigid clusters identical to those from topological rigidity, while nonzero modes rank protein motions by their hydrogen bond collective energy penalty. Our kinematic flexibility analysis bridges topological rigidity theory and ENM, enabling a detailed analysis of motion modes obtained from both approaches. Analysis of a large, structurally diverse data set revealed that collectivity of protein motions, reported by the Shannon entropy, is significantly reduced for rigidity theory compared to normal mode approaches. Strikingly, kinematic flexibility analysis suggests that the hydrogen bonding network encodes a protein-fold specific, spatial hierarchy of motions, which goes nearly undetected in ENM. This hierarchy reveals distinct motion regimes that rationalize experimental and simulated protein stiffness variations. Kinematic motion modes highly correlate with reported crystallographic B factors and molecular dynamics simulations of adenylate kinase. A formal expression for changes in free energy derived from the spectral decomposition indicates that motions across nearly 40% of modes obey enthalpy-entropy compensation. Taken together, our results suggest that hydrogen bond networks have evolved to modulate protein structure and dynamics, which can be efficiently probed by kinematic flexibility analysis.

All-atom normal mode dynamics of HIV-1 capsid

Dynamics of biomolecular assemblies offer invaluable insights into their functional mechanisms. For extremely large biomolecular systems, such as HIV-1 capsid that has nearly 5 millions atoms, obtaining its normal mode dynamics using even coarse-grained models can be a challenging task. In this work, we have successfully carried out a normal mode analysis of an entire HIV-1 capsid in full all-atom details. This is made possible through our newly developed BOSE (Block of Selected Elasticity) model that is founded on the principle of resonance discovered in our recent work. The resonance principle makes it possible to most efficiently compute the vibrations of a whole capsid at any given frequency by projecting the motions of component capsomeres into a narrow subspace. We have conducted also assessments of the quality of the BOSE modes by comparing them with benchmark modes obtained directly from the original Hessian matrix. Our all-atom normal mode dynamics study of the HIV-1 capsid reveals the dynamic role of the pentamers in stabilizing the capsid structure and is in agreement with experimental findings that suggest capsid disassembly and uncoating start when the pentamers become destabilized. Our results on the dynamics of hexamer pores suggest that nucleotide transport should take place mostly at hexamers near pentamers, especially at the larger hemispherical end.

Supramolecular assemblies are large biomolecular complexes composed of hundreds or even thousands of protein chains. They function as molecular machines or as large containers that store or facilitate the chemical reactions of other molecules. Whatever they do, their functional mechanisms are tightly linked to their structures and intrinsic dynamics. Recently, due to breakthroughs in experimental techniques, many supramolecular assemblies have been determined, such as the capsid of human immunodeficiency virus (HIV) that is composed of nearly 5 millions of atoms. Computational studies of these systems are challenging due to their extremely large sizes. In this work, we have successfully carried out a dynamics study of an entire HIV capsid in full all-atom details. Our study reveals new insights into the dynamics of the N-terminal loops, the stabilizing role of the pentamers, and where the nucleotide transport may take place.

Singular value decomposition for the correlation of atomic fluctuations with arbitrary angle

Many proteins exhibit a critical property called allostery, which enables intra-molecular transmission of information between distal sites. Microscopically, allosteric response is closely related to correlated atomic fluctuations. Conventional correlation analysis correlates the atomic fluctuations at two sites by taking the dot product (DP) between the fluctuations, which accounts only for the parallel and antiparallel components. Here, we present a singular value decomposition (SVD) method that analyzes the correlation coefficient of fluctuation dynamics with an arbitrary angle between the correlated directions. In a model allosteric system, the second PDZ domain (PDZ2) in the human PTP1E protein, approximately one third of the strong correlations have near-perpendicular directions, which are underestimated in the conventional method. The discrimination becomes more prominent for residue pairs with larger separation. The results of the proposed SVD method are more consistent with the experimentally determined PDZ2 dynamics than those of conventional method. In addition, the SVD method improved the prediction accuracy of the allosteric sites in a dataset of 23 known allosteric monomer proteins. The proposed method may inspire extended investigation not only into allostery, but also into protein dynamics and drug design.

Fast Normal Mode Computations of Capsid Dynamics Inspired by Resonance

Increasingly more and larger structural complexes are being determined experimentally. The sizes of these systems pose a formidable computational challenge to the study of their vibrational dynamics by normal mode analysis. To overcome this challenge, this work presents a novel resonance-inspired approach. Tests on large shell structures of protein capsids demonstrate there is a strong resonance between the vibrations of a whole capsid and those of individual capsomeres. We then show how this resonance can be taken advantage of to significantly speed up normal mode computations.

Multiscale virtual particle based elastic network model (MVP-ENM) for normal mode analysis of large-sized biomolecules

In this paper, a multiscale virtual particle based elastic network model (MVP-ENM) is proposed for the normal mode analysis of large-sized biomolecules. The multiscale virtual particle (MVP) model is proposed for the discretization of biomolecular density data. With this model, large-sized biomolecular structures can be coarse-grained into virtual particles such that a balance between model accuracy and computational cost can be achieved. An elastic network is constructed by assuming "connections" between virtual particles. The connection is described by a special harmonic potential function, which considers the influence from both the mass distributions and distance relations of the virtual particles. Two independent models, i.e., the multiscale virtual particle based Gaussian network model (MVP-GNM) and the multiscale virtual particle based anisotropic network model (MVP-ANM), are proposed. It has been found that in the Debye-Waller factor (B-factor) prediction, the results from our MVP-GNM with a high resolution are as good as the ones from GNM. Even with low resolutions, our MVP-GNM can still capture the global behavior of the B-factor very well with mismatches predominantly from the regions with large B-factor values. Further, it has been demonstrated that the low-frequency eigenmodes from our MVP-ANM are highly consistent with the ones from ANM even with very low resolutions and a coarse grid. Finally, the great advantage of MVP-ANM model for large-sized biomolecules has been demonstrated by using two poliovirus virus structures. The paper ends with a conclusion.

Combining protein sequence, structure, and dynamics: A novel approach for functional evolution analysis of PAS domain superfamily

PAS domains are widespread in archaea, bacteria, and eukaryota, and play important roles in various functions. In this study, we aim to explore functional evolutionary relationship among proteins in the PAS domain superfamily in view of the sequence-structure-dynamics-function relationship. We collected protein sequences and crystal structure data from RCSB Protein Data Bank of the PAS domain superfamily belonging to three biological functions (nucleotide binding, photoreceptor activity, and transferase activity). Protein sequences were aligned and then used to select sequence-conserved residues and build phylogenetic tree. Three-dimensional structure alignment was also applied to obtain structure-conserved residues. The protein dynamics were analyzed using elastic network model (ENM) and validated by molecular dynamics (MD) simulation. The result showed that the proteins with same function could be grouped by sequence similarity, and proteins in different functional groups displayed statistically significant difference in their vibrational patterns. Interestingly, in all three functional groups, conserved amino acid residues identified by sequence and structure conservation analysis generally have a lower fluctuation than other residues. In addition, the fluctuation of conserved residues in each biological function group was strongly correlated with the corresponding biological function. This research suggested a direct connection in which the protein sequences were related to various functions through structural dynamics. This is a new attempt to delineate functional evolution of proteins using the integrated information of sequence, structure, and dynamics.

Gaussian network model can be enhanced by combining solvent accessibility in proteins

Gaussian network model (GNM), regarded as the simplest and most representative coarse-grained model, has been widely adopted to analyze and reveal protein dynamics and functions. Designing a variation of the classical GNM, by defining a new Kirchhoff matrix, is the way to improve the residue flexibility modeling. We combined information arising from local relative solvent accessibility (RSA) between two residues into the Kirchhoff matrix of the parameter-free GNM. The undetermined parameters in the new Kirchhoff matrix were estimated by using particle swarm optimization. The usage of RSA was motivated by the fact that our previous work using RSA based linear regression model resulted out higher prediction quality of the residue flexibility when compared with the classical GNM and the parameter free GNM. Computational experiments, conducted based on one training dataset, two independent datasets and one additional small set derived by molecular dynamics simulations, demonstrated that the average correlation coefficients of the proposed RSA based parameter-free GNM, called RpfGNM, were significantly increased when compared with the parameter-free GNM. Our empirical results indicated that a variation of the classical GNMs by combining other protein structural properties is an attractive way to improve the quality of flexibility modeling.

The Renormalization Group and Its Applications to Generating Coarse-Grained Models of Large Biological Molecular Systems

Understanding the dynamics of biomolecules is the key to understanding their biological activities. Computational methods ranging from all-atom molecular dynamics simulations to coarse-grained normal-mode analyses based on simplified elastic networks provide a general framework to studying these dynamics. Despite recent successes in studying very large systems with up to a 100,000,000 atoms, those methods are currently limited to studying small- to medium-sized molecular systems due to computational limitations. One solution to circumvent these limitations is to reduce the size of the system under study. In this paper, we argue that coarse-graining, the standard approach to such size reduction, must define a hierarchy of models of decreasing sizes that are consistent with each other, i.e., that each model contains the information of the dynamics of its predecessor. We propose a new method, Decimate, for generating such a hierarchy within the context of elastic networks for normal-mode analysis. This method is based on the concept of the renormalization group developed in statistical physics. We highlight the details of its implementation, with a special focus on its scalability to large systems of up to millions of atoms. We illustrate its application on two large systems, the capsid of a virus and the ribosome translation complex. We show that highly decimated representations of those systems, containing down to 1% of their original number of atoms, still capture qualitatively and quantitatively their dynamics. Decimate is available as an OpenSource resource.

Comparative Normal Mode Analysis of the Dynamics of DENV and ZIKV Capsids

Key steps in the life cycle of a virus, such as the fusion event as the virus infects a host cell and its maturation process, relate to an intricate interplay between the structure and the dynamics of its constituent proteins, especially those that define its capsid, much akin to an envelope that protects its genomic material. We present a comprehensive, comparative analysis of such interplay for the capsids of two viruses from the flaviviridae family, Dengue (DENV) and Zika (ZIKV). We use for that purpose our own software suite, DD-NMA, which is based on normal mode analysis. We describe the elements of DD-NMA that are relevant to the analysis of large systems, such as virus capsids. In particular, we introduce our implementation of simplified elastic networks and justify their parametrization. Using DD-NMA, we illustrate the importance of packing interactions within the virus capsids on the dynamics of the E proteins of DENV and ZIKV. We identify differences between the computed atomic fluctuations of the E proteins in DENV and ZIKV and relate those differences to changes observed in their high resolution structures. We conclude with a discussion on additional analyses that are needed to fully characterize the dynamics of the two viruses.

Predicting the functional motions of p97 using symmetric normal modes

p97 is a protein complex of the AAA+ family. Although functions of p97 are well understood, the mechanism by which p97 performs its unfolding activities remains unclear. In this work, we present a novel way of applying normal mode analysis to study this six-fold symmetric molecular machine. By selecting normal modes that are axial symmetric and give the largest movements at D1 or D2 pore residues, we are able to predict the functional motions of p97, which are then validated by experimentally observed conformational changes. Our results shed light and provide new understandings on several key steps of the p97 functional process that were previously unclear or controversial, and thus are able to reconcile multiple previous findings. Specifically, our results reveal that (i) a venous valve-like mechanism is used at D2 pore to ensure a one-way exit-only traffic of substrates; (ii) D1 pore remains shut during the functional process; (iii) the "swing-up" motion of the N domain is closely coupled with the vertical motion of the D1 pore along the pore axis; (iv) because of the shut D1 pore and the one-way traffic at D2 pore, it is highly likely that substrates enter the chamber through the gaps at the D1/D2 interface. The limited chamber volume inside p97 suggests that a substrate may be pulling out from D2 while at the same time being pulling in at the interface; (v) lastly, p97 uses a series of actions that alternate between twisting and pulling to remove the substrate. Proteins 2016; 84:1823-1835. © 2016 Wiley Periodicals, Inc.

Theory of Acoustic Raman Modes in Proteins

We present a theoretical analysis that associates the resonances of extraordinary acoustic Raman (EAR) spectroscopy [Wheaton et al., Nat. Photonics 9, 68 (2015)] with the collective modes of proteins. The theory uses the anisotropic elastic network model to find the protein acoustic modes, and calculates Raman intensity by treating the protein as a polarizable ellipsoid. Reasonable agreement is found between EAR spectra and our theory. Protein acoustic modes have been extensively studied theoretically to assess the role they play in protein function; this result suggests EAR spectroscopy as a new experimental tool for studies of protein acoustic modes.

Biophysical and computational methods to analyze amino acid interaction networks in proteins

Globular proteins are held together by interacting networks of amino acid residues. A number of different structural and computational methods have been developed to interrogate these amino acid networks. In this review, we describe some of these methods, including analyses of X-ray crystallographic data and structures, computer simulations, NMR data, and covariation among protein sequences, and indicate the critical insights that such methods provide into protein function. This information can be leveraged towards the design of new allosteric drugs, and the engineering of new protein function and protein regulation strategies.