iGNM 2.0: the Gaussian network model database for biomolecular structural dynamics

Oliveira AS, Mulholland AJ, Serapian SA, Colombo G. Decrypting Allostery in Membrane-Bound K-Ras4B Using Complementary In Silico Approaches Based on Unbiased Molecular Dynamics Simulations

Protein functions are dynamically regulated by allostery, which enables conformational communication even between faraway residues, and expresses itself in many forms, akin to different "languages": allosteric control pathways predominating in an unperturbed protein are often unintuitively reshaped whenever biochemical perturbations arise (e.g., mutations). To accurately model allostery, unbiased molecular dynamics (MD) simulations require integration with a reliable method able to, e.g., detect incipient allosteric changes or likely perturbation pathways; this is because allostery can operate at longer time scales than those accessible by plain MD. Such methods are typically applied singularly, but we here argue their joint application─as a "multilingual" approach─could work significantly better. We successfully prove this through unbiased MD simulations (∼100 μs) of the widely studied, allosterically active oncotarget K-Ras4B, solvated and embedded in a phospholipid membrane, from which we decrypt allostery using four showcase "languages": Distance Fluctuation analysis and the Shortest Path Map capture allosteric hotspots at equilibrium; Anisotropic Thermal Diffusion and Dynamical Non-Equilibrium MD simulations assess perturbations upon, respectively, either superheating or hydrolyzing the GTP that oncogenically activates K-Ras4B. Chosen "languages" work synergistically, providing an articulate, mutually coherent, experimentally consistent picture of K-Ras4B allostery, whereby distinct traits emerge at equilibrium and upon GTP cleavage. At equilibrium, combined evidence confirms prominent allosteric communication from the membrane-embedded hypervariable region, through a hub comprising helix α5 and sheet β5, and up to the active site, encompassing allosteric "switches" I and II (marginally), and two proposed pockets. Upon GTP cleavage, allosteric perturbations mostly accumulate on the switches and documented interfaces.

Modeling protein-nucleic acid complexes with extremely large conformational changes using Flex-LZerD

Proteins and nucleic acids are key components in many processes in living cells, and interactions between proteins and nucleic acids are often crucial pathway components. In many cases, large flexibility of proteins as they interact with nucleic acids is key to their function. To understand the mechanisms of these processes, it is necessary to consider the 3D atomic structures of such protein-nucleic acid complexes. When such structures are not yet experimentally determined, protein docking can be used to computationally generate useful structure models. However, such docking has long had the limitation that the consideration of flexibility is usually limited to small movements or to small structures. We previously developed a method of flexible protein docking which could model ordered proteins which undergo large-scale conformational changes, which we also showed was compatible with nucleic acids. Here, we elaborate on the ability of that pipeline, Flex-LZerD, to model specifically interactions between proteins and nucleic acids, and demonstrate that Flex-LZerD can model more interactions and types of conformational change than previously shown.

Functional insight into Cordyceps militaris sugar transporters by structure modeling, network analysis and allosteric regulation

Insights into the structures, functions and dynamics of Cordyceps militaris (C. militaris) sugar transporters are necessary for understanding their versatile metabolic capability for fungal growth. The sequence-function relationship study of 85 C. militaris sugar transporters showed that there is a gap between phylogenetic-based subfamily classification and their functions. Beyond protein sequences, structural modeling and principal component analysis of the structural ensemble revealed the different folds of the Car and Org subfamilies. Performing channel detection and network analysis found that the Alp and Hex subfamilies can be specifically distinguished from others by the betweenness of channel residues. Signature dynamics analysis further suggested that the Hex subfamily demonstrates different dynamics, with high flexibility at the H1 region in TM11. Furthermore, the H1 region as an allosteric site was examined by network parameter calculations that guided allosteric pathways between this region and the channel cavity. Together with gene expression data of C. militaris, e.g., Hex06741 in the Hex subfamily, it was promisingly expressed when sugar utilization was altered. This work demonstrates an in silico framework for investigating C. militaris sugar transporters as an example case study of the allosteric activity of the Hex subfamily and can facilitate sugar transporter engineering design that can further optimize the preferable sugar utilization and fermentation process of C. militaris.

Domain-Based Protein Docking with Extremely Large Conformational Changes

Proteins are key components in many processes in living cells, and physical interactions with other proteins and nucleic acids often form key parts of their functions. In many cases, large flexibility of proteins as they interact is key to their function. To understand the mechanisms of these processes, it is necessary to consider the 3D structures of such protein complexes. When such structures are not yet experimentally determined, protein docking has long been present to computationally generate useful structure models. However, protein docking has long had the limitation that the consideration of flexibility is usually limited to very small movements or very small structures. Methods have been developed which handle minor flexibility via normal mode or other structure sampling, but new methods are required to model ordered proteins which undergo large-scale conformational changes to elucidate their function at the molecular level. Here, we present Flex-LZerD, a framework for docking such complexes. Via partial assembly multidomain docking and an iterative normal mode analysis admitting curvilinear motions, we demonstrate the ability to model the assembly of a variety of protein-protein and protein-nucleic acid complexes.

The effect of dimerization and ligand binding on the dynamics of Kaposi's sarcoma‐associated herpesvirus protease

The Kaposi's sarcoma‐associated herpesvirus protease is essential for virus maturation. This protease functions under allosteric regulation that establishes its enzymatic activity upon dimerization. It exists in equilibrium between an inactive monomeric state and an active, weakly associating, dimeric state that is stabilized upon ligand binding. The dynamics of the protease dimer and its monomer were studied using the Gaussian network model and the anisotropic network model , and its role in mediating the allosteric regulation is demonstrated. We show that the dimer is composed of five dynamical domains. The central domain is formed upon dimerization and composed of helix five of each monomer, in addition to proximal and distal domains of each monomer. Dimerization reduces the mobility of the central domains and increases the mobility of the distal domains, in particular the binding site within them. The three slowest ANM modes of the dimer assist the protease in ligand binding, motion of the conserved Arg142 and Arg143 toward the oxyanion, and reducing the activation barrier for the tetrahedral transition state by stretching the bond that is cleaved by the protease. In addition, we show that ligand binding reduces the motion of helices α1 and α5 at the interface and explain how ligand binding can stabilize the dimer.

Helical structure motifs made searchable for functional peptide design

The systematic design of functional peptides has technological and therapeutic applications. However, there is a need for pattern-based search engines that help locate desired functional motifs in primary sequences regardless of their evolutionary conservation. Existing databases such as The Protein Secondary Structure database (PSS) no longer serves the community, while the Dictionary of Protein Secondary Structure (DSSP) annotates the secondary structures when tertiary structures of proteins are provided. Here, we extract 1.7 million helices from the PDB and compile them into a database (Therapeutic Peptide Design database; TP-DB) that allows queries of compounded patterns to facilitate the identification of sequence motifs of helical structures. We show how TP-DB helps us identify a known purification-tag-specific antibody that can be repurposed into a diagnostic kit for Helicobacter pylori. We also show how the database can be used to design a new antimicrobial peptide that shows better Candida albicans clearance and lower hemolysis than its template homologs. Finally, we demonstrate how TP-DB can suggest point mutations in helical peptide blockers to prevent a targeted tumorigenic protein-protein interaction. TP-DB is made available at http://dyn.life.nthu.edu.tw/design/ .

Extracting the Dynamic Motion of Proteins Using Normal Mode Analysis

Normal mode analysis (NMA) is a technique for describing the conformational states accessible to a protein in a minimum energy conformation. NMA gives results similar to those produced by principal components analysis of a molecular dynamics simulation, but with only a fraction of the computational effort. Here, we provide a brief overview of the theory and describe three methods for carrying out NMA, including the use of one of the on-line services, the use of off-line software for calculating the projection of the modes calculated from one conformation onto another, and an all-atom NMA calculated using GROMACS. For all three methods, we will use the E1·2Ca²⁺ form of the Ca²⁺-ATPase as a concrete example.

Molecular interactions between C9ORF72 and SMCR8: A local energetic frustration perspective

The hexanucleotide repeat expansion in C9orf72 represents a major cause of amyotrophic lateral sclerosis (ALS) and frontotemporal dementia (FTD). C9orf72, together with SMCR8 and WDR41, can form a stable complex that regulates autophagy and membrane trafficking. Very recently, the cryo-EM structure of C9orf72-SMCR8-WDR41 helps in understanding the structure-function relationship of C9orf72. This protein complex is indispensable to several cellular processes and is strongly linked to familial ALS and FTD. Understanding the molecular basis of the C9orf72-SMCR8 protein-protein interaction is thus important to comprehend their function. To establish a basis for understanding the relationships between sequence, structure, and function of the C9orf72, this study reports a local frustration analysis on the C9orf72-SMCR8 complex structure. An analysis of local frustration profiles indicated that (1) the structural domains in C9orf72 are minimally-frustrated and relatively conserved, (2) high frustration patches on the protein-protein interface (3) increased frustration in the C-terminal helices involved in the dimerization of C9orf72 structures.

Structural analysis of mammalian protein phosphorylation at a proteome level

Phosphorylation is an essential post-translational modification for almost all cellular processes. Several global phosphoproteomics analyses have revealed phosphorylation profiles under different conditions. Beyond identification of phospho-sites, protein structures add another layer of information about their functionality. In this study, we systematically characterize phospho-sites based on their 3D locations in the protein and establish a location map for phospho-sites. More than 250,000 phospho-sites have been analyzed, of which 8,686 sites match at least one structure and are stratified based on their respective 3D positions. Core phospho-sites possess two distinct groups based on their dynamicity. Dynamic core phosphorylations are significantly more functional compared with static ones. The dynamic core and the interface phospho-sites are the most functional among all 3D phosphorylation groups. Our analysis provides global characterization and stratification of phospho-sites from a structural perspective that can be utilized for predicting functional relevance and filtering out false positives in phosphoproteomic studies.

Normal mode analysis of membrane protein dynamics using the vibrational subsystem analysis

The vibrational subsystem analysis is a useful approach that allows for evaluating the spectrum of modes of a given system by integrating out the degrees of freedom accessible to the environment. The approach could be utilized for exploring the collective dynamics of a membrane protein (system) coupled to the lipid bilayer (environment). However, the application to membrane proteins is limited due to high computational costs of modeling a sufficiently large membrane environment unbiased by end effects, which drastically increases the size of the investigated system. We derived a recursive formula for calculating the reduced Hessian of a membrane protein embedded in a lipid bilayer by decomposing the membrane into concentric cylindrical domains with the protein located at the center. The approach allows for the design of a time- and memory-efficient algorithm and a mathematical understanding of the convergence of the reduced Hessian with respect to increasing membrane sizes. The application to the archaeal aspartate transporter Glt_Ph illustrates its utility and efficiency in capturing the transporter's elevator-like movement during its transition between outward-facing and inward-facing states.

Predicting cryptic ligand binding sites based on normal modes guided conformational sampling

To greatly expand the druggable genome, fast and accurate predictions of cryptic sites for small molecules binding in target proteins are in high demand. In this study, we have developed a fast and simple conformational sampling scheme guided by normal modes solved from the coarse-grained elastic models followed by atomistic backbone refinement and side-chain repacking. Despite the observations of complex and diverse conformational changes associated with ligand binding, we found that simply sampling along each of the lowest 30 modes is near optimal for adequately restructuring cryptic sites so they can be detected by existing pocket finding programs like fpocket and concavity. We further trained machine-learning protocols to optimize the combination of the sampling-enhanced pocket scores with other dynamic and conservation scores, which only slightly improved the performance. As assessed based on a training set of 84 known cryptic sites and a test set of 14 proteins, our method achieved high accuracy of prediction (with area under the receiver operating characteristic curve >0.8) comparable to the CryptoSite server. Compared with CryptoSite and other methods based on extensive molecular dynamics simulation, our method is much faster (1-2 hours for an average-size protein) and simpler (using only pocket scores), so it is suitable for high-throughput processing of large datasets of protein structures at the genome scale.

Equally Weighted Multiscale Elastic Network Model and Its Comparison with Traditional and Parameter-Free Models

Dynamical properties of proteins play an essential role in their function exertion. The elastic network model (ENM) is an effective and efficient tool in characterizing the intrinsic dynamical properties encoded in biomacromolecule structures. The Gaussian network model (GNM) and anisotropic network model (ANM) are the two often-used ENM models. Here, we introduce an equally weighted multiscale ENM (equally weighted mENM) based on the original mENM (denoted as mENM), in which fitting weights of Kirchhoff/Hessian matrixes in mENM are removed since they neglect the details of pairwise interactions. Then, we perform its comparison with the mENM, traditional ENM, and parameter-free ENM (pfENM) in reproducing dynamical properties for the six representative proteins whose molecular dynamics (MD) trajectories are available in http://mmb.pcb.ub.es/MoDEL/. In the results, for B-factor prediction, mENM performs best, while the equally weighted mENM performs also well, better than the traditional ENM and pfENM models. As to the dynamical cross-correlation map calculation, mENM performs worst, while the results produced from the equally weighted mENM and pfENM models are close to those from MD trajectories with the latter a little better than the former. Furthermore, encouragingly, the equally weighted mANM displays the best performance in capturing the functional motional modes, followed by pfANM and traditional ANM models, while the mANM fails in all the cases. This work is helpful for strengthening the understanding of the elastic network model and provides a valuable guide for researchers to utilize the model to explore protein dynamics.

Identification of Allosteric Effects in Proteins by Elastic Network Models

Allostery is a fundamental regulatory mechanism in the majority of biological processes of molecular machines. Allostery is well-known as a dynamic-driven process, and thus, the molecular mechanism of allosteric signal transmission needs to be established. Elastic network models (ENMs) provide efficient methods for investigating the intrinsic dynamics and allosteric communication pathways in proteins. In this chapter, two ENM methods including Gaussian network model (GNM) coupled with Markovian stochastic model, as well as the anisotropic network model (ANM), were introduced to identify allosteric effects in hemoglobins. Techniques on model parameters, scripting and calculation, analysis, and visualization are shown step by step.

Resolving the paradox of ferroptotic cell death: Ferrostatin-1 binds to 15LOX/PEBP1 complex, suppresses generation of peroxidized ETE-PE, and protects against ferroptosis

Hydroperoxy-eicosatetraenoyl-phosphatidylethanolamine (HpETE-PE) is a ferroptotic cell death signal. HpETE-PE is produced by the 15-Lipoxygenase (15LOX)/Phosphatidylethanolamine Binding Protein-1 (PEBP1) complex or via an Fe-catalyzed non-enzymatic radical reaction. Ferrostatin-1 (Fer-1), a common ferroptosis inhibitor, is a lipophilic radical scavenger but a poor 15LOX inhibitor arguing against 15LOX having a role in ferroptosis. In the current work, we demonstrate that Fer-1 does not affect 15LOX alone, however, it effectively inhibits HpETE-PE production by the 15LOX/PEBP1 complex. Computational molecular modeling shows that Fer-1 binds to the 15LOX/PEBP1 complex at three sites and could disrupt the catalytically required allosteric motions of the 15LOX/PEBP1 complex. Using nine ferroptosis cell/tissue models, we show that HpETE-PE is produced by the 15LOX/PEBP1 complex and resolve the long-existing Fer-1 anti-ferroptotic paradox.

Distinguishing Enzymes and Non-enzymes Based on Structural Information with an Alignment Free Approach

Background:

Knowledge of protein functions is very crucial for the understanding of biological processes. Experimental methods for protein function prediction are powerless to treat the growing amount of protein sequence and structure data.

Objective:

To develop some computational techniques for the protein function prediction.

Method:

Based on the residue interaction network features and the motion mode information, an SVM model was constructed and used as the predictor. The role of these features was analyzed and some interesting results were obtained.

Results:

An alignment-free method for the classification of enzyme and non-enzyme is developed in this work. There is not any single feature that occupies a dominant position in the prediction process. The topological and the information-theoretic residue interaction network features have a better performance. The combination of the fast mode and the slow mode can get a better explanation for the classification result.

Conclusion:

The method proposed in this paper can act as a classifier for the enzymes and nonenzymes.

Engineering glutathione S-transferase with a point mutation at conserved F136 residue increases the xenobiotic-metabolizing activity

Glutathione S-transferases (GSTs) are multifunctional enzymes that play major roles in a wide range of biological processes, including cellular detoxification, biosynthesis, metabolism, and transport. The dynamic structural scaffold and diverse functional roles of GSTs make them important for enzyme engineering and for exploring novel biotechnological applications. The present study reported a significant gain-of-function activity in GST caused by a point mutation at the conserved F136 residue. The fluorescence quenching and kinetic data suggested that both binding affinity and catalytic efficiency of the mutant enzyme to the substrates 1-chloro-2,4-dinitrobenzene (CDNB), as well as the glutathione (GSH), is increased. Molecular docking showed that the mutation improves the binding interactions of the GSH with several binding-site residues. The simulation of molecular dynamics revealed that the mutant enzyme gained increased structural rigidity than the wild-type enzyme. The mutation also altered the residue interaction network (RIN) of the GSH-binding residues. These phenomena suggested that mutations led to conformational alterations and dominant differential motions in the enzyme that lead to increased rigidity and modifications in RIN. Collectively, engineering GST with a single point mutation at conserved F136 can significantly increase its xenobiotic activity by increasing the catalytic efficiency that may be exploited for biotechnological applications.

Rhapsody: predicting the pathogenicity of human missense variants

MOTIVATION

The biological effects of human missense variants have been studied experimentally for decades but predicting their effects in clinical molecular diagnostics remains challenging. Available computational tools are usually based on the analysis of sequence conservation and structural properties of the mutant protein. We recently introduced a new machine learning method that demonstrated for the first time the significance of protein dynamics in determining the pathogenicity of missense variants.

RESULTS

Here, we present a new interface (Rhapsody) that enables fully automated assessment of pathogenicity, incorporating both sequence coevolution data and structure- and dynamics-based features. Benchmarked against a dataset of about 20 000 annotated variants, the methodology is shown to outperform well-established and/or advanced prediction tools. We illustrate the utility of Rhapsody by in silico saturation mutagenesis studies of human H-Ras, phosphatase and tensin homolog and thiopurine S-methyltransferase.

AVAILABILITY AND IMPLEMENTATION

The new tool is available both as an online webserver at http://rhapsody.csb.pitt.edu and as an open-source Python package (GitHub repository: https://github.com/prody/rhapsody; PyPI package installation: pip install prody-rhapsody). Links to additional resources, tutorials and package documentation are provided in the 'Python package' section of the website.

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.

Maximizing cohesion and separation for detecting protein functional modules in protein-protein interaction networks

Protein Function Module (PFM) identification in Protein-Protein Interaction Networks (PPINs) is one of the most important and challenging tasks in computational biology. The quick and accurate detection of PFMs in PPINs can contribute greatly to the understanding of the functions, properties, and biological mechanisms in research on various diseases and the development of new medicines. Despite the performance of existing detection approaches being improved to some extent, there are still opportunities for further enhancements in the efficiency, accuracy, and robustness of such detection methods. Based on the uniqueness of the network-clustering problem in the context of PPINs, this study proposed a very effective and efficient model based on the Lin-Kernighan-Helsgaun algorithm for detecting PFMs in PPINs. To demonstrate the effectiveness and efficiency of the proposed model, computational experiments are performed using three different categories of species datasets. The computational results reveal that the proposed model outperforms existing detection techniques in terms of two key performance indices, i.e., the degree of polymerization inside PFMs (cohesion) and the deviation degree between PFMs (separation), while being very fast and robust. The proposed model can be used to help researchers decide whether to conduct further expensive and time-consuming biological experiments and to select target proteins from large-scale PPI data for further detailed research.

Allosteric regulation of glutamate dehydrogenase deamination activity

Glutamate dehydrogenase (GDH) is a key enzyme interlinking carbon and nitrogen metabolism. Recent discoveries of the GDH specific role in breast cancer, hyperinsulinism/hyperammonemia (HI/HA) syndrome, and neurodegenerative diseases have reinvigorated interest on GDH regulation, which remains poorly understood despite extensive and long standing studies. Notwithstanding the growing evidence of the complexity of allosteric network behind GDH regulation, identifications of allosteric factors and associated mechanisms are paramount to deepen our understanding of the complex dynamics that regulate GDH enzymatic activity. Combining structural analyses of cryo-electron microscopy data with molecular dynamic simulations, here we show that the cofactor NADH is a key player in the GDH regulation process. Our structural analysis indicates that, binding to the regulatory sites in proximity of the antenna region, NADH acts as a positive allosteric modulator by enhancing both the affinity of the inhibitor GTP binding and inhibition of GDH catalytic activity. We further show that the binding of GTP to the NADH-bound GDH activates a triangular allosteric network, interlinking the inhibitor with regulatory and catalytic sites. This allostery produces a local conformational rearrangement that triggers an anticlockwise rotational motion of interlinked alpha-helices with specific tilted helical extension. This structural transition is a fundamental switch in the GDH enzymatic activity. It introduces a torsional stress, and the associated rotational shift in the Rossmann fold closes the catalytic cleft with consequent inhibition of the deamination process. In silico mutagenesis examinations further underpin the molecular basis of HI/HA dominant mutations and consequent over-activity of GDH through alteration of this allosteric communication network. These results shed new light on GDH regulation and may lay new foundation in the design of allosteric agents.

ProSNEx: a web-based application for exploration and analysis of protein structures using network formalism

ProSNEx (Protein Structure Network Explorer) is a web service for construction and analysis of Protein Structure Networks (PSNs) alongside amino acid flexibility, sequence conservation and annotation features. ProSNEx constructs a PSN by adding nodes to represent residues and edges between these nodes using user-specified interaction distance cutoffs for either carbon-alpha, carbon-beta or atom-pair contact networks. Different types of weighted networks can also be constructed by using either (i) the residue-residue interaction energies in the format returned by gRINN, resulting in a Protein Energy Network (PEN); (ii) the dynamical cross correlations from a coarse-grained Normal Mode Analysis (NMA) of the protein structure; (iii) interaction strength. Upon construction of the network, common network metrics (such as node centralities) as well as shortest paths between nodes and k-cliques are calculated. Moreover, additional features of each residue in the form of conservation scores and mutation/natural variant information are included in the analysis. By this way, tool offers an enhanced and direct comparison of network-based residue metrics with other types of biological information. ProSNEx is free and open to all users without login requirement at http://prosnex-tool.com.

Breaking a single hydrogen bond in the mitochondrial tRNAPhe -PheRS complex leads to phenotypic pleiotropy of human disease

Various pathogenic variants in both mitochondrial tRNA^Phe and Phenylalanyl‐tRNA synthetase mitochondrial protein coding gene (FARS2) gene encoding for the human mitochondrial PheRS have been identified and associated with neurological and/or muscle‐related pathologies. An important Guanine‐34 (G34)A anticodon mutation associated with myoclonic epilepsy with ragged red fibers (MERRF) syndrome has been reported in hmit‐tRNA^Phe. The majority of G34 contacts in available aaRSs‐tRNAs complexes specifically use that base as an important tRNA identity element. The network of intermolecular interactions providing its specific recognition also largely conserved. However, their conservation depends also on the invariance of the residues in the anticodon binding domain (ABD) of human mitochondrial Phenylalanyl‐tRNA synthetase (hmit‐PheRS). A defect in recognition of the anticodon of tRNA^Phe may happen not only because of G34A mutation, but also due to mutations in the ABD. Indeed, a pathogenic mutation in FARS2 has been recently reported in a 9‐year‐old female patient harboring a p.Asp364Gly mutation. Asp364 is hydrogen bonded (HB) to G34 in WT hmit‐PheRS. Thus, there are two pathogenic variants disrupting HB between G34 and Asp364: one is associated with G34A mutation, and the other with Asp364Gly mutation. We have measured the rates of tRNA^Phe aminoacylation catalyzed by WT hmit‐PheRS and mutant enzymes. These data ranked the residues making a HB with G34 according to their contribution to activity and the signal transduction pathway in the hmit‐PheRS‐tRNA^Phe complex. Furthermore, we carried out extensive MD simulations to reveal the interdomain contact topology on the dynamic trajectories of the complex, and gaining insight into the structural and dynamic integrity effects of hmit‐PheRS complexed with tRNA^Phe.

Database

Structural data are available in PDB database under the accession number(s): 3CMQ, 3TUP, 5MGH, 5MGV.

Modulation of Toroidal Proteins Dynamics in Favor of Functional Mechanisms upon Ligand Binding

Toroidal proteins serve as molecular machines and play crucial roles in biological processes such as DNA replication and RNA transcription. Despite progress in the structural characterization of several toroidal proteins, we still lack a mechanistic understanding of the significance of their architecture, oligomerization states, and intermolecular interactions in defining their biological function. In this work, we analyze the collective dynamics of toroidal proteins with different oligomerization states, namely, dimeric and trimeric DNA sliding clamps, nucleocapsid proteins (4-, 5-, and 6-mers) and Trp RNA-binding attenuation proteins (11- and 12-mers). We observe common global modes, among which cooperative rolling stands out as a mechanism enabling DNA processivity, and clamshell motions as those underlying the opening/closure of the sliding clamps. Alterations in global dynamics due to complexation with DNA or the clamp loader are shown to assist in enhancing motions to enable robust function. The analysis provides new insights into the differentiation and enhancement of functional motions upon intersubunit and intermolecular interactions.

Molecular dynamics simulations and linear response theories jointly describe biphasic responses of myoglobin relaxation and reveal evolutionarily conserved frequent communicators

In this study, we provide a time-dependent mechanical model, taking advantage of molecular dynamics simulations, quasiharmonic analysis of molecular dynamics trajectories, and time-dependent linear response theories to describe vibrational energy redistribution within the protein matrix. The theoretical description explained the observed biphasic responses of specific residues in myoglobin to CO-photolysis and photoexcitation on heme. The fast responses were found to be triggered by impulsive forces and propagated mainly by principal modes <40 cm⁻¹. The predicted fast responses for individual atoms were then used to study signal propagation within the protein matrix and signals were found to propagate ~8 times faster across helices (4076 m/s) than within the helices, suggesting the importance of tertiary packing in the sensitivity of proteins to external perturbations. We further developed a method to integrate multiple intramolecular signal pathways and discover frequent “communicators”. These communicators were found to be evolutionarily conserved including those distant from the heme.

Shared Signature Dynamics Tempered by Local Fluctuations Enables Fold Adaptability and Specificity

Recent studies have drawn attention to the evolution of protein dynamics, in addition to sequence and structure, based on the premise structure-encodes-dynamics-encodes-function. Of interest is to understand how functional differentiation is accomplished while maintaining the fold, or how intrinsic dynamics plays out in the evolution of structural variations and functional specificity. We performed a systematic computational analysis of 26,899 proteins belonging to 116 CATH superfamilies. Characterizing cooperative mechanisms and convergent/divergent features that underlie the shared/differentiated dynamics of family members required a methodology that lends itself to efficient analyses of large ensembles of proteins. We therefore introduced, SignDy, an integrated pipeline for evaluating the signature dynamics of families based on elastic network models. Our analysis confirmed that family members share conserved, highly cooperative (global) modes of motion. Importantly, our analysis discloses a subset of motions that sharply distinguishes subfamilies, which lie in a low-to-intermediate frequency regime of the mode spectrum. This regime has maximal impact on functional differentiation of families into subfamilies, while being evolutionarily conserved among subfamily members. Notably, the high-frequency end of the spectrum also reveals evolutionary conserved features across and within subfamilies; but in sharp contrast to global motions, high-frequency modes are minimally collective. Modulation of robust/conserved global dynamics by low-to-intermediate frequency fluctuations thus emerges as a versatile mechanism ensuring the adaptability of selected folds and the specificity of their subfamilies. SignDy further allows for dynamics-based categorization as a new layer of information relevant to distinctive mechanisms of action of subfamilies, beyond sequence or structural classifications.

Essential site scanning analysis: A new approach for detecting sites that modulate the dispersion of protein global motions

Despite the wealth of methods developed for exploring the molecular basis of allostery in biomolecular systems, there is still a need for structure-based predictive tools that can efficiently detect susceptible sites for triggering allosteric responses. Toward this goal, we introduce here an elastic network model (ENM)-based method, Essential Site Scanning Analysis (ESSA). Essential sites are here defined as residues that would significantly alter the protein's global dynamics if bound to a ligand. To mimic the crowding induced upon substrate binding, the heavy atoms of each residue are incorporated as additional network nodes into the α-carbon-based ENM, and the resulting shifts in soft mode frequencies are used as a metric for evaluating the essentiality of each residue. Results on a dataset of monomeric proteins indicate the enrichment of allosteric and orthosteric binding sites, as well as global hinge regions among essential residues, highlighting the significant role of these sites in controlling the overall structural dynamics. Further integration of ESSA with information on predicted pockets and their local hydrophobicity density enables successful predictions of allosteric pockets for both ligand-bound and -unbound structures. ESSA can be efficiently applied to large multimeric systems. Three case studies, namely (i) G-protein binding to a GPCR, (ii) heterotrimeric assembly of the Ser/Thr protein phosphatase PP2A, and (iii) allo-targeting of AMPA receptor, demonstrate the utility of ESSA for identifying essential sites and narrowing down target allosteric sites identified by druggability simulations.

Dynamical comparison between Drosha and Dicer reveals functional motion similarities and dissimilarities

Drosha and Dicer are RNase III family members of classes II and III, respectively, which play a major role in the maturation of micro-RNAs. The two proteins share similar domain arrangement and overall fold despite no apparent sequence homology. The overall structural and catalytic reaction similarity of both proteins, on the one hand, and differences in the substrate and its binding mechanisms, on the other, suggest that both proteins also share dynamic similarities and dissimilarities. Since dynamics is essential for protein function, a comparison at their dynamics level is fundamental for a complete understanding of the overall relations between these proteins. In this study, we present a dynamical comparison between human Drosha and Giardia Dicer. Gaussian Network Model and Anisotropic Network Model modes of motion of the proteins are calculated. Dynamical comparison is performed using global and local dynamic programming algorithms for aligning modes of motion. These algorithms were recently developed based on the commonly used Needleman-Wunsch and Smith-Waterman algorithms for global and local sequence alignment. The slowest mode of Drosha is different from that of Dicer due to its more bended posture and allow the motion of the double-stranded RNA-binding domain toward and away from its substrate. Among the five slowest modes dynamics similarity exists only for the second slow mode of motion of Drosha and Dicer. In addition, high local dynamics similarity is observed at the catalytic domains, in the vicinity of the catalytic residues. The results suggest that the proteins exert a similar catalytic mechanism using similar motions, especially at the catalytic sites.

Pharmmaker: Pharmacophore modeling and hit identification based on druggability simulations

Recent years have seen progress in druggability simulations, that is, molecular dynamics simulations of target proteins in solutions containing drug‐like probe molecules to characterize their drug‐binding abilities, if any. An important consecutive step is to analyze the trajectories to construct pharmacophore models (PMs) to use for virtual screening of libraries of small molecules. While considerable success has been observed in this type of computer‐aided drug discovery, a systematic tool encompassing multiple steps from druggability simulations to pharmacophore modeling, to identifying hits by virtual screening of libraries of compounds, has been lacking. We address this need here by developing a new tool, Pharmmaker, building on the DruGUI module of our ProDy application programming interface. Pharmmaker is composed of a suite of steps: (Step 1) identification of high affinity residues for each probe molecule type; (Step 2) selecting high affinity residues and hot spots in the vicinity of sites identified by DruGUI; (Step 3) ranking of the interactions between high affinity residues and specific probes; (Step 4) obtaining probe binding poses and corresponding protein conformations by collecting top‐ranked snapshots; and (Step 5) using those snapshots for constructing PMs. The PMs are then used as filters for identifying hits in structure‐based virtual screening. Pharmmaker, accessible online at http://prody.csb.pitt.edu/pharmmaker/, can be used in conjunction with other tools available in ProDy.

LARMD: integration of bioinformatic resources to profile ligand-driven protein dynamics with a case on the activation of estrogen receptor

Protein dynamics is central to all biological processes, including signal transduction, cellular regulation and biological catalysis. Among them, in-depth exploration of ligand-driven protein dynamics contributes to an optimal understanding of protein function, which is particularly relevant to drug discovery. Hence, a wide range of computational tools have been designed to investigate the important dynamic information in proteins. However, performing and analyzing protein dynamics is still challenging due to the complicated operation steps, giving rise to great difficulty, especially for nonexperts. Moreover, there is a lack of web protocol to provide online facility to investigate and visualize ligand-driven protein dynamics. To this end, in this study, we integrated several bioinformatic tools to develop a protocol, named Ligand and Receptor Molecular Dynamics (LARMD, http://chemyang.ccnu.edu.cn/ccb/server/LARMD/ and http://agroda.gzu.edu.cn:9999/ccb/server/LARMD/), for profiling ligand-driven protein dynamics. To be specific, estrogen receptor (ER) was used as a case to reveal ERβ-selective mechanism, which plays a vital role in the treatment of inflammatory diseases and many types of cancers in clinical practice. Two different residues (Ile373/Met421 and Met336/Leu384) in the pocket of ERβ/ERα were the significant determinants for selectivity, especially Met336 of ERβ. The helix H8, helix H11 and H7-H8 loop influenced the migration of selective agonist (WAY-244). These computational results were consistent with the experimental results. Therefore, LARMD provides a user-friendly online protocol to study the dynamic property of protein and to design new ligand or site-directed mutagenesis.

An Efficient Timer and Sizer of Biomacromolecular Motions

Life ticks as fast as how proteins move. Computationally expensive molecular dynamics simulation has been the only theoretical tool to gauge the time and sizes of these motions, though barely to their slowest ends. Here, we convert a computationally cheap elastic network model (ENM) into a molecular timer and sizer to gauge the slowest functional motions of structured biomolecules. Quasi-harmonic analysis, fluctuation profile matching, and the Wiener-Khintchine theorem are used to define the "time periods," t, for anharmonic principal components (PCs), which are validated by nuclear magnetic resonance (NMR) order parameters. The PCs with their respective "time periods" are mapped to the eigenvalues (λ_ENM) of the corresponding ENM modes. Thus, the power laws t(ns) = 56.1λ_ENM^-1.6 and σ²(Ų) = 32.7λ_ENM^-3.0 can be established allowing the characterization of the timescales of NMR-resolved conformers, crystallographic anisotropic displacement parameters, and important ribosomal motions, as well as motional sizes of the latter.

Normal Mode Analysis as a Routine Part of a Structural Investigation

Normal mode analysis (NMA) is a technique that can be used to describe the flexible states accessible to a protein about an equilibrium position. These states have been shown repeatedly to have functional significance. NMA is probably the least computationally expensive method for studying the dynamics of macromolecules, and advances in computer technology and algorithms for calculating normal modes over the last 20 years have made it nearly trivial for all but the largest systems. Despite this, it is still uncommon for NMA to be used as a component of the analysis of a structural study. In this review, we will describe NMA, outline its advantages and limitations, explain what can and cannot be learned from it, and address some criticisms and concerns that have been voiced about it. We will then review the most commonly used techniques for reducing the computational cost of this method and identify the web services making use of these methods. We will illustrate several of their possible uses with recent examples from the literature. We conclude by recommending that NMA become one of the standard tools employed in any structural study.

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.

DynOmics: dynamics of structural proteome and beyond

DynOmics (dynomics.pitt.edu) is a portal developed to leverage rapidly growing structural proteomics data by efficiently and accurately evaluating the dynamics of structurally resolved systems, from individual molecules to large complexes and assemblies, in the context of their physiological environment. At the core of the portal is a newly developed server, ENM 1.0, which permits users to efficiently generate information on the collective dynamics of any structure in PDB format, user-uploaded or database-retrieved. ENM 1.0 integrates two widely used elastic network models (ENMs)—the Gaussian Network Model (GNM) and the Anisotropic Network Model (ANM), extended to take account of molecular environment. It enables users to assess potentially functional sites, signal transduction or allosteric communication mechanisms, and protein–protein and protein–DNA interaction poses, in addition to delivering ensembles of accessible conformers reconstructed at atomic details based on the global modes of motions predicted by the ANM. The ‘environment’ is defined in a flexible manner, from lipid bilayer and crystal contacts, to substrate or ligands bound to a protein, or surrounding subunits in a multimeric structure or assembly. User-friendly interactive features permit users to easily visualize how the environment alter the intrinsic dynamics of the query systems. ENM 1.0 can be accessed at http://enm.pitt.edu/ or http://dyn.life.nthu.edu.tw/oENM/.

Towards the functional high-resolution coordination chemistry of blood plasma human serum albumin

Human serum albumin (HSA) is a monomeric, globular, multi-carrier and the most abundant protein in the blood. HSA displays multiple ligand binding sites with extraordinary binding capacity for a wide range of ions and molecules. For decades, HSA's ability to bind to various ligands has led many scientists to study its physiological properties and protein structure; indeed, a better understanding of HSA-ligand interactions in human blood, at the atomic level, will likely foster the development of more potent, and overall more performant, diagnostic and therapeutic tools against serious human disorders such as diabetes, cardiovascular disorders, and cancer. Here, we present a concise overview of the current knowledge of HSA's structural characteristics, and its coordination chemistry with transition metal ions, within the scope and limitations of current techniques and biophysical methods to reach atomic resolution in solution and in blood serum. We also highlight the overwhelming need of a detailed atomistic understanding of HSA dynamic structures and interactions that are transient, weak, multi-site and multi-step, and allosterically affected by each other. Considering the fact that HSA is a current clinical tool for drug delivery systems and a potential contender as molecular cargo and nano-vehicle used in biophysical, clinical and industrial fields, we underline the emerging need for novel approaches to target the dynamic functional coordination chemistry of the human blood serum albumin in solution, at the atomic level.

The theory of frame ordering: observing motions in calmodulin complexes

Large scale functional motions of molecules are studied experimentally using numerous molecular and biophysics techniques, the data from which are subsequently interpreted using diverse models of Brownian molecular dynamics. To unify all rotational physics techniques and motional models, the frame order tensor - a universal statistical mechanics theory based on the rotational ordering of rigid body frames - is herein formulated. The frame ordering is the fundamental physics that governs how motions modulate rotational molecular physics and it defines the properties and maximum information content encoded in the observable physics. Using the tensor to link residual dipolar couplings and pseudo-contact shifts, two distinct information-rich and atomic-level biophysical measurements from the field of nuclear magnetic resonance spectroscopy, to a number of basic mechanical joint models, a highly dynamic state of calmodulin (CaM) bound to a target peptide in a tightly closed conformation was observed. Intra- and inter-domain motions reveal the CaM complex to be entropically primed for peptide release.

Characterization of Differential Dynamics, Specificity, and Allostery of Lipoxygenase Family Members

Accurate modeling of structural dynamics of proteins and their differentiation across different species can help us understand generic mechanisms of function shared by family members and the molecular basis of the specificity of individual members. We focused here on the family of lipoxygenases, enzymes that catalyze lipid oxidation, the mammalian and bacterial structures of which have been elucidated. We present a systematic method of approach for characterizing the sequence, structure, dynamics, and allosteric signaling properties of these enzymes using a combination of structure-based models and methods and bioinformatics tools applied to a data set of 88 structures. The analysis elucidates the signature dynamics of the lipoxygenase family and its differentiation among members, as well as key sites that enable its adaptation to specific substrate binding and allosteric activity.

Domain insertion permissibility-guided engineering of allostery in ion channels

Allostery is a fundamental principle of protein regulation that remains hard to engineer, particularly in membrane proteins such as ion channels. Here we use human Inward Rectifier K⁺ Channel Kir2.1 to map site-specific permissibility to the insertion of domains with different biophysical properties. We find that permissibility is best explained by dynamic protein properties, such as conformational flexibility. Several regions in Kir2.1 that are equivalent to those regulated in homologs, such as G-protein-gated inward rectifier K⁺ channels (GIRK), have differential permissibility; that is, for these sites permissibility depends on the structural properties of the inserted domain. Our data and the well-established link between protein dynamics and allostery led us to propose that differential permissibility is a metric of latent allosteric capacity in Kir2.1. In support of this notion, inserting light-switchable domains into sites with predicted latent allosteric capacity renders Kir2.1 activity sensitive to light.

Allostery is a fundamental principle of protein regulation that remains challenging to engineer. Here authors screen human Inward Rectifier K + Channel Kir2.1 for permissibility to domain insertions and propose that differential permissibility is a metric of latent allosteric capacity in Kir2.1.

Progress in Allosteric Database

An allosteric mechanism refers to the biological regulation process wherein macromolecules propagate the effect of ligand binding at one site to a spatially distant orthosteric locus, thus affecting activity. The theory has remained a trending topic in biology research for over 50 years, since the understanding of allostery is fundamental for gleaning numerous biological processes and developing new drug therapies. In the past two decades, the allosteric paradigm has evolved into more descriptive models, with ever-expanding amounts of experimental data pertaining to newly identified allosteric molecules. The AlloSteric Database (ASD, accessible at http://mdl.shsmu.edu.cn/ASD ), which is a comprehensive knowledge repository, has provided the public with integrated information encompassing allosteric proteins, modulators, sites, pathways, and networks to investigate allostery since 2009. In this chapter, we introduce the history and usage of the ASD and give attention to specific applications that have benefited from the ASD.

Coarse-Grained Protein Dynamics Studies Using Elastic Network Models

Elastic networks have been used as simple models of proteins to study their slow structural dynamics. They consist of point-like particles connected by linear Hookean springs and hence are convenient for linear normal mode analysis around a given reference structure. Furthermore, dynamic simulations using these models can provide new insights. As the computational cost associated with these models is considerably lower compared to that of all-atom models, they are also convenient for comparative studies between multiple protein structures. In this review, we introduce examples of coarse-grained molecular dynamics studies using elastic network models and their derivatives, focusing on the nonlinear phenomena, and discuss their applicability to large-scale macromolecular assemblies.

Modeling of Protein Structural Flexibility and Large-Scale Dynamics: Coarse-Grained Simulations and Elastic Network Models

Fluctuations of protein three-dimensional structures and large-scale conformational transitions are crucial for the biological function of proteins and their complexes. Experimental studies of such phenomena remain very challenging and therefore molecular modeling can be a good alternative or a valuable supporting tool for the investigation of large molecular systems and long-time events. In this minireview, we present two alternative approaches to the coarse-grained (CG) modeling of dynamic properties of protein systems. We discuss two CG representations of polypeptide chains used for Monte Carlo dynamics simulations of protein local dynamics and conformational transitions, and highly simplified structure-based elastic network models of protein flexibility. In contrast to classical all-atom molecular dynamics, the modeling strategies discussed here allow the quite accurate modeling of much larger systems and longer-time dynamic phenomena. We briefly describe the main features of these models and outline some of their applications, including modeling of near-native structure fluctuations, sampling of large regions of the protein conformational space, or possible support for the structure prediction of large proteins and their complexes.

Dynamical comparison between myoglobin and hemoglobin

Myoglobin and hemoglobin are globular hemeproteins, when the former is a monomer and the latter a heterotetramer. Despite the structural similarity of myoglobin to α and β subunits of hemoglobin, there is a functional difference between the two proteins, owing to the quaternary structure of hemoglobin. The effect of the quaternary structure of hemoglobin on the intrinsic dynamics of its subunits is explored by dynamical comparison of the two proteins. Anisotropic Network Model modes of motion were calculated for hemoglobin and myoglobin. Dynamical comparison between the proteins was performed using global and local Anisotropic Network Model mode alignment algorithms based on the algorithms of Smith-Waterman and Needleman-Wunsch for sequence comparison. The results indicate that the quaternary structure of hemoglobin substantially alters the intrinsic dynamics of its subunits, an effect that may contribute to the functional difference between the two proteins. Local dynamics similarity between the proteins is still observed at the major exit route of the ligand.

Pseudomonas aeruginosa utilizes host polyunsaturated phosphatidylethanolamines to trigger theft-ferroptosis in bronchial epithelium

Ferroptosis is a death program executed via selective oxidation of arachidonic acid-phosphatidylethanolamines (AA-PE) by 15-lipoxygenases. In mammalian cells and tissues, ferroptosis has been pathogenically associated with brain, kidney, and liver injury/diseases. We discovered that a prokaryotic bacterium, Pseudomonas aeruginosa, that does not contain AA-PE can express lipoxygenase (pLoxA), oxidize host AA-PE to 15-hydroperoxy-AA-PE (15-HOO-AA-PE), and trigger ferroptosis in human bronchial epithelial cells. Induction of ferroptosis by clinical P. aeruginosa isolates from patients with persistent lower respiratory tract infections was dependent on the level and enzymatic activity of pLoxA. Redox phospholipidomics revealed elevated levels of oxidized AA-PE in airway tissues from patients with cystic fibrosis (CF) but not with emphysema or CF without P. aeruginosa. We believe that the evolutionarily conserved mechanism of pLoxA-driven ferroptosis may represent a potential therapeutic target against P. aeruginosa-associated diseases such as CF and persistent lower respiratory tract infections.

Energy transport pathway in proteins: Insights from non-equilibrium molecular dynamics with elastic network model

Intra-molecular energy transport between distant functional sites plays important roles in allosterically regulating the biochemical activity of proteins. How to identify the specific intra-molecular signaling pathway from protein tertiary structure remains a challenging problem. In the present work, a non-equilibrium dynamics method based on the elastic network model (ENM) was proposed to simulate the energy propagation process and identify the specific signaling pathways within proteins. In this method, a given residue was perturbed and the propagation of energy was simulated by non-equilibrium dynamics in the normal modes space of ENM. After that, the simulation results were transformed from the normal modes space to the Cartesian coordinate space to identify the intra-protein energy transduction pathways. The proposed method was applied to myosin and the third PDZ domain (PDZ3) of PSD-95 as case studies. For myosin, two signaling pathways were identified, which mediate the energy transductions form the nucleotide binding site to the 50 kDa cleft and the converter subdomain, respectively. For PDZ3, one specific signaling pathway was identified, through which the intra-protein energy was transduced from ligand binding site to the distant opposite side of the protein. It is also found that comparing with the commonly used cross-correlation analysis method, the proposed method can identify the anisotropic energy transduction pathways more effectively.

Protein dynamic communities from elastic network models align closely to the communities defined by molecular dynamics

Dynamic communities in proteins comprise the cohesive structural units that individually exhibit rigid body motions. These can correspond to structural domains, but are usually smaller parts that move with respect to one another in a protein's internal motions, key to its functional dynamics. Previous studies emphasized their importance to understand the nature of ligand-induced allosteric regulation. These studies reported that mutations to key community residues can hinder transmission of allosteric signals among the communities. Usually molecular dynamic (MD) simulations (~ 100 ns or longer) have been used to identify the communities-a demanding task for larger proteins. In the present study, we propose that dynamic communities obtained from MD simulations can also be obtained alternatively with simpler models-the elastic network models (ENMs). To verify this premise, we compare the specific communities obtained from MD and ENMs for 44 proteins. We evaluate the correspondence in communities from the two methods and compute the extent of agreement in the dynamic cross-correlation data used for community detection. Our study reveals a strong correspondence between the communities from MD and ENM and also good agreement for the residue cross-correlations. Importantly, we observe that the dynamic communities from MD can be closely reproduced with ENMs. With ENMs, we also compare the community structures of stable and unstable mutant forms of T4 Lysozyme with its wild-type. We find that communities for unstable mutants show substantially poorer agreement with the wild-type communities than do stable mutants, suggesting such ENM-based community structures can serve as a means to rapidly identify deleterious mutants.

Dependence of prevalence of contiguous pathways in proteins on structural complexity

Allostery is a regulatory mechanism in proteins where an effector molecule binds distal from an active site to modulate its activity. Allosteric signaling may occur via a continuous path of residues linking the active and allosteric sites, which has been suggested by large conformational changes evident in crystal structures. An alternate possibility is that the signal occurs in the realm of ensemble dynamics via an energy landscape change. While the latter was first proposed on theoretical grounds, increasing evidence suggests that such a control mechanism is plausible. A major difficulty for testing the two methods is the ability to definitively determine that a residue is directly involved in allosteric signal transduction. Statistical Coupling Analysis (SCA) is a method that has been successful at predicting pathways, and experimental tests involving mutagenesis or domain substitution provide the best available evidence of signaling pathways. However, ascertaining energetic pathways which need not be contiguous is far more difficult. To date, simple estimates of the statistical significance of a pathway in a protein remain to be established. The focus of this work is to estimate such benchmarks for the statistical significance of contiguous pathways for the null model of selecting residues at random. We found that when 20% of residues in proteins are randomly selected, contiguous pathways at the 6 Å cutoff level were found with success rates of 51% in PDZ, 30% in p53, and 3% in MutS. The results suggest that the significance of pathways may have system specific factors involved. Furthermore, the possible existence of false positives for contiguous pathways implies that signaling could be occurring via alternate routes including those consistent with the energetic landscape model.