Conserved quantitative stability/flexibility relationships (QSFR) in an orthologous RNase H pair

Protein Function Analysis through Machine Learning

Machine learning (ML) has been an important arsenal in computational biology used to elucidate protein function for decades. With the recent burgeoning of novel ML methods and applications, new ML approaches have been incorporated into many areas of computational biology dealing with protein function. We examine how ML has been integrated into a wide range of computational models to improve prediction accuracy and gain a better understanding of protein function. The applications discussed are protein structure prediction, protein engineering using sequence modifications to achieve stability and druggability characteristics, molecular docking in terms of protein-ligand binding, including allosteric effects, protein-protein interactions and protein-centric drug discovery. To quantify the mechanisms underlying protein function, a holistic approach that takes structure, flexibility, stability, and dynamics into account is required, as these aspects become inseparable through their interdependence. Another key component of protein function is conformational dynamics, which often manifest as protein kinetics. Computational methods that use ML to generate representative conformational ensembles and quantify differences in conformational ensembles important for function are included in this review. Future opportunities are highlighted for each of these topics.

Comparisons of Ribonuclease HI Homologs and Mutants Uncover a Multistate Model for Substrate Recognition

Ribonuclease HI (RNHI) nonspecifically cleaves the RNA strand of RNA:DNA hybrid duplexes in a myriad of biological processes. Several RNHI homologs contain an extended domain, termed the handle region, which is critical to substrate binding. Nuclear magnetic resonance (NMR) spectroscopy and molecular dynamics (MD) simulations have suggested a kinetic model in which the handle region can exist in open (substrate-binding competent) or closed (substrate-binding incompetent) states in homologs containing arginine or lysine at position 88 (using sequence numbering of E. coli RNHI), while the handle region populates states intermediate between the open and closed conformers in homologs with asparagine at residue 88 [Stafford, K. A., et al., PLoS Comput. Biol. 2013, 9, 1-10]. NMR parameters characterizing handle region dynamics are highly correlated with enzymatic activity for RNHI homologs with two-state (open/closed) handle regions [Martin, J. A., et al., Biochemistry 2020, 59, 3201-3205]. The work presented herein shows that homologs containing asparagine 88 display distinct structural features compared with their counterparts containing arginine or lysine 88. Comparisons of RNHI homologs and site-directed mutants with asparagine 88 support a kinetic model for handle region dynamics that includes 12 unique transitions between eight conformations. Overall, these findings present an example of the structure-function relationships of enzymes and spotlight the use of NMR spectroscopy and MD simulations in uncovering fine details of conformational preferences.

Promiscuous Esterases Counterintuitively Are Less Flexible than Specific Ones

Understanding mechanisms of promiscuity is increasingly important from a fundamental and application point of view. As to enzyme structural dynamics, more promiscuous enzymes generally have been recognized to also be more flexible. However, examples for the opposite received much less attention. Here, we exploit comprehensive experimental information on the substrate promiscuity of 147 esterases tested against 96 esters together with computationally efficient rigidity analyses to understand the molecular origin of the observed promiscuity range. Unexpectedly, our data reveal that promiscuous esterases are significantly less flexible than specific ones, are significantly more thermostable, and have a significantly increased specific activity. These results may be reconciled with a model according to which structural flexibility in the case of specific esterases serves for conformational proofreading. Our results signify that an esterase sequence space can be screened by rigidity analyses for promiscuous esterases as starting points for further exploration in biotechnology and synthetic chemistry.

Can constraint network analysis guide the identification phase of KnowVolution? A case study on improved thermostability of an endo-β-glucanase

Cellulases are industrially important enzymes, e.g., in the production of bioethanol, in pulp and paper industry, feedstock, and textile. Thermostability is often a prerequisite for high process stability and improving thermostability without affecting specific activities at lower temperatures is challenging and often time-consuming. Protein engineering strategies that combine experimental and computational are emerging in order to reduce experimental screening efforts and speed up enzyme engineering campaigns. Constraint Network Analysis (CNA) is a promising computational method that identifies beneficial positions in enzymes to improve thermostability. In this study, we compare CNA and directed evolution in the identification of beneficial positions in order to evaluate the potential of CNA in protein engineering campaigns (e.g., in the identification phase of KnowVolution). We engineered the industrially relevant endoglucanase EGLII from Penicillium verruculosum towards increased thermostability. From the CNA approach, six variants were obtained with an up to 2-fold improvement in thermostability. The overall experimental burden was reduced to 40% utilizing the CNA method in comparison to directed evolution. On a variant level, the success rate was similar for both strategies, with 0.27% and 0.18% improved variants in the epPCR and CNA-guided library, respectively. In essence, CNA is an effective method for identification of positions that improve thermostability.

Quantifying the Relationship between Conformational Dynamics and Enzymatic Activity in Ribonuclease HI Homologues

Ribonuclease HI (RNHI), a ubiquitous, non-sequence-specific endonuclease, cleaves the RNA strand in RNA/DNA hybrids. RNHI functions in replication and genome maintenance, and retroviral reverse transcriptases contain an essential ribonuclease H domain. Nuclear magnetic resonance (NMR) spectroscopy combined with molecular dynamics (MD) simulations suggests a model in which the extended handle region domain of Escherichia coli RNHI populates (substrate-binding-competent) "open" and (substrate-binding-incompetent) "closed" states, while the thermophilic Thermus thermophilus RNHI mainly populates the closed state at 300 K [Stafford, K. A., Robustelli, P., and Palmer, A. G., III (2013) PLoS Comput. Biol. 9, 1-10]. In addition, an in silico-designed mutant E. coli Val98Ala RNHI was predicted to populate primarily the closed state. The work presented here validates this model and confirms the predicted properties of the designed mutant. MD simulations suggest that the conformational preferences of the handle region correlate with the conformations of Trp85, Thr92, and Val101. NMR residual dipolar coupling constants, three-bond scalar coupling constants, and chemical shifts experimentally define the conformational states of these residues and hence of the handle domain. These NMR parameters correlate with the Michaelis constants for RNHI homologues, confirming the important role of the handle region in the modulation of substrate recognition and illustrating the power of NMR spectroscopy in dissecting the conformational preferences underlying enzyme function.

Rigidity and flexibility characteristics of DD[E/D]-transposases Mos1 and Sleeping Beauty

DD[E/D]-transposases catalyze the multistep reaction of cut-and-paste DNA transposition. Structurally, several DD[E/D]-transposases have been characterized, revealing a multi-domain structure with the catalytic domain possessing the RNase H-like structural motif that brings three catalytic residues (D, D, and E or D) into close proximity for the catalysis. However, the dynamic behavior of DD[E/D]-transposases during transposition remains poorly understood. Here, we analyze the rigidity and flexibility characteristics of two representative DD[E/D]-transposases Mos1 and Sleeping Beauty (SB) using the minimal distance constraint model (mDCM). We find that the catalytic domain of both transposases is globally rigid, with the notable exception of the clamp loop being flexible in the DNA-unbound form. Within this globally rigid structure, the central β-sheet of the RNase H-like motif is much less rigid in comparison to its surrounding α-helices, forming a cage-like structure. The comparison of the original SB transposase to its hyperactive version SB100X reveals the region where the change in flexibility/rigidity correlates with increased activity. This region is found to be within the RNase H-like structural motif and comprise the loop leading from beta-strand B3 to helix H1, helices H1 and H2, which are located on the same side of the central beta-sheet, and the loop between helix H3 and beta-strand B5. We further identify the RKEN214-217DAVQ mutations of the set of hyperactive mutations within the catalytic domain of SB transposase to be the driving factor that induces change in residue-pair rigidity correlations within SB transposase. Given that a signature RNase H-like structural motif is found in DD[E/D]-transposases and, more broadly, in a large superfamily of polynucleotidyl transferases, our results are relevant to these proteins as well.

Ensemble- and Rigidity Theory-Based Perturbation Approach To Analyze Dynamic Allostery

Allostery describes the functional coupling between sites in biomolecules. Recently, the role of changes in protein dynamics for allosteric communication has been highlighted. A quantitative and predictive description of allostery is fundamental for understanding biological processes. Here, we integrate an ensemble-based perturbation approach with the analysis of biomolecular rigidity and flexibility to construct a model of dynamic allostery. Our model, by definition, excludes the possibility of conformational changes, evaluates static, not dynamic, properties of molecular systems, and describes allosteric effects due to ligand binding in terms of a novel free-energy measure. We validated our model on three distinct biomolecular systems: eglin c, protein tyrosine phosphatase 1B, and the lymphocyte function-associated antigen 1 domain. In all cases, it successfully identified key residues for signal transmission in very good agreement with the experiment. It correctly and quantitatively discriminated between positively or negatively cooperative effects for one of the systems. Our model should be a promising tool for the rational discovery of novel allosteric drugs.

Mutations in Antibody Fragments Modulate Allosteric Response Via Hydrogen-Bond Network Fluctuations

A mechanical perturbation method that locally restricts conformational entropy along the protein backbone is used to identify putative allosteric sites in a series of antibody fragments. The method is based on a distance constraint model that integrates mechanical and thermodynamic viewpoints of protein structure wherein mechanical clamps that mimic substrate or cosolute binding are introduced. Across a set of six single chain-Fv fragments of the anti-lymphotoxin-β receptor antibody, statistically significant responses are obtained by averaging over 10 representative structures sampled from a molecular dynamics simulation. As expected, the introduced clamps locally rigidify the protein, but long-ranged increases in both rigidity and flexibility are also frequently observed. Expanding our analysis to every molecular dynamics frame demonstrates that the allosteric responses are modulated by fluctuations within the hydrogen-bond network where the native ensemble is comprised of conformations that both are, and are not, affected by the perturbation in question. Population shifts induced by the mutations alter the allosteric response by adjusting which hydrogen-bond networks are the most probable. These effects are compared using response maps that track changes across each single chain-Fv fragment, thus providing valuable insight into how sensitive allosteric mechanisms are to mutations.

Novel Ricin Subunit Antigens With Enhanced Capacity to Elicit Toxin-Neutralizing Antibody Responses in Mice

RiVax is a candidate ricin toxin subunit vaccine antigen that has proven to be safe in human phase I clinical trials. In this study, we introduced double and triple cavity-filling point mutations into the RiVax antigen with the expectation that stability-enhancing modifications would have a beneficial effect on overall immunogenicity of the recombinant proteins. We demonstrate that 2 RiVax triple mutant derivatives, RB (V81L/C171L/V204I) and RC (V81I/C171L/V204I), when adsorbed to aluminum salts adjuvant and tested in a mouse prime-boost-boost regimen were 5- to 10-fold more effective than RiVax at eliciting toxin-neutralizing serum IgG antibody titers. Increased toxin neutralizing antibody values and seroconversion rates were evident at different antigen dosages and within 7 days after the first booster. Quantitative stability/flexibility relationships analysis revealed that the RB and RC mutations affect rigidification of regions spanning residues 98-103, which constitutes a known immunodominant neutralizing B-cell epitope. A more detailed understanding of the immunogenic nature of RB and RC may provide insight into the fundamental relationship between local protein stability and antibody reactivity.

Dynamics and thermodynamic properties of CXCL7 chemokine

Chemokines form a family of signaling proteins mainly responsible for directing the traffic of leukocytes, where their biological activity can be modulated by their oligomerization state. We characterize the dynamics and thermodynamic stability of monomer and homodimer structures of CXCL7, one of the most abundant platelet chemokines, using experimental methods that include circular dichroism (CD) and nuclear magnetic resonance (NMR) spectroscopy, and computational methods that include the anisotropic network model (ANM), molecular dynamics (MD) simulations and the distance constraint model (DCM). A consistent picture emerges for the effects of dimerization and Cys5-Cys31 and Cys7-Cys47 disulfide bonds formation. The presence of disulfide bonds is not critical for maintaining structural stability in the monomer or dimer, but the monomer is destabilized more than the dimer upon removal of disulfide bonds. Disulfide bonds play a key role in shaping the characteristics of native state dynamics. The combined analysis shows that upon dimerization flexibly correlated motions are induced between the 30s and 50s loop within each monomer and across the dimer interface. Interestingly, the greatest gain in flexibility upon dimerization occurs when both disulfide bonds are present, and the homodimer is least stable relative to its two monomers. These results suggest that the highly conserved disulfide bonds in chemokines facilitate a structural mechanism that is tuned to optimally distinguish functional characteristics between monomer and dimer.

Rigidity Emerges during Antibody Evolution in Three Distinct Antibody Systems: Evidence from QSFR Analysis of Fab Fragments

The effects of somatic mutations that transform polyspecific germline (GL) antibodies to affinity mature (AM) antibodies with monospecificity are compared among three GL-AM Fab pairs. In particular, changes in conformational flexibility are assessed using a Distance Constraint Model (DCM). We have previously established that the DCM can be robustly applied across a series of antibody fragments (VL to Fab), and subsequently, the DCM was combined with molecular dynamics (MD) simulations to similarly characterize five thermostabilizing scFv mutants. The DCM is an ensemble based statistical mechanical approach that accounts for enthalpy/entropy compensation due to network rigidity, which has been quite successful in elucidating conformational flexibility and Quantitative Stability/Flexibility Relationships (QSFR) in proteins. Applied to three disparate antibody systems changes in QSFR quantities indicate that the VH domain is typically rigidified, whereas the VL domain and CDR L2 loop become more flexible during affinity maturation. The increase in CDR H3 loop rigidity is consistent with other studies in the literature. The redistribution of conformational flexibility is largely controlled by nonspecific changes in the H-bond network, although certain Arg to Asp salt bridges create highly localized rigidity increases. Taken together, these results reveal an intricate flexibility/rigidity response that accompanies affinity maturation.

Antibodies are protective proteins used by the immune system to recognize and neutralize foreign objects through interactions with a specific part of the target, called an antigen. Antibody structures are Y-shaped, contain multiple protein chains, and include two antigen-binding sites. The binding sites are located at the end of the Fab fragments, which are the upward facing arms of the Y-structure. The Fab fragments maintain binding affinity by themselves, and are thus often used as surrogates to student antibody-antigen interactions. High affinity antibodies are produced during the course of an immune response by successive mutations to germline gene-encoded antibodies. Germline antibodies are more likely to be polyspecific, whereas the affinity maturation process yields monoclonal antibodies that bind specifically to the target antigen. In this work, we use a computational Distance Constraint Model to characterize how mechanical properties change as three disparate germline antibodies are converted to affinity mature. Our results reveal a rich set of mechanical responses throughout the Fab structure. Nevertheless, increased rigidity in the VH domain is consistently observed, which is consistent with the transition from polyspecificity to monospecificity. That is, flexible antibody structures are able to recognize multiple antigens, while increased affinity and specificity is achieved—in part—by structural rigidification.

Flexibility Correlation between Active Site Regions Is Conserved across Four AmpC β-Lactamase Enzymes

β-lactamases are bacterial enzymes that confer resistance to β-lactam antibiotics, such as penicillins and cephalosporins. There are four classes of β-lactamase enzymes, each with characteristic sequence and structure properties. Enzymes from class A are the most common and have been well characterized across the family; however, less is known about how physicochemical properties vary across the C and D families. In this report, we compare the dynamical properties of four AmpC (class C) β-lactamases using our distance constraint model (DCM). The DCM reliably predicts thermodynamic and mechanical properties in an integrated way. As a consequence, quantitative stability/flexibility relationships (QSFR) can be determined and compared across the whole family. The DCM calculates a large number of QSFR metrics. Perhaps the most useful is the flexibility index (FI), which quantifies flexibility along the enzyme backbone. As typically observed in other systems, FI is well conserved across the four AmpC enzymes. Cooperativity correlation (CC), which quantifies intramolecular couplings within structure, is rarely conserved across protein families; however, it is in AmpC. In particular, the bulk of each structure is composed of a large rigid cluster, punctuated by three flexibly correlated regions located at the active site. These regions include several catalytic residues and the Ω-loop. This evolutionary conservation combined with active their site location strongly suggests that these coupled dynamical modes are important for proper functioning of the enzyme.

Structures of mesophilic and extremophilic citrate synthases reveal rigidity and flexibility for function

Citrate synthase (CS) catalyses the entry of carbon into the citric acid cycle and is highly-conserved structurally across the tree of life. Crystal structures of dimeric CSs are known in both "open" and "closed" forms, which differ by a substantial domain motion that closes the substrate-binding clefts. We explore both the static rigidity and the dynamic flexibility of CS structures from mesophilic and extremophilic organisms from all three evolutionary domains. The computational expense of this wide-ranging exploration is kept to a minimum by the use of rigidity analysis and rapid all-atom simulations of flexible motion, combining geometric simulation and elastic network modeling. CS structures from thermophiles display increased structural rigidity compared with the mesophilic enzyme. A CS structure from a psychrophile, stabilized by strong ionic interactions, appears to display likewise increased rigidity in conventional rigidity analysis; however, a novel modified analysis, taking into account the weakening of the hydrophobic effect at low temperatures, shows a more appropriate decreased rigidity. These rigidity variations do not, however, affect the character of the flexible dynamics, which are well conserved across all the structures studied. Simulation trajectories not only duplicate the crystallographically observed symmetric open-to-closed transitions, but also identify motions describing a previously unidentified antisymmetric functional motion. This antisymmetric motion would not be directly observed in crystallography but is revealed as an intrinsic property of the CS structure by modeling of flexible motion. This suggests that the functional motion closing the binding clefts in CS may be independent rather than symmetric and cooperative.

Redistribution of flexibility in stabilizing antibody fragment mutants follows Le Châtelier's principle

Le Châtelier's principle is the cornerstone of our understanding of chemical equilibria. When a system at equilibrium undergoes a change in concentration or thermodynamic state (i.e., temperature, pressure, etc.), La Châtelier's principle states that an equilibrium shift will occur to offset the perturbation and a new equilibrium is established. We demonstrate that the effects of stabilizing mutations on the rigidity ⇔ flexibility equilibrium within the native state ensemble manifest themselves through enthalpy-entropy compensation as the protein structure adjusts to restore the global balance between the two. Specifically, we characterize the effects of mutation to single chain fragments of the anti-lymphotoxin-β receptor antibody using a computational Distance Constraint Model. Statistically significant changes in the distribution of both rigidity and flexibility within the molecular structure is typically observed, where the local perturbations often lead to distal shifts in flexibility and rigidity profiles. Nevertheless, the net gain or loss in flexibility of individual mutants can be skewed. Despite all mutants being exclusively stabilizing in this dataset, increased flexibility is slightly more common than increased rigidity. Mechanistically the redistribution of flexibility is largely controlled by changes in the H-bond network. For example, a stabilizing mutation can induce an increase in rigidity locally due to the formation of new H-bonds, and simultaneously break H-bonds elsewhere leading to increased flexibility distant from the mutation site via Le Châtelier. Increased flexibility within the VH β4/β5 loop is a noteworthy illustration of this long-range effect.

A case study comparing quantitative stability-flexibility relationships across five metallo-β-lactamases highlighting differences within NDM-1

The Distance Constraint Model (DCM) is an ensemble-based biophysical model that integrates thermodynamic and mechanical viewpoints of protein structure. The DCM outputs a large number of structural characterizations that collectively allow for Quantified Stability-Flexibility Relationships (QSFR) to be identified and compared across protein families. Using five metallo-β-lactamases (MBLs) as a representative set, we demonstrate how QSFR properties are both conserved and varied across protein families. Similar to our characterizations on other protein families, the backbone flexibility of the five MBLs are overall visually conserved, yet there are interesting specific quantitative differences. For example, the plasmid-encoded NDM-1 enzyme, which leads to a fast spreading drug-resistant version of Klebsiella pneumoniae, has several regions of significantly increased rigidity relative to the other four. In addition, the set of intramolecular couplings within NDM-1 are also atypical. While long-range couplings frequently vary significantly across protein families, NDM-1 is distinct because it has limited correlated flexibility, which is isolated within the active site S3/S4 and S11/H6 loops. These loops are flexibly correlated in the other members, suggesting it is important to function, but the others also have significant amounts of correlated flexibility throughout the rest of their structures.

Efficient and robust analysis of biomacromolecular flexibility using ensembles of network topologies based on fuzzy noncovalent constraints

We describe an approach (ENT(FNC)) for performing rigidity analyses of biomacromolecules on ensembles of network topologies (ENT) generated from a single input structure. The ENT is based on fuzzy noncovalent constraints, which considers thermal fluctuations of biomacromolecules without actually sampling conformations. Definitions for fuzzy noncovalent constraints were derived from persistency data from molecular dynamics (MD) simulations. A very good agreement between local flexibility and rigidity characteristics from ENT(FNC) and MD simulations-generated ensembles is found. Regarding global characteristics, convincing results were obtained when relative thermostabilities of citrate synthase and lipase A structures were computed. The ENT(FNC) approach significantly improves the robustness of rigidity analyses, is highly efficient, and does not require a protein-specific parameterization. Its low computational demand makes it especially valuable for the analysis of large data sets, e.g., for data-driven protein engineering.

Variations within class-A β-lactamase physiochemical properties reflect evolutionary and environmental patterns, but not antibiotic specificity

The bacterial enzyme β-lactamase hydrolyzes the β-lactam ring of penicillin and chemically related antibiotics, rendering them ineffective. Due to rampant antibiotic overuse, the enzyme is evolving new resistance activities at an alarming rate. Related, the enzyme's global physiochemical properties exhibit various amounts of conservation and variability across the family. To that end, we characterize the extent of property conservation within twelve different class-A β-lactamases, and conclusively establish that the systematic variations therein parallel their evolutionary history. Large and systematic differences within electrostatic potential maps and pairwise residue-to-residue couplings are observed across the protein, which robustly reflect phylogenetic outgroups. Other properties are more conserved (such as residue pKa values, electrostatic networks, and backbone flexibility), yet they also have systematic variations that parallel the phylogeny in a statistically significant way. Similarly, the above properties also parallel the environmental condition of the bacteria they are from in a statistically significant way. However, it is interesting and surprising that the only one of the global properties (protein charge) parallels the functional specificity patterns; meaning antibiotic resistance activities are not significantly constraining the global physiochemical properties. Rather, extended spectrum activities can emerge from the background of nearly any set of electrostatic and dynamic properties.

Protein mechanics: how force regulates molecular function

BACKGROUND

Regulation of proteins is ubiquitous and vital for any organism. Protein activity can be altered chemically, by covalent modifications or non-covalent binding of co-factors. Mechanical forces are emerging as an additional way of regulating proteins, by inducing a conformational change or by partial unfolding.

SCOPE

We review some advances in experimental and theoretical techniques to study protein allostery driven by mechanical forces, as opposed to the more conventional ligand driven allostery. In this respect, we discuss recent single molecule pulling experiments as they have substantially augmented our view on the protein allostery by mechanical signals in recent years. Finally, we present a computational analysis technique, Force Distribution Analysis, that we developed to reveal allosteric pathways in proteins.

MAJOR CONCLUSIONS

Any kind of external perturbation, being it ligand binding or mechanical stretching, can be viewed as an external force acting on the macromolecule, rendering force-based experimental or computational techniques, a very general approach to the mechanics involved in protein allostery.

GENERAL SIGNIFICANCE

This unifying view might aid to decipher how complex allosteric protein machineries are regulated on the single molecular level.

Constraint Network Analysis (CNA): a Python software package for efficiently linking biomacromolecular structure, flexibility, (thermo-)stability, and function

For deriving maximal advantage from information on biomacromolecular flexibility and rigidity, results from rigidity analyses must be linked to biologically relevant characteristics of a structure. Here, we describe the Python-based software package Constraint Network Analysis (CNA) developed for this task. CNA functions as a front- and backend to the graph-based rigidity analysis software FIRST. CNA goes beyond the mere identification of flexible and rigid regions in a biomacromolecule in that it (I) provides a refined modeling of thermal unfolding simulations that also considers the temperature-dependence of hydrophobic tethers, (II) allows performing rigidity analyses on ensembles of network topologies, either generated from structural ensembles or by using the concept of fuzzy noncovalent constraints, and (III) computes a set of global and local indices for quantifying biomacromolecular stability. This leads to more robust results from rigidity analyses and extends the application domain of rigidity analyses in that phase transition points ("melting points") and unfolding nuclei ("structural weak spots") are determined automatically. Furthermore, CNA robustly handles small-molecule ligands in general. Such advancements are important for applying rigidity analysis to data-driven protein engineering and for estimating the influence of ligand molecules on biomacromolecular stability. CNA maintains the efficiency of FIRST such that the analysis of a single protein structure takes a few seconds for systems of several hundred residues on a single core. These features make CNA an interesting tool for linking biomacromolecular structure, flexibility, (thermo-)stability, and function. CNA is available from http://cpclab.uni-duesseldorf.de/software for nonprofit organizations.

Global and local indices for characterizing biomolecular flexibility and rigidity

Understanding flexibility and rigidity characteristics of biomolecules is a prerequisite for understanding biomolecular structural stability and function. Computational methods have been implemented that directly characterize biomolecular flexibility and rigidity by constraint network analysis. For deriving maximal advantage from these analyses, their results need to be linked to biologically relevant characteristics of a structure. Such links are provided by global and local measures ("indices") of biomolecular flexibility and rigidity. To date, more than 14 indices are available with sometimes overlapping or only vague definitions. We present concise definitions of these indices, analyze the relation between, and the scope and limitations of them, and compare their informative value. For this, we probe the structural stability of the calcium binding protein α-lactalbumin as a showcase, both in the "ground state" and after perturbing the system by changing the network topology. In addition, we introduce three indices for the first time that extend the application domain of flexibility and rigidity analyses. The results allow us to provide guidelines for future studies suggesting which of these indices could best be used for analyzing, understanding, and quantifying structural features that are important for biomolecular stability and function. Finally, we make suggestions for proper index notations in future studies to prevent the misinterpretation and to facilitate the comparison of results obtained from flexibility and rigidity analyses.

Changes in Lysozyme Flexibility upon Mutation Are Frequent, Large and Long-Ranged

We investigate changes in human c-type lysozyme flexibility upon mutation via a Distance Constraint Model, which gives a statistical mechanical treatment of network rigidity. Specifically, two dynamical metrics are tracked. Changes in flexibility index quantify differences within backbone flexibility, whereas changes in the cooperativity correlation quantify differences within pairwise mechanical couplings. Regardless of metric, the same general conclusions are drawn. That is, small structural perturbations introduced by single point mutations have a frequent and pronounced affect on lysozyme flexibility that can extend over long distances. Specifically, an appreciable change occurs in backbone flexibility for 48% of the residues, and a change in cooperativity occurs in 42% of residue pairs. The average distance from mutation to a site with a change in flexibility is 17–20 Å. Interestingly, the frequency and scale of the changes within single point mutant structures are generally larger than those observed in the hen egg white lysozyme (HEWL) ortholog, which shares 61% sequence identity with human lysozyme. For example, point mutations often lead to substantial flexibility increases within the β-subdomain, which is consistent with experimental results indicating that it is the nucleation site for amyloid formation. However, β-subdomain flexibility within the human and HEWL orthologs is more similar despite the lowered sequence identity. These results suggest compensating mutations in HEWL reestablish desired properties.

The functional importance of protein dynamics is universally accepted, making the study of dynamical similarities and differences among proteins of the same function an intriguing problem. While some metrics are likely to be conserved across family, differences are also very common. In previous works we have used a Distance Constraint Model to quantify flexibility differences across sets of orthologous proteins, which reproduce this diversity. In the same manner, this work investigates changes occurring upon individual point mutations. Somewhat surprisingly, the small structural perturbations caused by mutation lead to changes throughout the protein. These changes can be quite large, actually surpassing the scale for differences between ortholog pairs. Moreover, changes in flexibility frequently occur at sites far from the mutation site. These results underscore the sensitivity of protein dynamics in connection with allostery, and help explain why differences across protein families are so common.

Calculating ensemble averaged descriptions of protein rigidity without sampling

Previous works have demonstrated that protein rigidity is related to thermodynamic stability, especially under conditions that favor formation of native structure. Mechanical network rigidity properties of a single conformation are efficiently calculated using the integer body-bar Pebble Game (PG) algorithm. However, thermodynamic properties require averaging over many samples from the ensemble of accessible conformations to accurately account for fluctuations in network topology. We have developed a mean field Virtual Pebble Game (VPG) that represents the ensemble of networks by a single effective network. That is, all possible number of distance constraints (or bars) that can form between a pair of rigid bodies is replaced by the average number. The resulting effective network is viewed as having weighted edges, where the weight of an edge quantifies its capacity to absorb degrees of freedom. The VPG is interpreted as a flow problem on this effective network, which eliminates the need to sample. Across a nonredundant dataset of 272 protein structures, we apply the VPG to proteins for the first time. Our results show numerically and visually that the rigidity characterizations of the VPG accurately reflect the ensemble averaged properties. This result positions the VPG as an efficient alternative to understand the mechanical role that chemical interactions play in maintaining protein stability.

Flexibility analysis of biomacromolecules with application to computer-aided drug design

Flexibility characteristics of biomacromolecules can be efficiently determined down to the atomic level by a graph-theoretical technique as implemented in the FIRST (Floppy Inclusion and Rigid Substructure Topology) and ProFlex software packages. The method has been successfully applied to a series of protein and nucleic acid structures. Here, we describe practical guidelines for setting up and performing a flexibility analysis, discuss current bottlenecks of the approach, and provide sample applications as to how this technique can support computer-aided drug design approaches.

Ensemble properties of network rigidity reveal allosteric mechanisms

The distance constraint model (DCM) is a unique computational modeling paradigm that integrates mechanical and thermodynamic descriptions of macromolecular structure. That is, network rigidity calculations are used to account for nonadditivity within entropy components, thus restoring the utility of free-energy decomposition. The DCM outputs a large number of structural characterizations that collectively allow for quantified stability-flexibility relationships (QSFR) to be identified. In this review, we describe the theoretical underpinnings of the DCM and introduce several common QSFR metrics. Application of the DCM across protein families highlights the sensitivity within the set of protein structure residue-to-residue couplings. Further, we have developed a perturbation method to identify putative allosteric sites, where large changes in QSFR upon rigidification (mimicking ligand-binding) detect sites likely to invoke allosteric changes.

How do thermophilic proteins and proteomes withstand high temperature?

We attempt to understand the origin of enhanced stability in thermophilic proteins by analyzing thermodynamic data for 116 proteins, the largest data set achieved to date. We compute changes in entropy and enthalpy at the convergence temperature where different driving forces are maximally decoupled, in contrast to the majority of previous studies that were performed at the melting temperature. We find, on average, that the gain in enthalpy upon folding is lower in thermophiles than in mesophiles, whereas the loss in entropy upon folding is higher in mesophiles than in thermophiles. This implies that entropic stabilization may be responsible for the high melting temperature, and hints at residual structure or compactness of the denatured state in thermophiles. We find a similar trend by analyzing a homologous set of proteins classified based only on the optimum growth temperature of the organisms from which they were extracted. We find that the folding free energy at the temperature of maximal stability is significantly more favorable in thermophiles than in mesophiles, whereas the maximal stability temperature itself is similar between these two classes. Furthermore, we extend the thermodynamic analysis to model the entire proteome. The results explain the high optimal growth temperature in thermophilic organisms and are in excellent quantitative agreement with full thermal growth rate data obtained in a dozen thermophilic and mesophilic organisms.

Nonadditivity in conformational entropy upon molecular rigidification reveals a universal mechanism affecting folding cooperativity

Previously, we employed a Maxwell counting distance constraint model (McDCM) to describe α-helix formation in polypeptides. Unlike classical helix-coil transition theories, the folding mechanism derives from nonadditivity in conformational entropy caused by rigidification of molecular structure as intramolecular cross-linking interactions form along the backbone. For example, when a hydrogen bond forms within a flexible region, both energy and conformational entropy decrease. However, no conformational entropy is lost when the region is already rigid because atomic motions are not constrained further. Unlike classical zipper models, the same mechanism also describes a coil-to-β-hairpin transition. Special topological features of the helix and hairpin structures allow the McDCM to be solved exactly. Taking full advantage of the fact that Maxwell constraint counting is a mean field approximation applied to the distribution of cross-linking interactions, we present an exact transfer matrix method that does not require any special topological feature. Upon application of the model to proteins, cooperativity within the folding transition is yet again appropriately described. Notwithstanding other contributing factors such as the hydrophobic effect, this simple model identifies a universal mechanism for cooperativity within polypeptide and protein-folding transitions, and it elucidates scaling laws describing hydrogen-bond patterns observed in secondary structure. In particular, the native state should have roughly twice as many constraints as there are degrees of freedom in the coil state to ensure high fidelity in two-state folding cooperativity, which is empirically observed.

Nonadditivity in the alpha-helix to coil transition

The Lifson-Roig Model (LRM) and all its variants describe the α-helix to coil transition in terms of additive component-free energies within a free energy decomposition scheme, and these contributions are interpreted through sequence-context dependent nucleation and propagation parameters. Although this phenomenological approach is able to adequately fit experimental data on helix content and heat capacity, the number of required parameters increases dramatically with additional sequence variation. Moreover, due to nonadditive competing microscopic effects that are difficult to disentangle within a LRM, large uncertainties within the parameters emerge. We offer an alternative view that removes the need for sequence-context parameterization by focusing on individual microsopic interactions within a free energy decomposition and explicitly account for nonadditivity in conformational entropy through network rigidity using a Distance Constraint Model (DCM). We apply a LRM and a DCM to previously published experimental heat capacity and helix content data for a series of heterogeneous polypeptides. Both models describe the experimental data well, and the parameters from both models are consistent with prior work. However, the number of DCM parameters is independent of sequence-variability, the parameter values exhibit better transferability, and the helix nucleation is predicted by the DCM explicitly through the nonadditive nature of conformational entropy. The importance of these results is that the DCM offers a system-independent approach for modeling stability within polypeptides and proteins, where the demonstrated accuracy for the α-helix to coil transition over a series of heterogeneous polypeptides described here is one case in point.

Allosteric response is both conserved and variable across three CheY orthologs

A computational method to identify residues likely to initiate allosteric signals has been developed. The method is based on differences within stability and flexibility profiles between wild-type and perturbed structures as computed by a distance constraint model. Application of the approach to three bacterial chemotaxis protein Y (CheY) orthologs provides a comparison of allosteric response across protein family divergence. Interestingly, we observe a rich mixture of both conservation and variability within the identified allosteric sites. While similarity within the overall response parallels the evolutionary relationships, >50% of the best scoring putative sites are only identified in a single ortholog. These results suggest that detailed descriptions of intraprotein communication are substantially more variable than structure and function, yet do maintain some evolutionary relationships. Finally, structural clusters of large response identify four allosteric hotspots, including the β4/α4 loop known to be critical to relaying the CheY phosphorylation signal.

Ensemble-based methods for describing protein dynamics

Molecular dynamics (MD) simulation is a natural approach for studying protein dynamics, and coupled with the ideas of multiscale modeling, MD proves to be the gold standard in computational biology to investigate mechanistic details related to protein function. In principle, if MD trajectories are long enough, the ensemble of protein conformations generated allows thermodynamic and kinetic properties to be predicted. We know from experiments that proteins exhibit a high degree of fidelity in function, and that empirical kinetic models are successful in describing kinetics, suggesting that the ensemble of conformations cluster into well-defined thermodynamic states, which are frequently metastable. The experimental evidence suggest that more efficient computational models that retain only essential properties of the protein can be constructed to faithfully reproduce the relatively few observed thermodynamic states, and perhaps describe transition states if the model is sufficiently detailed. Indeed, there are many so-called ensemble-based methods that attempt to generate more complete ensembles than MD can provide by focusing on the most important driving forces through simplified representations of how elements within the protein interact. Although coarse-graining is employed in MD and other approaches, such as in elastic network models, the key distinguishing factor of ensemble-based methods is that they are meant to efficiently generate a large ensemble of conformations without solving explicit equations of motion. This review highlights three types of ensemble-based methods, illustrated by 'COREX' and the Wako-Saito-Munoz-Eaton (WSME) model, the Framework Rigidity Optimized Dynamic Algorithm (FRODA) and the distance constraint model (DCM).

Helix/coil nucleation: a local response to global demands

A complete description of protein structure and function must include a proper treatment of mechanisms that lead to cooperativity. The helix/coil transition serves as a simple example of a cooperative folding process, commonly described by a nucleation-propagation mechanism. The prevalent view is that coil structure must first form a short segment of helix in a localized region despite paying a free energy cost (nucleation). Afterward, helical structure propagates outward from the nucleation site. Both processes entail enthalpy-entropy compensation that derives from the loss in conformational entropy on helix formation with concomitant gain in favorable interactions. Nucleation-propagation models inherently assume that cooperativity arises from a sequential series of local events. An alternative distance constraint model asserts there is a direct link between available degrees of freedom and cooperativity through the nonadditivity in conformational entropy. That is, helix nucleation is a concerted manifestation of rigidity propagating through atomic structure. The link between network rigidity and nonadditivity of conformational entropy is shown in this study by solving the distance constraint model using a simple global constraint counting approximation. Cooperativity arises from competition between excess and deficiency in available degrees of freedom in the coil and helix states respectively.

Thermostability in rubredoxin and its relationship to mechanical rigidity

The source of increased stability in proteins from organisms that thrive in extreme thermal environments is not well understood. Previous experimental and theoretical studies have suggested many different features possibly responsible for such thermostability. Many of these thermostabilizing mechanisms can be accounted for in terms of structural rigidity. Thus a plausible hypothesis accounting for this remarkable stability in thermophilic enzymes states that these enzymes have enhanced conformational rigidity at temperatures below their native, functioning temperature. Experimental evidence exists to both support and contradict this supposition. We computationally investigate the relationship between thermostability and rigidity using rubredoxin as a case study. The mechanical rigidity is calculated using atomic models of homologous rubredoxin structures from the hyperthermophile Pyrococcus furiosus and mesophile Clostridium pasteurianum using the FIRST software. A global increase in structural rigidity (equivalently a decrease in flexibility) corresponds to an increase in thermostability. Locally, rigidity differences (between mesophilic and thermophilic structures) agree with differences in protection factors.

Unifying mechanical and thermodynamic descriptions across the thioredoxin protein family

We compare various predicted mechanical and thermodynamic properties of nine oxidized thioredoxins (TRX) using a Distance Constraint Model (DCM). The DCM is based on a nonadditive free energy decomposition scheme, where entropic contributions are determined from rigidity and flexibility of structure based on distance constraints. We perform averages over an ensemble of constraint topologies to calculate several thermodynamic and mechanical response functions that together yield quantitative stability/flexibility relationships (QSFR). Applied to the TRX protein family, QSFR metrics display a rich variety of similarities and differences. In particular, backbone flexibility is well conserved across the family, whereas cooperativity correlation describing mechanical and thermodynamic couplings between the residue pairs exhibit distinctive features that readily standout. The diversity in predicted QSFR metrics that describe cooperativity correlation between pairs of residues is largely explained by a global flexibility order parameter describing the amount of intrinsic flexibility within the protein. A free energy landscape is calculated as a function of the flexibility order parameter, and key values are determined where the native-state, transition-state, and unfolded-state are located. Another key value identifies a mechanical transition where the global nature of the protein changes from flexible to rigid. The key values of the flexibility order parameter help characterize how mechanical and thermodynamic response is linked. Variation in QSFR metrics and key characteristics of global flexibility are related to the native state X-ray crystal structure primarily through the hydrogen bond network. Furthermore, comparison of three TRX redox pairs reveals differences in thermodynamic response (i.e., relative melting point) and mechanical properties (i.e., backbone flexibility and cooperativity correlation) that are consistent with experimental data on thermal stabilities and NMR dynamical profiles. The results taken together demonstrate that small-scale structural variations are amplified into discernible global differences by propagating mechanical couplings through the H-bond network.

Constraint counting on RNA structures: linking flexibility and function

RNA structures are highly flexible biomolecules that can undergo dramatic conformational changes required to fulfill their diverse functional roles. Constraint counting on a topological network representation of an RNA structure can provide very efficiently detailed insights into the intrinsic flexibility characteristics of the biomolecule. In the network, vertices represent atoms and edges represent covalent and strong non-covalent bonds and angle constraints. Initially, the method has been successfully applied to identify rigid and flexible regions in proteins. Here, we present recent progress in extending the approach to RNA structures. As a case study, we analyze stability characteristics of the ribosomal exit tunnel and relate these findings to the tunnel's active role in co-translational processes.

Conformational Entropy of an Ideal Cross-Linking Polymer Chain

We present a novel analytical method to calculate conformational entropy of ideal cross-linking polymers from the configuration integral by employing a Mayer series expansion. Mayer-functions describing chemical bonds within the chain and for cross-links are sharply peaked over the temperature range of interest, and, are well approximated as statistically weighted Dirac delta-functions that enforce distance constraints. All geometrical deformations consistent with a set of distance constraints are integrated over. Exact results for a contiguous series of connected loops are employed to substantiate the validity of a previous phenomenological distance constraint model that describes protein thermodynamics successfully based on network rigidity.

Hydrogen bond networks determine emergent mechanical and thermodynamic properties across a protein family

Background

Gram-negative bacteria use periplasmic-binding proteins (bPBP) to transport nutrients through the periplasm. Despite immense diversity within the recognized substrates, all members of the family share a common fold that includes two domains that are separated by a conserved hinge. The hinge allows the protein to cycle between open (apo) and closed (ligated) conformations. Conformational changes within the proteins depend on a complex interplay of mechanical and thermodynamic response, which is manifested as an increase in thermal stability and decrease of flexibility upon ligand binding.

Results

We use a distance constraint model (DCM) to quantify the give and take between thermodynamic stability and mechanical flexibility across the bPBP family. Quantitative stability/flexibility relationships (QSFR) are readily evaluated because the DCM links mechanical and thermodynamic properties. We have previously demonstrated that QSFR is moderately conserved across a mesophilic/thermophilic RNase H pair, whereas the observed variance indicated that different enthalpy-entropy mechanisms allow similar mechanical response at their respective melting temperatures. Our predictions of heat capacity and free energy show marked diversity across the bPBP family. While backbone flexibility metrics are mostly conserved, cooperativity correlation (long-range couplings) also demonstrate considerable amount of variation. Upon ligand removal, heat capacity, melting point, and mechanical rigidity are, as expected, lowered. Nevertheless, significant differences are found in molecular cooperativity correlations that can be explained by the detailed nature of the hydrogen bond network.

Conclusion

Non-trivial mechanical and thermodynamic variation across the family is explained by differences within the underlying H-bond networks. The mechanism is simple; variation within the H-bond networks result in altered mechanical linkage properties that directly affect intrinsic flexibility. Moreover, varying numbers of H-bonds and their strengths control the likelihood for energetic fluctuations as H-bonds break and reform, thus directly affecting thermodynamic properties. Consequently, these results demonstrate how unexpected large differences, especially within cooperativity correlation, emerge from subtle differences within the underlying H-bond network. This inference is consistent with well-known results that show allosteric response within a family generally varies significantly. Identifying the hydrogen bond network as a critical determining factor for these large variances may lead to new methods that can predict such effects.

New insight into long-range nonadditivity within protein double-mutant cycles

Additivity principles in chemistry, biochemistry, and biophysics have been used extensively for decades. Nevertheless, it is well known that additivity frequently breaks down in complex biomacromolecules. Nonadditivity within protein double mutant free energy cycles of spatially close residue pairs is a generally well-understood phenomenon, whereas a robust description of nonadditivity extending over large distances remains to be developed. Here, we test the hypothesis that the mutational effects tend to be nonadditive if two structurally well-separated mutated residues belong to the same rigid cluster within the wild type protein, and additive if they are located within different clusters. We find the hypothesis to be statistically significant with P-values that range from 10(-5) to 10(-6). To the best of our knowledge, this result represents the first demonstration of a statistically significant preponderance for nonadditivity over long distances. These findings provide new insight into the origins of long-range nonadditivity in double mutant cycles, which complements the conventional wisdom that nonadditivity arises in double mutations involving contacting residues. Consequently, these results should have far-reaching implications for a proper understanding of protein stability, structure/function analyses, and protein design.

Designable structures are easy to unfold

We study the structural stability of models of proteins for which the selected folds are unusually stable to mutation, that is, designable. A two-dimensional hydrophobic-polar lattice model was used to determine designable folds and these folds were investigated through Langevin dynamics. We find that the phase diagram of these proteins depends on their designability. In particular, highly designable folds are found to be weaker, i.e., easier to unfold, than low designable ones. We expect this to be related to protein flexibility.

Elucidating quantitative stability/flexibility relationships within thioredoxin and its fragments using a distance constraint model

Numerous quantitative stability/flexibility relationships, within Escherichia coli thioredoxin (Trx) and its fragments are determined using a minimal distance constraint model (DCM). A one-dimensional free energy landscape as a function of global flexibility reveals Trx to fold in a low-barrier two-state process, with a voluminous transition state. Near the folding transition temperature, the native free energy basin is markedly skewed to allow partial unfolded forms. Under native conditions the skewed shape is lost, and the protein forms a compact structure with some flexibility. Predictions on ten Trx fragments are generally consistent with experimental observations that they are disordered, and that complementary fragments reconstitute. A hierarchical unfolding pathway is uncovered using an exhaustive computational procedure of breaking interfacial cross-linking hydrogen bonds that span over a series of fragment dissociations. The unfolding pathway leads to a stable core structure (residues 22-90), predicted to act as a kinetic trap. Direct connection between degree of rigidity within molecular structure and non-additivity of free energy is demonstrated using a thermodynamic cycle involving fragments and their hierarchical unfolding pathway. Additionally, the model provides insight about molecular cooperativity within Trx in its native state, and about intermediate states populating the folding/unfolding pathways. Native state cooperativity correlation plots highlight several flexibly correlated regions, giving insight into the catalytic mechanism that facilitates access to the active site disulfide bond. Residual native cooperativity correlations are present in the core substructure, suggesting that Trx can function when it is partly unfolded. This natively disordered kinetic trap, interpreted as a molten globule, has a wide temperature range of metastability, and it is identified as the "slow intermediate state" observed in kinetic experiments. These computational results are found to be in overall agreement with a large array of experimental data.