Towards a free energy-based elastic network model and its application to the SARS-COV2 binding to ACE2.

Classical normal mode analysis (cNMA) is a standard method for studying the equilibrium vibrations of macromolecules. A major limitation of cNMA is that it requires a cumbersome step of energy minimization that also alters the input structure significantly. Variants of normal mode analysis (NMA) exist that perform NMA directly on PDB structures without energy minimization, while maintaining most of the accuracy of cNMA. Spring-based NMA (sbNMA) is such a model. sbNMA uses an all-atom force field as cNMA does, which includes bonded terms such as bond stretching, bond angle bending, torsional, improper, and non-bonded terms such as van der Waals interactions. Electrostatics was not included in sbNMA because it introduced negative spring constants. In this work, we present a way to incorporate most of the electrostatic contributions in normal mode computations, which marks another significant step toward a free-energy-based elastic network model (ENM) for NMA. The vast majority of ENMs are entropy models. One significance of having a free energy-based model for NMA is that it allows one to study the contributions of both entropy and enthalpy. As an application, we apply this model to study the binding stability between SARS-COV2 and angiotensin converting enzyme 2 (or ACE2). Our results show that the stability at the binding interface is contributed nearly equally by hydrophobic interactions and hydrogen bonds.

The amounts of thermal vibrations and static disorder in protein X-ray crystallographic B-factors.

Crystallographic B-factors provide direct dynamical information on the internal mobility of proteins that is closely linked to function, and are also widely used as a benchmark in assessing elastic network models. A significant question in the field is: what is the exact amount of thermal vibrations in protein crystallographic B-factors? This work sets out to answer this question. First, we carry out a thorough, statistically sound analysis of crystallographic B-factors of over 10 000 structures. Second, by employing a highly accurate all-atom model based on the well-known CHARMM force field, we obtain computationally the magnitudes of thermal vibrations of nearly 1000 structures. Our key findings are: (i) the magnitude of thermal vibrations, surprisingly, is nearly protein-independent, as a corollary to the universality for the vibrational spectra of globular proteins established earlier; (ii) the magnitude of thermal vibrations is small, less than 0.1 Å2 at 100 K; (iii) the percentage of thermal vibrations in B-factors is the lowest at low resolution and low temperature (<10%) but increases to as high as 60% for structures determined at high resolution and at room temperature. The significance of this work is that it provides for the first time, using an extremely large dataset, a thorough analysis of B-factors and their thermal and static disorder components. The results clearly demonstrate that structures determined at high resolution and at room temperature have the richest dynamics information. Since such structures are relatively rare in the PDB database, the work naturally calls for more such structures to be determined experimentally.

Effect of solvent on protein structure and dynamics.

Understanding how much solvents influence the structures and dynamics of proteins is important to understand functional mechanisms of solvated proteins. We propose a solvated potential model that approximates the potential energy of a solvated protein by projecting the solvent information into the protein structure. Using the model, we derive three properties of the solvent. First, the influence of the solvent on protein structure and dynamics, mostly by the bulk solvent, decays drastically (near-exponentially) as the distances of the solvent from the protein increase. Using this decay pattern, we suggest the economical size of solvent boxes in molecular dynamics simulations. Second, the hydration shell regulates the protein dynamics by effecting extra interactions within the protein structure. Lastly, the lowest frequency modes are determined mostly by protein structures.

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.

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.

Slow normal modes of proteins are accurately reproduced across different platforms.

The Protein data bank (PDB) (Berman et al 2000 Nucl. Acids Res. 28 235-42) contains the atomic structures of over 105 biomolecules with better than 2.8 Å resolution. The listing of the identities and coordinates of the atoms comprising each macromolecule permits an analysis of the slow-time vibrational response of these large systems to minor perturbations. 3D video animations of individual modes of oscillation demonstrate how regions interdigitate to create cohesive collective motions, providing a comprehensive framework for and familiarity with the overall 3D architecture. Furthermore, the isolation and representation of the softest, slowest deformation coordinates provide opportunities for the development of mechanical models of enzyme function. The eigenvector decomposition, therefore, must be accurate, reliable as well as rapid to be generally reported upon. We obtain the eigenmodes of a 1.2 Å 34 kDa PDB entry using either exclusively heavy atoms or partly or fully reduced atomic sets; Cartesian or internal coordinates; interatomic force fields derived either from a full Cartesian potential, a reduced atomic potential or a Gaussian distance-dependent potential; and independently developed software. These varied technologies are similar in that each maintains proper stereochemistry either by use of dihedral degrees of freedom which freezes bond lengths and bond angles, or by use of a full atomic potential that includes realistic bond length and angle restraints. We find that the shapes of the slowest eigenvectors are nearly identical, not merely similar.

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.

The effective degeneracy of protein normal modes.

Normal modes are frequently computed and used to portray protein dynamics and interpret protein conformational changes. In this work, we investigate the nature of normal modes and find that the normal modes of proteins, especially those at the low frequency range (0-600 cm(-1)), are highly susceptible to degeneracy. Two or more modes are degenerate if they have the same frequency and consequently any orthogonal transformation of them also is a valid representation of the mode subspace. Thus, degenerate modes can no longer characterize unique directions of motions as regular modes do. Though the normal modes of proteins are usually of different frequencies, the difference in frequency between neighboring modes is so small that, under even slight structural uncertainty that unavoidably exists in structure determination, it can easily vanish and as a result, a mode becomes effectively degenerate with its neighboring modes. This can be easily observed in that some modes seem to disappear and their matching modes cannot be found when the structure used to compute the modes is modified only slightly. We term this degeneracy the effective degeneracy of normal modes. This work is built upon our recent discovery that the vibrational spectrum of globular proteins is universal. The high density of modes observed in the vibrational frequency spectra of proteins renders their normal modes highly susceptible to degeneracy, under even the smallest structural uncertainty. Indeed, we find the degree of degeneracy of modes is proportional to the density of modes in the vibrational spectrum. This means that for modes at the same frequency, degeneracy is more severe for larger proteins. Degeneracy exists also in the modes of coarse-grained models, but to a much lesser extent than those of all-atom models. In closing, we discuss the implications of the effective degeneracy of normal modes: how it may significantly affect the ways in which normal modes are used in various normal modes-based applications.

Universality of vibrational spectra of globular proteins.

It is shown that the density of modes of the vibrational spectrum of globular proteins is universal, i.e. regardless of the protein in question, it closely follows one universal curve. The present study, including 135 proteins analyzed with a full atomic empirical potential (CHARMM22) and using the full complement of all atoms Cartesian degrees of freedom, goes far beyond previous claims of universality, confirming that universality holds even in the frequency range that is well above 100 cm(-1) (300-4000 cm(-1)), where peaks and turns in the density of states are faithfully reproduced from one protein to the next. We also characterize fluctuations of the spectral density from the average, paving the way to a meaningful discussion of rare, unusual spectra and the structural reasons for the deviations in such 'outlier' proteins. Since the method used for the derivation of the vibrational modes (potential energy formulation, set of degrees of freedom employed, etc) has a dramatic effect on the spectral density, another significant implication of our findings is that the universality can provide an exquisite tool for assessing and improving the quality of potential functions and the quality of various models used for NMA computations. Finally, we show that the input configuration also affects the density of modes, thus emphasizing the importance of simplified potential energy formulations that are minimized at the outset. In summary, our findings call for a serious two-way dialogue between theory and experiment: experimental spectra of proteins could now guide the fine tuning of theoretical empirical potentials, and the various features and peaks observed in theoretical studies--being universal, and hence now rising in importance--would hopefully spur experimental confirmation.

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

Dynamics can provide deep insights into the functional mechanisms of proteins and protein complexes. For large protein complexes such as GroEL/GroES with more than 8,000 residues, obtaining a fine-grained all-atom description of its normal mode motions can be computationally prohibitive and is often unnecessary. For this reason, coarse-grained models have been used successfully. However, most existing coarse-grained models use extremely simple potentials to represent the interactions within the coarse-grained structures and as a result, the dynamics obtained for the coarse-grained structures may not always be fully realistic. There is a gap between the quality of the dynamics of the coarse-grained structures given by all-atom models and that by coarse-grained models. In this work, we resolve an important question in protein dynamics computations--how can we efficiently construct coarse-grained models whose description of the dynamics of the coarse-grained structures remains as accurate as that given by all-atom models? Our method takes advantage of the sparseness of the Hessian matrix and achieves a high efficiency with a novel iterative matrix projection approach. The result is highly significant since it can provide descriptions of normal mode motions at an all-atom level of accuracy even for the largest biomolecular complexes. The application of our method to GroEL/GroES offers new insights into the mechanism of this biologically important chaperonin, such as that the conformational transitions of this protein complex in its functional cycle are even more strongly connected to the first few lowest frequency modes than with other coarse-grained models.

Conventional NMA as a better standard for evaluating elastic network models.

Normal mode analysis (NMA) is an important tool for studying protein dynamics. Because of the complexity of conventional NMA that uses an all-atom model and a semi-empirical force field, many simplified NMA models have been developed, some of which are known as elastic network models. The quality of these simplified NMA models was assessed mostly by evaluating their predictions against experimental B-factors, and rarely by comparing them with the original NMA. In this work, we take the effort to create a publicly accessible dataset of proteins with their minimized structures, NMA modes, and mean-square fluctuations. Then, for the first time, we evaluate the quality of individual normal modes of several widely used elastic network models by comparing them with the conventional NMA. Our results demonstrate that the conventional NMA presents a better and more complete evaluation measure of the quality of elastic network models. This realization should be very helpful in improving current or designing new, higher quality elastic network models. Moreover, using the conventional NMA as the standard of evaluation, a number of interesting and significant insights into the elastic network models are gained.

Quantitative delineation of how breathing motions open ligand migration channels in myoglobin and its mutants.

Ligand migration and binding are central to the biological functions of many proteins such as myoglobin (Mb) and it is widely thought that protein breathing motions open up ligand channels dynamically. However, how a protein exerts its control over the opening and closing of these channels through its intrinsic dynamics is not fully understood. Specifically, a quantitative delineation of the breathing motions that are needed to open ligand channels is lacking. In this work, we present and apply a novel normal mode-based method to quantitatively delineate what and how breathing motions open ligand migration channels in Mb and its mutants. The motivation behind this work springs from the observation that normal mode motions are closely linked to the breathing motions that are thought to open ligand migration channels. In addition, the method provides a direct and detailed depiction of the motions of each and every residue that lines a channel and can identify key residues that play a dominating role in regulating the channel. The all-atom model and the full force-field employed in the method provide a realistic energetics on the work cost required to open a channel, and as a result, the method can be used to efficiently study the effects of mutations on ligand migration channels and on ligand entry rates. Our results on Mb and its mutants are in excellent agreement with MD simulation results and experimentally determined ligand entry rates.

The performance of fine-grained and coarse-grained elastic network models and its dependence on various factors.

In a recent work we developed a method for deriving accurate simplified models that capture the essentials of conventional all-atom NMA and identified two best simplified models: ssNMA and eANM, both of which have a significantly higher correlation with NMA in mean square fluctuation calculations than existing elastic network models such as ANM and ANMr2, a variant of ANM that uses the inverse of the squared separation distances as spring constants. Here, we examine closely how the performance of these elastic network models depends on various factors, namely, the presence of hydrogen atoms in the model, the quality of input structures, and the effect of crystal packing. The study reveals the strengths and limitations of these models. Our results indicate that ssNMA and eANM are the best fine-grained elastic network models but their performance is sensitive to the quality of input structures. When the quality of input structures is poor, ANMr2 is a good alternative for computing mean-square fluctuations while ANM model is a good alternative for obtaining normal modes.

Bridging between normal mode analysis and elastic network models.

Normal mode analysis (NMA) has been a powerful tool for studying protein dynamics. Elastic network models (ENM), through their simplicity, have made normal mode computations accessible to a much broader research community and for many more biomolecular systems. The drawback of ENMs, however, is that they are less accurate than NMA. In this work, through steps of simplification that starts with NMA and ends with ENMs we build a tight connection between NMA and ENMs. In the process of bridging between the two, we have also discovered several high-quality simplified models. Our best simplified model has a mean correlation with the original NMA that is as high as 0.88. In addition, the model is force-field independent and does not require energy minimization, and thus can be applied directly to experimental structures. Another benefit of drawing the connection is a clearer understanding why ENMs work well and how it can be further improved. We discovered that ANM can be greatly enhanced by including an additional torsional term and a geometry term.

A natural unification of GNM and ANM and the role of inter-residue forces.

The Gaussian network model (GNM) and anisotropic network model (ANM) are two elastic network models that have been widely used to study protein fluctuation dynamics. Both models have strengths and weaknesses. Attempts have been made in the past to unify the two models but they had severe limitations. This work presents a novel theoretical result that shows how GNM and ANM can be unified through taking into account the effect of inter-residue forces. The unified model, called the force spring model, or FSM, is reduced to ANM when all the inter-residue forces are set to zero. Moreover, the unification reveals the role of inter-residue forces in protein fluctuation dynamics. Specifically, the effect of inter-residue forces is closely examined by studying the changes in mean-square fluctuations when inter-residue forces are present.

Generalized spring tensor models for protein fluctuation dynamics and conformation changes.

BACKGROUND: In the last decade, various coarse-grained elastic network models have been developed to study the large-scale motions of proteins and protein complexes where computer simulations using detailed all-atom models are not feasible. Among these models, the Gaussian Network Model (GNM) and Anisotropic Network Model (ANM) have been widely used. Both models have strengths and limitations. GNM can predict the relative magnitudes of protein fluctuations well, but due to its isotropy assumption, it cannot be applied to predict the directions of the fluctuations. In contrast, ANM adds the ability to do the latter, but loses a significant amount of precision in the prediction of the magnitudes. RESULTS: In this book chapter, we present a single model, called generalized spring tensor model (STeM), that is able to predict well both the magnitudes and the directions of the fluctuations. Specifically, STeM performs equally well in B-factor predictions as GNM and has the ability to predict the directions of fluctuations as ANM. This is achieved by employing a physically more realistic potential, the Go-like potential. The potential, which is more sophisticated than that of either GNM or ANM, though adds complexity to the derivation process of the Hessian matrix (which fortunately has been done once for all and the MATLAB code is freely available electronically at http://www.cs.iastate.edu/~gsong/STeM ), causes virtually no performance slowdown. In addition, we show that STeM can be further extended to an all-atom model and protein fluctuation dynamics computed by all-atom STeM matches closely with that by Normal Mode Analysis (NMA). CONCLUSIONS: Derived from a physically more realistic potential, STeM proves to be a natural solution in which advantages that used to exist in two separate models, namely GNM and ANM, are achieved in one single model. It thus lightens the burden to work with two separate models and to relate the modes of GNM with those of ANM at times. By examining the contributions of different interaction terms in the Go potential to the fluctuation dynamics, STeM reveals, (i) a physical explanation for why the distance-dependent, inverse distance square (i.e., 1/r (2)) spring constants perform better than the uniform ones, and (ii), the importance of three-body and four-body interactions to properly modeling protein dynamics.STeM is not limited to coarse-grained protein models that use a single bead, usually the alpha carbon, to represent each residue. The core idea of STeM, deriving the Hessian matrix directly from a physically realistic potential, can be extended to all-atom models as well. We did this and discovered that all-atom STeM model represents a highly close approximation of NMA, yet without the need for energy minimization.

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.

An efficient algorithm calculating common solvent accessible volume.

The solvent accessible surface area and the solvent accessible volume are measurements commonly used in implicit solvent models to include the effect of forces exerted by solvents on the protein surfaces (or the atoms on protein surfaces). The two measurements have limitations in describing interactions between proteins (or proteins' atoms) mediated/bridged by solvents. This is because describing the interactions between proteins should be able to capture the chain of protein-solvent-protein interactions while the solvent accessible surface area or the solvent accessible volume can capture only protein-solvent interactions. If we represent the solvent as a continuous medium, we can consider an atom of a protein can effectively interact with the solvent within a certain distance from its surface (or its own solvent-interacting sphere). In this case, the protein-solvent-protein interactions can be measured by the amount of solvent interacting with two proteins' atoms at the same time (or the volume shared by the two atoms' solvent-interacting spheres excluding the volumes occupied by proteins' atoms). We call the shared volume as the common solvent accessible volume (CSAV); there has been no method developed to determine the CSAV. In this work, we propose a new sweep-line-based method that efficiently calculates the common solvent accessible volume. The performance and accuracy of the proposed sweep-line-based method are compared with those of the naïve voxel-based method. The proposed method takes log-linear time to the number of atoms involved in a CSAV calculation and linear time to the resolution. Our results, tested with 52 protein structures of various sizes, show that the proposed sweep-line-based method is superior to the voxel-based method in both computational efficiency and accuracy.