Identification of ligand binding sites of proteins using the Gaussian Network Model

Dynamically driven correlations in elastic net models reveal sequence of events and causality in proteins

An explicit analytic solution is given for the Langevin equation applied to the Gaussian Network Model of a protein subjected to both a random and a deterministic periodic force. Synchronous and asynchronous components of time correlation functions are derived and an expression for phase differences in the time correlations of residue pairs is obtained. The synchronous component enables the determination of dynamic communities within the protein structure. The asynchronous component reveals causality, where the time correlation function between residues i and j differs depending on whether i is observed before j or vice versa, resulting in directional information flow. Driver and driven residues in the allosteric process of cyclophilin A and human NAD-dependent isocitrate dehydrogenase are determined by a perturbation-scanning technique. Factors affecting phase differences between fluctuations of residues, such as network topology, connectivity, and residue centrality, are identified. Within the constraints of the isotropic Gaussian Network Model, our results show that asynchronicity increases with viscosity and distance between residues, decreases with increasing connectivity, and decreases with increasing levels of eigenvector centrality.

Estimating the Similarity between Protein Pockets

With the exponential increase in publicly available protein structures, the comparison of protein binding sites naturally emerged as a scientific topic to explain observations or generate hypotheses for ligand design, notably to predict ligand selectivity for on- and off-targets, explain polypharmacology, and design target-focused libraries. The current review summarizes the state-of-the-art computational methods applied to pocket detection and comparison as well as structural druggability estimates. The major strengths and weaknesses of current pocket descriptors, alignment methods, and similarity search algorithms are presented. Lastly, an exhaustive survey of both retrospective and prospective applications in diverse medicinal chemistry scenarios illustrates the capability of the existing methods and the hurdle that still needs to be overcome for more accurate predictions.

Prediction of Allosteric Communication Pathways in Proteins

MOTIVATION

Allostery in proteins is an essential phenomenon in biological processes. In this paper, we present a computational model to predict paths of maximum information transfer between active and allosteric sites. In this information theoretic study, we use mutual information as the measure of information transfer, where transition probability of information from one residue to its contacting neighbors is proportional to the magnitude of mutual information between the two residues. Starting from a given residue and using a Hidden Markov Model, we successively determine the neighboring residues that eventually lead to a path of optimum information transfer. The Gaussian approximation of mutual information between residue pairs is adopted. The limits of validity of this approximation are discussed in terms of a nonlinear theory of mutual information and its reduction to the Gaussian form.

RESULTS

Predictions of the model are tested on six widely studied cases, CheY Bacterial Chemotaxis, B-cell Lymphoma extra-large Bcl-xL, Human proline isomerase cyclophilin A (CypA), Dihydrofolate reductase DHFR, HRas GTPase, and Caspase-1. The communication transmission rendering the propagation of local fluctuations from the active sites throughout the structure in multiple paths correlate well with the known experimental data. Distinct paths originating from the active site may likely represent a multi functionality such as involving more than one allosteric site and/or preexistence of some other functional states. Our model is computationally fast and simple, and can give allosteric communication pathways, which are crucial for the understanding and control of protein functionality.

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.

Gaussian network model revisited: Effects of mutation and ligand binding on protein behavior

The coarse-grained Gaussian Network model, GNM, considers only the alpha carbons of the folded protein. Therefore it is not directly applicable to the study of mutation or ligand binding problems where atomic detail is required. This shortcoming is improved by including all atom pairs within the coordination shell of each other into the Kirchoff Adjacency Matrix. Counting all contacts rather than only alpha carbon contacts diminishes the magnitude of fluctuations in the system. But more importantly, it changes the graph-like connectivity structure, i.e., the Kirchoff Adjacency Matrix of the protein. This change depends on amino acid type which introduces amino acid specific and position specific information into the classical coarse-grained GNM which was originally modelled in analogy with the phantom network model of rubber elasticity. With this modification, it is now possible to explain the consequences of mutation and ligand binding on residue fluctuations, their pair-correlations and mutual information (MI) shared by each pair. We refer to the new model as 'all-atom GNM'. Using examples from published data we show that the all-atom GNM gives B-factors that are in better agreement with experiment, can explain effects of mutation on long range communication in PDZ domains and can predict effects of GDP and GTP binding on the dimerization of KRAS.

Synchronous and Asynchronous Response in Dynamically Perturbed Proteins

We present a dynamic perturbation-response model of proteins based on the Gaussian Network Model, where a residue is perturbed periodically, and the dynamic response of other residues is determined. The model shows that periodic perturbation causes a synchronous response in phase with the perturbation and an asynchronous response that is out of phase. The asynchronous component results from the viscous effects of the solvent and other dispersive factors in the system. The model is based on the solution of the Langevin equation in the presence of solvent, noise, and perturbation. We introduce several novel ideas: The concept of storage and loss compliance of the protein and their dependence on structure and frequency; the amount of work lost and the residues that contribute significantly to the lost work; new dynamic correlations that result from perturbation; causality, that is, the response of j when i is perturbed is not equal to the response of i when j is perturbed. As examples, we study two systems, namely, bovine rhodopsin and the class of nanobodies. The general results obtained are (i) synchronous and asynchronous correlations depend strongly on the frequency of perturbation, their magnitude decreases with increasing frequency, (ii) time-delayed mean-squared fluctuations of residues have only synchronous components. Asynchronicity is present only in cross correlations, that is, correlations between different residues, (iii) perturbation of loop residues leads to a large dissipation of work, (iv) correlations satisfy the hypothesis of pre-existing pathways according to which information transfer by perturbation rides on already existing equilibrium correlations in the system, (v) dynamic perturbation can introduce a selective response in the system, where the perturbation of each residue excites different sets of responding residues, and (vi) it is possible to identify nondissipative residues whose perturbation does not lead to dissipation in the protein. Despite its simplicity, the model explains several features of allosteric manipulation.

Comparative effects of oncogenic mutations G12C, G12V, G13D, and Q61H on local conformations and dynamics of K-Ras

K-Ras is the most frequently mutated protein in human cancers. However, until very recently, its oncogenic mutants were viewed as undruggable. To develop inhibitors that directly target oncogenic K-Ras mutants, we need to understand both their mutant-specific and pan-mutant dynamics and conformations. Recently, we have investigated how the most frequently observed K-Ras mutation in cancer patients, G12D, changes its local dynamics and conformations (Vatansever et al., 2019). Here, we extend our analysis to study and compare the local effects of other frequently observed oncogenic mutations, G12C, G12V, G13D and Q61H. For this purpose, we have performed Molecular Dynamics (MD) simulations of each mutant when active (GTP-bound) and inactive (GDP-bound), analyzed their trajectories, and compared how each mutant changes local residue conformations, inter-protein distance distributions, local flexibility and residue pair correlated motions. Our results reveal that in the four active oncogenic mutants we have studied, the α2 helix moves closer to the C-terminal of the α3 helix. However, P-loop mutations cause α3 helix to move away from Loop7, and only G12 mutations change the local conformational state populations of the protein. Furthermore, the motions of coupled residues are mutant-specific: G12 mutations lead to new negative correlations between residue motions, while Q61H destroys them. Overall, our findings on the local conformational states and protein dynamics of oncogenic K-Ras mutants can provide insights for both mutant-selective and pan-mutant targeted inhibition efforts.

Interpreting the Dynamics of Binding Interactions of snRNA and U1A Using a Coarse-Grained Model

The binding interactions of small nuclear RNAs (snRNA) and the associated protein factors are critical to the function of spliceosomes in alternatively splicing primary RNA transcripts. Although molecular dynamics simulations are a powerful tool to interpret the mechanism of biological processes, the atomic-level simulations are, however, too expensive and with limited accuracy for the large-size systems, such as snRNA-protein complexes. We extend the coarse-grained Gaussian network model, which models the RNA-protein complexes as a harmonic chain of Cα, P, and O4' atoms, to investigating the impact of the snRNA-binding interaction on the dynamic stability of the human U1A protein, which is a major component of the spliceosomal U1 small nuclear ribonucleoprotein particle. The results reveal that the first and third loops and the C-terminal helix regions of the U1A domain undergo a significant loss of flexibility upon the RNA binding due to the forming of mostly electrostatic and hydrogen bond interactions with RNA 5' stem and loop. By examining the residues whose mutations significantly change the binding free energy between U1A and snRNA, the Gaussian network model-based calculations show that not only the residues at the binding sites that are traditionally considered to play a major role in U1A-RNA association but also those residues that are far away from the RNA-binding interface can participate in the long-range allosteric signal transmission; these calculations are quantitatively consistent with the data observed in the recent snRNA binding experiments. The study demonstrates a useful avenue to utilize the simplified elastic network model to investigate the dynamics characteristics of the biologically important macromolecular interactions.

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.

An efficient perturbation method to predict the functionally key sites of glutamine binding protein

Glutamine-Binding Protein (GlnBP) of Escherichia coli, an important member of the periplasmic binding protein family, is responsible for the first step in the active transport of glutamine across the cytoplasmic membrane. In this work, the functionally key regulation sites of GlnBP were identified by utilizing a perturbation method proposed by our group, in which the residues whose perturbations markedly change the binding free energy between GlnBP and glutamine are considered to be functionally key residues. The results show that besides the substrate binding sites, some other residues distant from the binding pocket, including the ones in the hinge regions between the two domains, the front- and back- door channels and the exposed region, are important for the function of glutamine binding and transport. The predicted results are well consistent with the theoretical and experimental data, which indicates that our method is an effective approach to identify the key residues important for both ligand binding and long-range allosteric signal transmission. This work can provide some insights into the function performance of GlnBP and the physical mechanism of its allosteric regulation.

Identification of Hot Spots in Protein Structures Using Gaussian Network Model and Gaussian Naive Bayes

Residue fluctuations in protein structures have been shown to be highly associated with various protein functions. Gaussian network model (GNM), a simple representative coarse-grained model, was widely adopted to reveal function-related protein dynamics. We directly utilized the high frequency modes generated by GNM and further performed Gaussian Naive Bayes (GNB) to identify hot spot residues. Two coding schemes about the feature vectors were implemented with varying distance cutoffs for GNM and sliding window sizes for GNB based on tenfold cross validations: one by using only a single high mode and the other by combining multiple modes with the highest frequency. Our proposed methods outperformed the previous work that did not directly utilize the high frequency modes generated by GNM, with regard to overall performance evaluated using F1 measure. Moreover, we found that inclusion of more high frequency modes for a GNB classifier can significantly improve the sensitivity. The present study provided additional valuable insights into the relation between the hot spots and the residue fluctuations.

Explaining the autoinhibition of the SMYD enzyme family: A theoretical study

The SMYD enzymes (SMYD1-5) are lysine methyltransferases that have diverse biological functions including gene expression and regulation of skeletal and cardiac muscle development and function. Recently, they have gained more attention as potential drug targets because of their involvement in cardiovascular diseases and in the progression of different cancer types. Their activity has been suggested to be regulated by a posttranslational mechanism and by autoinhibition. The later relies on a hinge-like movement of the N- and C-lobes to adopt an open or closed conformation, consequently, determining the accessibility of the active site and substrate specificity. In this study we aim to investigate and explain the possibility of the regulatory autoinhibition process of the SMYD enzymes by a thorough computational exploration of their dynamic, energetic, and structural changes by using extended molecular dynamics simulations; normal mode analysis (NMA); and energy correlations. Three SMYD models (SMYD1-3) were used in this study. Our results showed an obvious hinge-like motion between the N- and C-lobes. Also, we identified interaction energy pathways within the 3D structures of the proteins, and hot spots on their surfaces that could be of particular importance for the regulation of their activities via allosteric means. These results can help in a better understanding of the nature of these promising drug targets; and in designing selective drugs that can interfere with (inhibit) the function of a specific SMYD member by disrupting its dynamical and conformational behaviour without disrupting the function of the entire SMYD proteins.

Regulation of ryanodine receptor RyR2 by protein-protein interactions: prediction of a PKA binding site on the N-terminal domain of RyR2 and its relation to disease causing mutations

Protein-protein interactions are the key processes responsible for signaling and function in complex networks. Determining the correct binding partners and predicting the ligand binding sites in the absence of experimental data require predictive models. Hybrid models that combine quantitative atomistic calculations with statistical thermodynamics formulations are valuable tools for bioinformatics predictions. We present a hybrid prediction and analysis model for determining putative binding partners and interpreting the resulting correlations in the yet functionally uncharacterized interactions of the ryanodine RyR2 N-terminal domain. Using extensive docking calculations and libraries of hexameric peptides generated from regulator proteins of the RyR2 channel, we show that the residues 318-323 of protein kinase A, PKA, have a very high affinity for the N-terminal of RyR2. Using a coarse grained Elastic Net Model, we show that the binding site lies at the end of a pathway of evolutionarily conserved residues in RyR2. The two disease causing mutations are also on this path. The program for the prediction of the energetically responsive residues by the Elastic Net Model is freely available on request from the corresponding author.

Sequence-based Gaussian network model for protein dynamics

Motivation: Gaussian network model (GNM) is widely adopted to analyze and understand protein dynamics, function and conformational changes. The existing GNM-based approaches require atomic coordinates of the corresponding protein and cannot be used when only the sequence is known.

Results: We report, first of its kind, GNM model that allows modeling using the sequence. Our linear regression-based, parameter-free, sequence-derived GNM (L-pfSeqGNM) uses contact maps predicted from the sequence and models local, in the sequence, contact neighborhoods with the linear regression. Empirical benchmarking shows relatively high correlations between the native and the predicted with L-pfSeqGNM B-factors and between the cross-correlations of residue fluctuations derived from the structure- and the sequence-based GNM models. Our results demonstrate that L-pfSeqGNM is an attractive platform to explore protein dynamics. In contrast to the highly used GNMs that require protein structures that number in thousands, our model can be used to study motions for the millions of the readily available sequences, which finds applications in modeling conformational changes, protein–protein interactions and protein functions.

Contact:zerozhua@126.com

Supplementary information:Supplementary data are available at Bioinformatics online.

A fast approximate method of identifying paths of allosteric communication in proteins

Fluctuations of the distance between a pair of residues i and j may be correlated with the fluctuations of the distance between another pair k and l. In this case, information may be transmitted among these four residues. Allosteric activity is postulated to proceed through such correlated paths. In this short communication a fast method for calculating correlations among all possible pairs ij and kl leading to a pathway of correlated residues of a protein is proposed. The method is based on the alpha carbon centered Gaussian Network Model. The model is applied to Glutamine Amidotransferase and pathways of allosteric activity are identified and compared with literature.

A quantitative measure of electrostatic perturbation in holo and apo enzymes induced by structural changes

Biological pathways are subject to subtle manipulations that achieve a wide range of functional variation in differing physiological niches. In many instances, changes in the structure of an enzyme on ligand binding germinate electrostatic perturbations that form the basis of its changed catalytic or transcriptional efficiency. Computational methods that seek to gain insights into the electrostatic changes in enzymes require expertise to setup and computing prowess. In the current work, we present a fast, easy and reliable methodology to compute electrostatic perturbations induced by ligand binding (MEPP). The theoretical foundation of MEPP is the conserved electrostatic potential difference (EPD) in cognate pairs of active site residues in proteins with the same functionality. Previously, this invariance has been used to unravel promiscuous serine protease and metallo-β-lactamase scaffolds in alkaline phosphatases. Given that a similarity in EPD is significant, we expect differences in the EPD to be significant too. MEPP identifies residues or domains that undergo significant electrostatic perturbations, and also enumerates residue pairs that undergo significant polarity change. The gain in a certain polarity of a residue with respect to neighboring residues, or the reversal of polarity between two residues might indicate a change in the preferred ligand. The methodology of MEPP has been demonstrated on several enzymes that employ varying mechanisms to perform their roles. For example, we have attributed the change in polarity in residue pairs to be responsible for the loss of metal ion binding in fructose 1,6-bisphosphatases, and corroborated the pre-organized state of the active site of the enzyme with respect to functionally relevant changes in electric fields in ketosteroid isomerases.

Inter-domain movements in polyketide synthases: a molecular dynamics study

Insights into the structure and dynamics of modular polyketide synthases (PKS) are essential for understanding the mechanistic details of the biosynthesis of a large number of pharmaceutically important secondary metabolites. The crystal structures of the KS-AT di-domain from erythromycin synthase have revealed the relative orientation of various catalytic domains in a minimal PKS module. However, the relatively large distance between catalytic centers of KS and AT domains in the static structure has posed certain intriguing questions regarding mechanistic details of substrate transfer during polyketide biosynthesis. In order to investigate the role of inter-domain movements in substrate channeling, we have carried out a series of explicit solvent MD simulations for time periods ranging from 10 to 15 ns on the KS-AT di-domain and its sub-fragments. Analyses of these MD trajectories have revealed that both the catalytic domains and the structured inter-domain linker region remain close to their starting structures. Inter-domain movements at KS-linker and linker-AT interfaces occur around hinge regions which connect the structured linker region to the catalytic domains. The KS-linker interface was found to be more flexible compared to the linker-AT interface. However, inter-domain movements observed during the timescale of our simulations do not significantly reduce the distance between catalytic centers of KS and AT domains for facilitating substrate channeling. Based on these studies and prediction of intrinsic disorder we propose that the intrinsically unstructured linker stretch preceding the ACP domain might be facilitating movement of ACP domains to various catalytic centers.