Information

Why are Ramachandran angles of first and the last amino acid not necessary to define the full 3D structure of a protein chain?


I have come across an online ppt slide of the bioinformatic algorithm where it is said that first and the last amino acid Ramachandran angle is not necessary to tell all its internal coordinates. Is it really not necessary for any simulations while predicting the 3D structure of proteins? For instance, if we have a protein chain of N amino acid sequence, there will be 2N dihedral angles. Of which we need to specify only the 2N-2 (apart from the first and the last amino acid dihedral angles) dihedral angles. Why the two angles are ignored? Let's assume all the bond angles and bond lengths are provided for the entire chain.

I have seen this statement in the following ppt, slide number -9, titled representation of internal coordinates, under bullet point 3.

https://www.cs.umb.edu/~nurith/cs612/Manipulation.pdf


Think about the case of a chain of N=3 points in space. There are only angles associated with the middle point, the two endpoints don't have angles associated with them because they only have one incident segment rather than two.


The NKCC and NCC Genes

E. Predicted but not Demonstrated Topologies of SLC12A1, 2 and 3 Proteins

A protein topology predicted in silico is halfway from the peptide sequence to the real three-dimensional structure of the protein (von Heijne, 2006 ). Hence, computer algorithms developed to predict protein topology or structure based on physicochemical properties of amino acid sequences as well as by comparison with known protein structures (e.g. threading and homology modeling) are invaluable tools to infer topology and/or function–structure relationships.

Most of the SLC12A proteins appear to share similar predicted structures with several transmembrane domains and long intracellular N- or C-termini. This assumption is based on the estimated hydrophilicity/hydrophobicity profiles of deduced SLC12A protein sequences according to Kyte-Doolittle’s algorithm ( Kyte and Doolittle, 1982 ). A key feature of this algorithm is the so-called “window size”, i.e. the number of amino acids examined at a time to determine a point of hydrophobic character ( Kyte and Doolittle, 1982 ). Hence, it is critical to choose a window size that corresponds to the expected size of the structural motif under investigation (i.e. a window size of 19–21 (about the size of a membrane spanning α-helix) will make hydrophobic, membrane-spanning domains stand out on the Kyte-Doolittle scale (typically >1.6)). However, windows sizes ranging from 11 to 15 amino acids were used to generate hydropathy plots predicting 12 transmembrane (TM) domains in mammalian members of the SLC12A family ( Caron et al., 2000 Delpire et al., 1994 Gamba et al., 1994 Gillen et al., 1996 Hiki et al., 1999 Moore-Hoon and Turner, 1998 Payne and Forbush, 1994 Payne et al., 1996 Yerby et al., 1997 ). Although alternative topological models for members of the SLC12A family have been proposed ( Park and Saier, 1996 ) and several transport protein families include members that probably have more or less than 12 TM domains ( Espanol and Saier, 1995 Paulsen and Skurray, 1993 ), it is accepted that the SLC12A family are proteins of 12 TM domains.

It is now clear that the most important factor in determining membrane insertion is the hydrophobicity of 19–21 amino acid sequences ( Zhao and London, 2006 ). This concept is better represented by using the experimentally determined transfer-free energies (ΔG) for each amino acid (i.e. a thermodynamic scale of hydrophobicity) originally proposed by Wimley and White ( Wimley and White, 1996 ). Hence, the hydrophobicity plot of Wimley-White (also known as the octanol plot) identifies the position of transmembrane α-helices in protein sequences with less ambiguity than the Kyte-Doolittle plot. As shown in Fig. 11.3 , octanol plots obtained for SLC12A1 (NKCC2), SLC12A2 (NKCC1) and SLC12A3 (NCC) are different to the ones originally proposed for these gene products using the Kyte-Doolittle algorithm with a window size of 11–15 ( Delpire et al., 1994 Gamba et al., 1994 Payne and Forbush, 1994 Yerby et al., 1997 ). However, the octanol plot correlates very well with the Kyte-Doolittle plot if the latter is constructed using a window size of 19–21 amino acids ( Fig. 11.3 ).

Figure 11.3 . Kyte-Doolittle and White-Wimley plots of NKCC2 and NKCC1 protein sequences. A. Predicted NKCC protein topology. Putative transmembrane domains (TM) are indicated as gray boxes across the lipid bilayer. The position of NKCC2 amino acids predicted to be localized at the TM domains is numbered underneath each potential TM domain. The continuous gray line represents the amino acid chain of the NKCC2 proteins. Colored dots located at the cytoplasmic N-terminal and C-terminal portions of NKCC2s represent the location of residues predicted to be phosphorylated (blue: Ser, green: Thr and black: Tyr) and potential N-glycosylation sites (red dots). The potential site for tyrosine sulfination at the N-terminus of NKCC2 is indicated with an arrowhead. Phosphorylation and sulfination sites on NKCC2 proteins were predicted using NetPhos ( www.cbs.dtu.dk/services/NetPhos ) and Sulfinator ( www.expasy.ch/tools/sulfinator ), respectively. B. Hydropathy plots of hNKCC2A (ABU69043) (top), rNKCC2A (ABU63482) (center) and hNKCC1a (AAC50561) (bottom). These analyses were performed using a window size of 19 residues. Window sizes of 19 or 21 make hydrophobic, membrane-spanning domains stand out clearly (typically, values &gt 1.6 on the Kyte and Doolittle scale). Under these conditions, hNKCC2 proteins are predicted to have 10 TM regions: 174–198, 208–228, 259–279, 298–318, 323–349, 380–402, 413–441, 489–512, 551–579 and 604–627. Each TM is ∼20 residues in length and highly identical among species. All predicted TMs in NKCC2s have energetic preferences for being in the lipid environment as characterized by the total free energy (ΔG) above zero in the White-Wimley interface hydropathy plot. The mean charges of the amino acids are calculated by giving the residues D (Asp) and E (Glu) a charge of −1, K (Lys) and R (Arg) a charge of +1, and the residue H (His) a charge of +0.5. The represented data were obtained using jEMBOSS for Linux (emboss.sourceforge.net/Jemboss), TMap, TMPredProtScale (at the ExPASy molecular biology server) and PROTEUS Structure Prediction Server v2.0 ( wks16338.biology.ualberta.ca/proteus ).

Prediction algorithms based solely on hydrophobicity plots ( Kyte and Doolittle, 1982 ) or thermodynamic scales of hydrophobicity ( Wimley and White, 1996 ) are somewhat incomplete and innacurate. The fact that ∼5% of the transmembrane α-helices in the known structures are very short (<15 residues) and only partially span the membrane, together with the lack of critical thermodynamic data, have made transmembrane prediction algorithms somewhat unsatisfactory. It was not until recently that the free energy contributions from individual amino acids in different positions along the membrane was reported ( Hessa et al., 2007 ). Hence, the accuracy of algorithms predicting TM helixes has been recently improved by the development of new tools such as MemBrain ( Shen and Chou, 2008 ), TopPred ΔG ( Hessa et al., 2007 ), SCAMPI ( Bernsel et al., 2008 ), ZPRED ( Granseth et al., 2006 ) and PRO/PRODIV-TMHMM ( Viklund and Elofsson, 2004 ). Most of these algorithms are part of TOPCONS protein topology prediction server (topcons.cbr.su.se). By using MemBrain or SCAMPI, human SLC12A1, SLC12A2 and SLC12A3 proteins (i.e. NKCC2, NKCC1 and NCC) it can be predicted that these proteins may have 13 TM domains, whereas PRODIV, PRO or OCTUPUS predicts 12 TM domains ( Fig. 11.4 ). It should be mentioned that the model having 13 TM domains places the N- and C-termini in different compartments (inside and outside, respectively), which is not supported by current experimental evidence.

Figure 11.4 . Consensus prediction of membrane protein topology. The topological information of hNKCC2A, hNKCC1a and hNCCa proteins (GenBank ABU69043, AAC50561 and AAC50355, respectively) was generated by using five different algorithms: SCAMPI, OCTUPUS, ZPRED, PRO/PRODIV-TMHMM ( topcons.cbr.su.se/ ) and MemBrain, an algorithm used for predicting ends of TM domains that are shorter than 15 residues. A. Predicted topology of NKCC and NCC proteins according to the algorithms used (MemBrain (red), SCAMPI (blue), PRO/PRODIV and TOPCONS (green)). Predicted TM domains are indicated as gray boxes across the lipid bilayer. The NKCC/NCC amino acid location predicted to be located in each TM is numbered underneath each transmembrane domain and varies according to the algorithm used. The continuous gray line represents the amino acid chain of the NKCC/NCC proteins whereas dotted lines represent the potential topologies according to the algorithm used. The cytoplasmic N-terminal and C-terminal portions of NKCCs/NCC are indicated. B. Predicted total free energy (ΔG) values of each residue in hNKCC2A (top), hNKCC1a (center) and hNCCa (bottom) protein sequences.


Background

Protein taxonomy [1–5] reveals that crystallographic protein structures have surprisingly little conformational diversity. It might be that the majority of different conformations have already been found [6, 7]. This apparent convergence in protein structure provides the rationale for the development of comparative modelling or threading techniques [8–12]. These approaches try to predict the tertiary structure of a folded protein using libraries of known protein structures as templates. According to the community-wide Critical Assessment for Structural Prediction (CASP) tests [13], at the moment this kind of methods have the best predictive power to determine a folded conformation.

In the loop regions, comparative modelling approaches still continue lacking in their precision [14, 15]. It is not uncommon that there are gaps in the loop regions that need to be filled by various insertion techniques. The success in loop modelling is also often limited to super-secondary structures where α-helices and β-strands are connected to each other by relatively short twists and turns [16, 17]. In the case of a very short loop, with no more than three residues, the shape can be determined by a combination of geometrical considerations and stereochemical constraints [18]. In the case of longer loops, both template based and template independent methods are being developed to predict their shapes [19–21]. The underlying assumption is that the number of loop conformations which can be accommodated by a given sequence should be limited. The different fragments which are already available in the Protein Data Bank (PDB) [22] database could then be used like Lego bricks, as structural building blocks in constructing the loops. A given amino acid sequence is simply divided into short fragments, and the shape of the ensuing loop is deduced using homologically related fragments that have known structures. The entire protein is then assembled by joining these fragments together. For the process of joining the fragments, both all-atom energy functions and comparisons with closely homologous template structures in the Protein Data Bank can be utilised [8, 9, 12, 14].

In the present article we propose a new systematic, purely quantitative method to identify and classify the modular building blocks of PDB loops we identify a loop following the DSSP [23] convention. Our approach is based on a first-principles energy function [24–29]. It is built on the concept of universality [30–36] to model the fragments of even long protein loops in terms of different parameterisations of a unique kink that solves a variant [37, 38] of the discrete nonlinear Schrödinger (DNLS) equation [39, 40]. Our starting point is the observation made in [41] that over 92 % of loops in those PDB structures that have been measured with better than 2.0 Å resolution, can be composed from 200 different parameterisations of the kink profile, with better than 0.65 Ångström RMSD (root-mean-square-distance) accuracy. Here we refine this observation, with the aim to develop it into a systematic loop fragment classification scheme. For this we consider only those ultrahigh precision PDB structures that have been measured with better than 1.0 Å resolution. This ensures that the B-factors in the loop regions are small, and in particular that the structures have not been subjected to extensive refinement procedures. Indeed, two loop fragments should be considered different only, when the average interatomic distance is larger than the average Debye-Waller B-factor fluctuation distance. If the B-factors are large, any systematic attempt to identify and/or distinguish two fragments becomes ambiguous. In the case of these intra-high resolution structures we can aim for the RMSD precision of 0.2 Å. We estimate this to be the extent of zero point fluctuations i.e. a distance around 0.2 Å corresponds to the intrinsic uncertainty in the determination of heavy atom positions along the protein backbone. Thus any difference less than 0.2 Å between average atomic coordinates is essentially undetectable. By explicit constructions, we show how in the case of this subset of ultrahigh resolution PDB protein structures, the loops can be systematically modeled using combinations of the unique kink of the generalised DNLS equation. As such, our approach provides a foundation for a general approach to classify loops in high precision crystallographic PDB structures, in terms of an energy function based first-principles mathematical concept.


Results

Global better than local model for residue couplings

Mutual information does not sufficiently correlate with residue proximity.

We first attempted the prediction of residue-residue proximity relationships using the straightforward local mutual information (MI) measure. MI(i,j) for each residue pair i, j is a difference entropy which compares the experimentally observed co-occurrence frequencies fij(Ai,Aj) of amino-acid pairs Ai, Aj in positions i, j of the alignment to the distribution fi(Ai)fj(Aj) that has no residue pair couplings (details in Text S1): (1) Contact maps constructed from residue pairs assigned high MI values, and thus interpreted as predicted contacts, differ substantially from the correct contact maps deduced from native structures, consistent with the work of Fodor et al. [9] (Figure S1). Visual inspection of MI-predicted contacts as lines connecting residue pairs superimposed on the observed crystal structure confirms that the contacts predicted from MI are often incorrect and/or unevenly distributed (Figure 3, left, blue lines). Presumably this arises due to the local nature of MI, which is independently calculated for each residue pair i,j. Plausibly, the key confounding factor is the transitivity of pair correlations, where the simplest case involves residue triplets for example, if residue B co-varies with both A and C, because B is spatially close to both A and C, then A and C may co-vary even without physical proximity (A–C is a transitive pair correlation). Any local measure of correlation, not just mutual information, is limited by this transitivity effect.

Extraction of evolutionary information about residue coupling and predicted contacts from multiple sequence alignments works much better using the global statistical model (right, Direct Information, DI, Equation 3) than the local statistical model (left, Mutual Information, MI, Equation 1). Predicted contacts for DI (red lines connecting the residues predicted to be coupled from sequence information) are better positioned in the experimentally observed structure (grey ribbon diagram), than those for MI (left, blue lines), shown here for the RAS protein (upper) and ELAV4 protein (lower). The DI residue pairs are also more evenly distributed along the chain and overlap more accurately with the contacts in the observed structure (red stars [predicted, grey circles [observed] in contact map center, upper right triangle) than those using MI (blue [predicted], grey circles [observed] center, lower left triangle). Details of contact maps for all proteins comparing predicted and observed contacts are in Figures S1 and S2, Text S1.

Effective residue couplings from a global maximum entropy model.

To disentangle such direct and indirect correlation effects, we use a global statistical model to compute a set of direct residue couplings that best explains all pair correlations observed in the multiple sequence alignment (see Methods and Text S1) [15], [47]. More precisely, we seek a general model, P(A1…AL), for the probability of a particular amino acid sequence A1…AL of length L to be a member of the iso-structural family under consideration, such that the implied probabilities Pij(Ai,Aj) for pair occurrences (marginals) are consistent with the data. In other words, we require Pij(Ai,Aj)∼fij(Ai,Aj), where fij(Ai,Aj) are the observed pair frequencies of amino acids at positions i and j in the known sequences in the family and the marginals Pij(Ai,Aj) are calculated by summing P(A1…AL) over all amino acid types at all sequence positions other than i and j. As specification of residue pair properties (ignoring higher order terms) leaves the amino acid sequence underdetermined, there are many probability models that would be consistent with the observed pair frequencies. One can therefore impose an additional condition, the maximum entropy condition, which requires a maximally even distribution of the probabilities - while still requiring consistency with data. Probability distributions that are solutions of this constrained optimization problem are of the form [11], [45], [49]: (2) Here Ai and Aj are particular amino acids at sequence positions i and j, and Z is the normalization constant. The Lagrange multipliers eij(Ai,Aj) and hi(Ai) constrain the agreement of the probability model with pair and single residue occurrences, respectively. This global statistical model is analogous to statistical physics expressions for the probability of the configuration of a multiple particle system, such as in the Ising or Potts models. In this analogy, a sequence position i corresponds to a particle, such as a spin, and can be in one of 21 states (Ai = 1..21) and, the Hamiltonian (the expression in curly brackets) consists of a sum of particle-particle coupling energies eij(Ai,Aj) and single particle coupling energies to external fields hi(Ai).

For our protein sequence problem, the eij(Ai,Aj) in equation 2 are essential residue couplings that are used in the prediction of folding constraints and the hi(Ai) are single residue terms that reflect consistency with observed single residue frequencies. These parameters are thus optimal with respect to the two key conditions, (1) consistency with observed data (pair and single residue frequencies) and (2) maximum entropy of the global probability over the set of all possible sequences. In practice, once these parameters are determined by matrix inversion (Equations M4, M5), one can directly compute the effective pair probabilities Pij Dir (Ai,Aj) (Equation M6), and from these the effective residue couplings (‘direct information’, in analogy to the term ‘mutual information’) DIij by summing over all possible amino acid pairs Ai,Aj at positions i,j: (3) The crucial difference between this expression for direct information DIij (Equation 3) and the equation for mutual information MIij (Equation 1) is to replace pair probabilities estimated based on local frequency counts fij(Ai,Aj), by the doubly constrained pair probabilities Pij Dir (Ai,Aj), which are globally consistent over all pairs i,j.

Global maximum entropy statistical model reveals residue proximity.

We now examine whether the residue coupling scores DIij (Equation 3 Equation 22, Text S1) from the maximum entropy model provide information about spatial proximity. Are residue pairs with higher DIij scores more likely to be close to each other in 3D structure? Examination of contact maps displaying residue pairs with highly ranked DIij values, overlaid onto contact maps for an observed (crystal) structure, reveals a surprisingly accurate match. The high-scoring residue pairs are often close in the observed structure, and these pairs are well distributed throughout the protein sequence and structure, in contrast to pairs with high-scoring MIij values, (Figure 3, Figure S2). This remarkable level of correct contact prediction holds for all of our test cases (Table 1, Table S1) in the four main fold classes.

Others have shown that given sufficient correct (true positive) contacts combined with a lack of incorrect (false positive) contacts, predicted contacts can be implemented as residue-residue distance restraints to fold proteins from the main four fold categories with up to ∼200 residues to under 3 Å Cα-RMSD error from the crystal structure [50] and, in later work, up to 365 residues with accuracy under 3 Å Cα-RMSD error [50], [51]. We were therefore encouraged to use our blindly predicted proximity relations as residue-residue distance restraints to fold proteins de novo from extended polypeptide chains.

Protein all-atom structures inferred from evolutionary constraints

In spite of elegant analyses using subsets of real contacts [50], [51], it is not a priori obvious to what extent accuracy of contact prediction translates to accuracy of 3D structure prediction and, in particular, how robust such prediction is to the presence of false positives. We therefore decided to assess the accuracy of contact prediction by the very stringent criterion of accuracy of predicted 3D structures.

Generating model structures.

Starting from an extended polypeptide chain with the amino acid sequence of a protein from the family (Table S1) we used well-established distance geometry algorithms, as used for structure determination by nuclear magnetic resonance (NMR) spectroscopy [52] (Text S1). The distance constraints were constructed using residue pairs with high DI scores pairs and secondary structure constraints predicted from sequence (Text S1, Appendix A1, Table S2). The protocol generates initial 3D conformations and then applies simulated annealing [48] (steps outlined in Text S1 and Appendix A2). We reasoned that the number of distance constraints (NC) needed should scale monotonically with the protein length L, as seen in fold reconstruction from observed contact maps [50], [51]. To explore the variability of predicted structure using a given set of distance restraints, we generated 20 candidate structures for a range of NC values which started at NC = 30 and incremented in steps of 10 to the nearest multiple of 10 to L, e.g., from NC = 30 to NC = 160 for the Hras proteins which has 160 core residues in the PFAM alignment. Thus, in total we generate on the order of 2*L candidate three-dimensional structures for each protein family as prediction candidates, more precisely, between 400 and 560, depending on the size of the protein (Table 1, Appendix A3). In practice, a smaller number of candidate structures may be sufficient. Each candidate is an all-atom structure prediction for a particular reference protein of interest chosen from the family. The model structures satisfy a maximal fraction of the predicted distance constraints and meet the conditions of good stereochemistry and consistency with non-bonded intermolecular potentials. The top predicted structure for each protein is selected by blind ranking of these candidate structures using objective, primarily geometric, criteria (Figure 2, Figure S2, Appendix A3).

3D structure inference for small and larger proteins of diverse fold types

To evaluate the information content of residue pair correlations with respect to protein fold prediction, we apply the method to increasingly difficult cases. We start with small single-domain proteins and move on to larger, more difficult targets, eventually covering a set of well-studied protein domains of wide-ranging biological interest, from different fold classes. We report detailed results for four example families, and summary results for 11 further test families, and provide detailed 3D views of all 15 test protein families in Figure S3 and detailed 3D coordinates and Pymol session files for interactive inspection in Appendices A3 and A4, http://cbio.mskcc.org/foldingproteins.

Small: an RNA binding domain (RRM).

The blind prediction of the 71-residue RRM domain of the human Elav4 protein (Uniprot ID: Elav4_human) is a typical example of a smaller protein. The distance constraints are derived from a rich corpus of 25K example proteins in the PFAM family. The highest ranking predicted structure has a (excellent) low 2.9 Å Cα -RMSD deviation from the crystal structure over 67 out of 71 residues, a TM score of 0.57 and GDT_TS 54.6, indicating overall good structural similarity to the observed crystal structure, [53], [54], (Figure 2 top, Table 1). It has correct topography of the five β-strands and two α-helices, marred only by a missing H-bond pattern between strands 1 and 3, at least partly due to the truncation of the strand 1, a consequence of the short length of the sequence in the PFAM alignment. Strands 2 and 3 align with only 1.6 Å Cα-RMSD deviation over the length of the predicted strands and are positioned well enough for hydrogen bonding, with some correct registration. Interestingly, the 4 th β-strand (penultimate) missed by the secondary structure prediction method is placed in the correct region in 3D: this is one of several examples in which residue coupling information overrides incorrect local prediction. The predicted top-ranked domain of Elav4 very likely lies within the refinement basin of the native structure.

Medium size: Ras oncogene (G-domain), an α/β domain with an GTPase active site.

The G-domain family in PFAM, with Human Ras proto-oncogene protein (Uniprot-ID: hras_human) chosen as the protein of interest, has a core multiple sequence alignment (MSA) of 161 residues. The structure has an α/β fold with a 6-stranded β-sheet, surrounded by 5 α-helices, one of which (α-2) is involved in the GTPase switch transition after GTP hydrolysis. The highest ranked, blindly predicted structure is 3.6 Å Cα-RMSD from the crystal structure, over 161 residues (Figure 2 middle) and has a high TM score of 0.7 (range 0.0–1.0, with 1.0 implying 100% of residues are within a set distance from the correct position [53]). The six β-strands and five α-helices are placed in the correct spatial positions and are correctly threaded (Appendices A3 and A4). The 6 β-strands, which make 5 β-strand pairs are not within hydrogen boding distance for all backbone bonding, but the correct register can be easily predicted for 26/30 of the residue pairs, Text S1. The accuracy of overall topography of the highest-ranked structures is remarkable (Table 1) and, as far as we know, currently not achievable for proteins of this size by any de novo structure prediction method [27].

Larger: trypsin, an enzyme with a two-domain β-barrel structure.

The largest (non-membrane) protein family tested in the blind test is the trypsin-fold serine protease family, with rat trypsin chosen as a representative protein. Its size, at 223 amino acids, is significantly larger than proteins that can be predicted by other de novo computational methods. Trypsin consists of β-strands in two structurally isomorphous β-barrel domains. The highest-ranked predicted structure has 4.3 Å Cα-RMSD error over 186 out of 223 residues (Figure 2 bottom, Table 1, Appendices A3 and A4). The overall distribution of secondary structure elements in space is approximately correct and our method correctly predicts 5 disulfide bonded cysteine pairs, which lie within our alignment, Text S1. The topography of the first β-barrel (domain 1) is good and plausibly within refinement range of the observed structure. Five correct pairs of β-strands are identified (one absent) and 70% of hydrogen bonding paired residues are predicted with correct register, Text S1. However, domain 2 has a number of incorrect loop progressions (see Pymol session in Appendix A3), and possibly (by inspection) is not within refinement range of the correct structure. Predicting the structure of proteins in the trypsin family is particularly challenging, as the structure is known to undergo a conformational change after cleavage of the activation peptide [55] and, as the N-terminal and C-terminal peptide cross from one domain to the other.

Inferring the residue configuration in the active site of trypsin.

In spite of the limited quality of structure prediction in domain 2 of trypsin, it is interesting that the top-ranked structures place the Cα atoms of the highly conserved active site triad residues Ser-His-Asp in correct relative spatial proximity, i.e., within 0.64 3 Å Cα-RMSD (and 1.3 Å all atom-RMSD) error, after superimposition of the three residues of the catalytic site with the same three residues of the experimental structure (Figure S4). This may reflect strong evolutionary constraints near functional sites and may imply that the configuration of resides around an active site can be predicted more accurately than other detailed aspects of the 3D structure. The ability to predict active site constellations at this level of accuracy would be particularly interesting for the design of drugs on predicted structural templates.

Exploration: rhodopsin, an α-helical transmembrane protein.

Rhodopsin is the first membrane protein predicted using this method. This important class of membrane proteins has 7 helices and the PFAM family from which the distance restraints are inferred contains many subfamilies of class A G-protein coupled receptors [56]. For the highest ranked predicted rhodopsin structure (4.84 Å Cα-RMSD error from a representative crystal structure over 171 residues), the overall topography of the helices is accurate (TM score 0.5), with most of the positional deviation arising for helices 1 and 7, which are misaligned relative to the direction perpendicular to the membrane surface, (Table 1, Figure S3). The predicted structure with the highest TM score (0.55), and 4.29 Å Cα-RMSD over 180 residues, also misaligns the terminal helices but does recapitulate a network of close distances (<4.5 Å) between the side chains of Arg135 (helix III) and Glu247, Thr251 (helix VI) as well as other well-known inter-helical proximities such as Asn78 (helix II) to Trp161 (helix IV) and Ser127 (helix III) [57]. Given that the current version of the method has no information about membrane orientation for membrane proteins, this constitutes an excellent starting point for future application of the method to 3D structure prediction for membrane proteins.

Ranking inferred structures.

To arrive at useful and objective blind predictions, the set of inferred structures for each family is ranked by objective criteria based on physical principles and a priori knowledge of general principles of protein structure. In the current implementation, we use consistency with the well-established empirical observation of right-handed chain twist in α-helices and right-handed inter-strand twist for β-strand pairs [58] (Text S1). The virtual dihedrals of the α-helices and the predicted β-twists in the candidate structures were combined together as a score, weighted by the relative numbers of residues in β-strands and α-helices for each protein, see scores for all structures in Appendix A5. We found these geometric criteria effective in eliminating artifacts that appear to arise from the fact that distance constraints do not have any chiral information, such that the starting structures prior to refinement using molecular dynamics, while consistent with distance constraints, may have incorrect chirality, either globally or locally. We also eliminated candidate structures with knots (as with the top ranked trypsin prediction) according to the method of Mirny et al. [59].

The highest-ranked all-atom model structure is taken as the top blindly predicted structure (Table 1, Table S1). Lower ranked structures are expected to have lower accuracy of 3D structure, but this has to be tested after blind prediction by comparison with known structures. As a test of the entire procedure and the ranking criteria, we assessed our blind predictions by comparing the ranking score of the predicted structures with the experimentally observed structure, from X-ray crystallography, of the chosen reference protein, (Text S1, Figure 4A, Figure S5 and Appendix A5). For proteins such as RAS and Trypsin (Figure 4B), the objective criteria successfully ranks those predicted structures with the lowest Cα-RMSD error to a crystal structure as highest scoring. As we remove obviously knotted proteins [59] we would miss genuinely knotted proteins [60] which are, however, rarely observed.

A. The overall performance of the de novo structure prediction reported here based on contacts inferred from evolutionary information (EICs), ranges from good to excellent for the 15 test proteins (on left: 3D structure type [α = α-helix-containing, β = β-strand-containing, 7tm-α = containing seven trans-membrane helices] in parentheses: size of protein domain/number of residues used for Cα-RMSD error calculation on bar: Uniprot database ID). Larger bars mean better performance, i.e., lower Cα-RMSD co-ordinate error. Left: performance for the top ranked structure for each target protein out of 400–560 (depending on the size of the protein, 20 structures per NC bin, NC in steps of 10, details in Appendix A3 and A6) candidate structures in blind prediction mode right: performance of the best structure, in hindsight, out of 20 candidate structures generated, for 20 sets of constraints ranging from 10∶200, in steps of 10. This reflects what would be achievable with better ranking criteria or independent post-prediction validation of structure quality (Table 1 details of blind ranking scores in Web Appendix A5). Other well-accepted methods for error assessment, such as GDT-TS and TM score are useful for comparison purposes (Table S1, Web Appendix A6). B. Ranking score of each candidate structure (quantifying expected structure quality) versus Cα-RMSD error. Ideally, higher-ranking scores correspond to lower error. The distribution of the candidate structures (black dots) for Elav4, Ras and Trypsin shows, in retrospect, that the ranking criteria used here are relatively useful and help in anticipating which structures are likely to be best (plots for all tested proteins in Figure S5). In blind prediction mode, a list of predicted candidate 3D structure has to be ranked by objective and automated criteria, with a single top ranked structure or a set of top ranked structures nominated as preferred predictions.

Assessment of prediction accuracy: 3D structures

Summary of blinded 3D accuracy for 15 test proteins of known structure.

We were surprised at the extent and high value of the information in the derived distance constraints about the 3D fold of examples from all major fold classes containing various proportions of α-helices and β-sheets. This high information content in residue couplings, derived from the maximum entropy statistical model, extends, so far, to proteins as large as G-domains, like H-ras, with 161 residues, and serine proteases, like trypsin, with 223 residues, as well as the rhodopsin family, a trans-membrane protein, with 258 aligned residues. This size has so far been out of range for state-of-the-art de novo prediction methods even when three-dimensional fragments are used [22], [61]. In general we find that predicted α/β folds, among the 15 proteins investigated in detail, produce the most accurate overall topography (Table 1, Table S1, Figure S5.). We anticipate that these results will likely extend to many protein families and that accurate structures can be generated for many of these using distance constraints derived from evolutionary information and predicted secondary structure alone, followed by energy refinement. For 12 out of the set of 15 protein families (Table 1), the top blindly ranked structures have coordinate errors from 2.7 Å–4.8 Å for at least 75% of the residues, using the accepted practice of omitting a moderate fraction of badly fitting residues in order to avoid exaggerated influence from outliers resulting from the square in the definition of Cα-RMSD (using the MaxCluster suite [62]). For most practical purposes, one might consider these to be within the basin of attraction within which one is highly likely to be able to identify the particular correct fold, which we estimate roughly to have a radius of about 5 Å Cα-RMSD. The partial exceptions are rhodopsin (OPSD) for which the relatively low 4.8 Å Cα-RMSD error is limited to 171 out of 258 residues (66%) and PCBP1 at 4.7 Å for 46/63 residues (73%). For these proteins, the agreement is limited to a smaller, though still sizable, fraction of the protein and it is less likely that the correct overall fold would be recognized. The major exception is SPTB2 at 4.0 Å for 47/108 residues (44%), which we consider not satisfactory. The TM scores customary in CASP reflect these differences and it is plausible that the top-ranked predictions for 11 out of the 15 test proteins would be considered excellent for de novo modeled structures of this size (Table S1) [27], [61], [63].

Detailed examination of the close contacts of top ranked predicted structures reveals interesting violations, (Figure 5). For Ras and Trypsin false positive DI constraints (between Ser145 and Asp57 for Ras, and Ser127 and Ala37 for trypsin) are not satisfied in the top predicted structures thereby improving the accuracy. Conversely, a contact is made the N-terminal β-strand and the C-terminal helix in RAS and C-terminal β-strand in ELAV4, despite the fact that no constraints are used in the vicinity of these contacts (grey circles, Figure 5).

The top blindly ranked structures are evaluated in terms of quality of contact prediction (NC = 40 for Elav4, NC = 130 for Ras, NC = 160 for Trypsin). The predicted constraints (red stars) are correct when they coincide with contacts derived from the observed structure (grey circles) and otherwise incorrect (false positives, red on white). The contacts derived from the predicted 3D structure (dark blue) are in good general agreement with those from the observed structure (grey). The cooperative nature of the folding prediction process permits favorable situations, in which contacts regions not touched by a predicted constraint (red) are still predicted correctly (black circle for RAS, dark blue on grey, no red) and false positive constraints are not strong enough to lead to incorrect contacts (left black circle Elav4, red star, no dark blue or grey). However, in unfavorable situations missing constraints may imply that contact regions are fully or partially missed (black circle, trypsin) or mostly missed (right black circle for Elav4, grey adjacent to and wider than dark blue).

Best 3D prediction accuracy in top 400 candidate structures.

To assess the potential of the method and with a view toward future improvements of ranking criteria for sets of candidate structures, one can ask the question, from hindsight, which of, say, 400 candidate structure has the highest accuracy. This question is analogous to protein structure prediction reports that discuss the relationship (scatter plots) of, e.g., model energy against model error. Here, the best candidate structures by TM score, selected from among 400 candidate structures for each protein (NC = 10–200), have TM scores from 0.5 to 0.76 and typically a lower error than the blindly top ranked structure, ranging from 2.8 Å to 4.6 Å Cα-RMSD for all 15 families, covering at least 80% of the residues, with the exception of OPSD where we achieve 4.3 Å for 180/258 residues (66%), (Figure 4B, Table1, Table S1). The fact that in most cases better 3D structures are found in the top 400 candidates is a non-trivial positive indication, as the conformational search space of protein folds is so large, that random methods, or moderately effective methods, would have an exceedingly low probability of achieving errors in this low range in as few as 400 structures. However, some of the structures generated here among the top 400 appear topologically incorrect, with the polypeptide chain passing through loops in a way that is, according to visual intuition, atypical of fully correct structures. Such topologically incorrectly structures would not be within a basin of attraction of conventional energy refinement, e.g., by simulated annealing. This indicates that neither low Cα-RMSD as a measure of overall accuracy, nor the more recently developed template modeling (TM) score, nor the global distance test - total score (GDT-TS), is fully informative indicators of structure quality. These classic structure comparison metrics need to be supplemented by more sophisticated measures, which quantify topographical differences in chain progression in 3D space, a direction for future work [64], [65], together with an analysis of violations of constraints in the spirit of Miller et al. [3]. In any case, the encouragingly high accuracy of the folds we generate amongst a relatively small number of candidates imply that improved ranking criteria may lead to a better set of top-ranked, fully blinded predictions.

Current technical limits of 3D prediction accuracy.

As an estimate of the accuracy maximally achievable by this method and its particular implementation, we performed reference calculations using artificial, fully correct, distance constraints derived from the experimentally observed structure. With this ideal set of constraints, we can construct protein structure models at an error of not lower than about 2.0 Å Cα-RMSD (Text S1, Table S3, larger values for some of the larger proteins). This places a lower bound on the expected error, inherent in the distance geometry and refinement part of the method and this error will scale to some extent with the length of the protein as others have noted [50]. That we achieve candidate structures close to these bounds with predicted distance constraints is consistent with the notion that the inferred residue couplings contain almost all the information required to find the native protein structure, at least for the 15 protein families examined here. This technical lower limit also represents a challenge for generic methods improvement for computation of accurate all-atom structures from distance constraints.

Assessment of prediction accuracy

Accuracy of contact prediction.

The accuracy of prediction of 3D structures crucially depends on the accuracy of contact prediction and the choice of distance constraints from a set of predicted contacts. Note that residue-residue proximity is a different requirement than residue-residue contact, as residues may be near each other in space without any of their atoms, being in inter-atomic contact (defined as inter-atomic distance near the minimum of non-bonded inter-atomic potentials (‘van der Waals’), say, about 3.5 Å). Here, we use the term inter-residue contact interchangeably with inter-residue proximity, i.e. minimum atom distance of less than 5 Angstroms. We assess the accuracy of contact prediction in terms of the number of true positives and false positives among predicted contacts, i.e., those that agree and those that disagree with the contacts observed in known 3D protein structures.

We find that the highest scoring pairs provide remarkably accurate information about residue-residue proximity (Figure 6A, Figures S6 and S7). For example, the rate of true positives is above 0.8 for the first 50 pairs for HRAS and still above 0.5 for the first 200 pairs for other proteins, it is lower but still relatively high, e.g., above 0.7 and 0.4 for the first 50 and 200 for ELAV4. These results are consistent with our parallel evaluation of contact prediction accuracy for a large number of bacterial protein domains [47] and represent a significant improvement over local methods of contact prediction from correlated mutations or co-evolution. Not surprisingly, there is a general trend for a higher rate of true positive contact prediction to results in better predicted 3D structures, The predicted structures of proteins such as Ras and CheY with a high proportion of true positive predicted contacts tend to be more accurate than those with lower rates, for example the KH domain of PCBP1 and the calponin homology domain of SPTB2. However, this relationship between the proportion of true positives and the accuracy of the best-predicted structures is not as simple as one might have expected, Figures S6, S8 and S9. For instance the thioredoxin predicted structures are on the whole more accurate than the predicted the lectin domain (A8MVQ9_HUMAN) structures despite the fact that thioredoxin has a lower true positive rate than lectin domain for its predicted contacts. Since the quality of 3D structures could depend also on the distribution of the contacts through the chain, for each protein we also calculated the distance of a experimental contact to the nearest predicted contact and this ‘spread’ showed a good correlation with the Cα-RMSD accuracy achieved, (Figure S10 and Text S1).

Evaluation of accuracy in terms of predicted contacts (A) and predicted 3D structures (B). (A) The two global models, the Bayesian network model (BNM, green [13]) and direct information model (DI, red, this work and [47]) have a consistently high rate of correctly predicted contacts (true positives) among the top NC ranked residue pairs two local models, mutual information (MI, green, equation 1) and SCA (black, [66]) have a consistently lower rate of true positives. Here, local refers to statistical independence of each pair i,j, while global refers to statistical consistency of all pairs. In (B), only the predicted 3D structures (green, BNM red, EIC) for the global models agree well with the observed structure (grey) Cα-RMSDs is calculated over the number or residues in parentheses (Pymol sessions for all structures in Web Appendix A4). Attempts to generate 3D structures for the two local methods MI and SCA failed (not shown). Comparing (A) and (B) confirms that a higher rate of true positives for contact prediction leads to better 3D structures and that for DI one needs at least a true positive rate of about 0.5 for about 100 predicted contacts, depending on size and other details of particular protein families. Interestingly, a false positive rate as high as about 0.3–0.5 can still be consistent with good 3D structure prediction.

Comparison of contact prediction accuracy between global and local models.

How well do other contact prediction methods work? The two global models, the Bayesian Network Model (BNM, [13], [46]) and the DI model (this work and [15] have a consistently high rate of correctly predicted contacts (true positive rate) among the top NC ranked residue pairs in comparison two local models, MI (Equation 1) and statistical coupling analysis (SCA, [66]), both have a lower rate of true positives (Figure 6A, Figures S6, S7, S11, S12, S13, S14, and S15). The relatively high accuracy of contact prediction in the BNM model encouraged us to generate predicted 3D structures based on the BNM ranked residue pairs as the basis for inferred distance constraints, following the protocol developed for the DI model. For ten test proteins, folded all-atom 3D structures for BNM agree well with the observed structure (green structures in Figure 6B and data not shown). On the whole, the Cα-RMSD errors are somewhat higher for the structures from the BNM model than those for the DI model (red structures in Figure 6B). In particular, using the notation [protein identifier/error for BNM/error for DI], we have: [RASH/5.6 Å/2.8 Å], [ELAV4/3.8 Å/2.6 Å], [YES/4.6 Å/3.6 Å] [CADH/4.7 Å/3.9 Å] and trypsin did not reach an accuracy lower than 12 Å Cα-RMSD with the BNM constraints (Figure 6B and data not shown). On the other hand, the BNM and the DI predictions for OMPR were in the same accuracy range when compared to the experimental structure, as the BNM result was over 74 atoms as opposed to 63 atoms for the DI method [OMPR/4.4 Å/4.0 Å].

These results confirm that in general a higher rate of true positives for contact prediction leads to better 3D structure prediction and, that for the global methods one needs at least a true positive rate of about 0.5 and on the order of about 100 predicted contacts, depending on size and other details of particular protein families. Interestingly, a false positive rate as high as about 0.3–0.5 can still be consistent with good 3D structure prediction. Clearly, the global statistical models provide a substantial increase in the accuracy of prediction of residue contacts and of 3D structures.

Information requirements for improved prediction of 3D structures

Requirement of sufficient sequence range coverage by the multiple sequence alignment.

Among the test set of twelve protein families, the lowest accuracy was obtained for the SPBT2 and rhodopsin proteins, (see Table 1, Table S1, Figure S3). In these cases a significant number of key residues are not included in the PFAM hidden Markov model (HMM) and thus were excluded from our analysis. If the alignment covers only part of the structure, the statistical model of the sequence is restricted to this part of the structure and does not provide information for non-covered regions. Since regions not covered by the PFAM alignments are often at the N-terminus or C-terminus of the protein and these are in contact in many protein structures, this will significantly harm the accuracy of prediction that is possible. Our analysis also shows that prediction is less likely to be accurate even within the covered region when ends of the alignment are absent. How much additional sequence information is required to build an alignment for the entire protein sequence in each case? This question is non-trivial as the diversity sampled at each sequence position by evolution varies greatly. Indeed the strength of structural evolutionary constraints may diminish towards the protein termini, analogous to the ‘frayed ends’ observed in many NMR-determined structures.

Correct folding with a surprisingly small number of distance constraints.

What is the minimum number of predicted distance constraints needed to generate an approximate 3D fold? An important parameter of our folding protocol is the number of inferred distance constraints, NC, used to generate candidate structures. While residues with the highest ranked pair correlations are usually close in 3D structure (Figures S6 and S7) the reliability decreases with decreasing value of DIij. We assessed the accuracy of the predicted protein folds for 15 evaluation families as a function of NC (Figures 7A and S16, Table S1).

A. How many distance constraints are needed for fold prediction? What fraction of false positives can be tolerated? With increasing number of predicted essential distance constraints (NC, horizontal axis), 3D prediction error decreases rapidly, as assessed by Cα-RMSD between the best of 20 (in each NC bin) predicted structures and the observed structure (here, for the 15 test proteins, using Pymol). Remarkably, as few as ∼NRES/2 (∼L/2) distance constraints dij (with chain distance |i−j|>5) suffice for good quality predictions below 5 Å Cα-RMSD, where NRES is the number of amino acid residues in the protein multiple sequence alignment. We therefore routinely generated candidate protein structures for up to NC = NRES distance constraints for blinded ranking (and for up to NC = 200 for other tests). Eventually the number of false positives does degrade prediction quality, e.g., for the 58 residue protein BPTI once NC is about 80 (1.5 NRES) the prediction quality is lost. In practice, we do not recommend using NC>NRES, i.e, more than about one constraint dij with |i−j|>5, per residue. B. When would it have been possible to fold from sequence? The increase in the number of sequences available in public databases (here, from successive archival releases of the PFAM collection of protein family alignments) is one of two key elements in the ability to predict protein folds from correlated mutations. Nevertheless plotting the numbers of sequences and dates shows that it would have been possible to calculate the structures up to 10 years ago for some proteins and that amazingly few sequences are sufficient. For example, although the retrospective prediction error (vertical axis, Cα-RMSD, using Pymol) for the best 3D structure (of 400 candidates each) in four protein families (Ras, SH3 domain (YES_human) and RnaseH from Ecoli) has decreased over time, the decrease is not strictly monotonic, as the result of non-systematic growth of the database. The point at which a predicted protein structure from a particular family reaches below 4 Å Cα-RMSD varies considerably. For example, while RnaseH required about 6000 sequence to dip below 4 Å error, reached around 2008, the structure of CheY could have been predicted to 3.3 Å Cα-RMSD, with only the 600 sequences available in 1999.

Going from 10 to typically 200 distance constraints, we find that the prediction error drops sharply as EIC constraints are added, until false positives gradually start to degrade the prediction quality. We conclude that one needs about 0.5 to 0.75 predicted constraints per residue, or about 25–35% of the total number of contacts, to achieve reasonable 3D structure prediction. This number is close to those reported by other groups, who used fully correct close residue pairs to impose inexact distances as constraints [50], [51], [67]. For instance, Elav4 (length 71) folds to below 5 Å Cα-RMSD with only 20 constraints, whilst Trypsin (length 223) takes 130 constraints. However, the number of constraints per residue to reach below 5 Å Cα-RMSD is not constant (column 15 Table S1), and proteins such as OMPR at 0.66 constraints per residue, and Ras at 0.25 constraints per residue show that this will depend on other factors, such as type of fold and false positive rates. While the accuracy of structure prediction for some proteins clearly decreases as the number of false positives, for example Cadh1, Elav4 and Yes, other proteins, such as Ras and CheY stay the same or even improve in accuracy as the false positive proportion increases, (Figure S8). This result underlines the necessity of using the constraints to attempt to fold the proteins, in order to assay the quality of predicted contacts, rather than relying on true positive rates of contact prediction alone.

Increasing prediction accuracy over time, but lower than expected numbers of sequences needed.

Since we not require today's standard of high performance computing, we wondered how long ago it would have been possible to make good structural predictions. How does the accuracy of predicted folds depend on the number of sequences in the multiple sequence alignment and their evolutionary diversity? To start to explore these questions we computed the accuracy of folding using distance constraints for four representative proteins, using alignments from 20 different releases of PFAM [1] covering the last 13 years. For each multiple sequence alignment we calculated 20 structures for a range of constraints from 30–200, (Figure 7B). During this period the available sequence information has increased dramatically as the result of new sequencing technology and large-scale genome projects, so we examined the best structure attained as a function of the number of sequences. Although there is a clear overall trend for the Cα-RMSD of predicted structures to drop monotonically as the number of sequences in the family increases (for example, RnaseH, 4 Å Cα-RMSD threshold was reached in 2009 when the number of sequences reached 5000), not all protein families behave the same way. The predicted Ras structures reached under a 4 Å Cα-RMSD in 2002 with as few as 1200 sequences, then, surprisingly, rose again as more sequences were included, to finally dip to 2.5 Å Cα-RMSD in 2009. Similarly, although the predicted structures of CheY and the SH3 domain from the Yes protein improve with the number of sequences available, predicted structures had Cα-RMSD in errors as low as 3.3 Å and 4.7 Å respectively in 1999, with ∼600 sequences for both. (Figure 7B). Most surprisingly, a predicted OMPR structure with an error under 5 Å Cα-RMSD would have been possibly using as few as 170 sequences (1999 PFAM release).

Hence our results highlight the overall relationship of accuracy of the predicted fold to the number of sequences available. However, this relationship is not straightforward. The distribution of sequences in the sequence space of a particular family will doubtless have an effect. In our current implementation of the algorithm, sequences with over 70% residue identity to family neighbors are down-weighted (Text S1). Therefore the effective number of sequences used for the DI coupling calculation is far less than the size of the family. Approximately only 12–40% of sequences available in the family are actually used for the calculation (Table S1). This reduction in the effective number of sequences varies substantially between families, highlighting the different distributions over sequence space covered by individual families (column 18 in Table S1). We speculate that future work will improve our understanding of which, as well as how many sequences are optimal for contact inference from evolutionary information.


Background

Protein structure validation methods like MolProbity [1] and Procheck [2] help crystallographers to find and fix potential problems that are incurred during fitting and refinement. These methods are commonly based on a priori chemical knowledge and utilise various well tested and broadly accepted stereochemical paradigms. Likewise, template based structure prediction and analysis packages [3] and molecular dynamics force fields [4] are customarily built on such paradigms. Among these, the Ramachandran map [5,6] has a central role. It is widely deployed both to various analyses of the protein structures, and as a tool in protein visualisation. The Ramachandran map describes the statistical distribution of the two dihedral angles φ and ψ that are adjacent to the Cα carbons along the protein backbone. A comparison between the observed values of the individual dihedrals in a given protein with the statistical distribution of the Ramachandran map is an appraised method to validate the backbone geometry.

In the case of side chain atoms, visual analysis methods like the Ramachandran map have been introduced. For example, the Janin map [7] can be used to compare observed side chain dihedrals such as χ1 and χ2 in a given protein, against their statistical distribution, in a manner which is analogous to the Ramachandran map.

Crystallographic refinement and validation programs like Phenix [8], Refmac [9] and others, often utilize the statistical data obtained from the Engh and Huber library [10,11]. This library is built using small molecular structures that have been determined with a very high resolution. At the level of entire proteins, side chain restraints are commonly derived from analysis of high resolution crystallographic structures [12,13] in Protein Data Bank (PDB) [14]. A backbone independent rotamer library [15] makes no reference to backbone conformation. But the possibility that the side-chain rotamer population depends on the local protein backbone conformation, was considered already by Chandrasekaran and Ramachandran [16]. Subsequently both secondary structure dependent [17], see also [7] and [15], and backbone dependent rotamer libraries [18,19] have been developed. We note that the subject remains under active investigation [20-25].

The information content in the secondary structure dependent libraries and the backbone independent libraries essentially coincide [13]. Both kinds of libraries are used extensively during crystallographic protein structure model building and refinement. But for the prediction of side-chain conformations, for example in the case of homology modeling and protein design, there can be an advantage to use the more revealing backbone dependent rotamer libraries.

In x-ray crystallographical protein structure experiments, the skeletonisation of the electron density map is a common technique to interpret the data and to build the initial model [26]. The Cα atoms are located at the branch points between the backbone and the side chain. As such they are subject to relatively stringent stereochemical constraints this is the reason why model building often starts with the initial identification of the skeletal Cα trace. The central role of the Cα atoms is widely exploited in structural classification schemes such as CATH [27] and SCOP [28], in various threading modeling techniques such as I-Tasser [29] and homology base approaches including SWISS-MODEL [30] and other related methods [31], in de novo approaches [32], and in the development of coarse grained energy functions for folding prediction [33]. As a consequence the so-called Cα-trace problem has become the subject of extensive investigations [34-38]. The resolution of the problem would consist of an accurate main chain and/or all-atom model of the folded protein, based on the knowledge of the positions of the central Cα atoms only. Both knowledge-based approaches such and MAXSPROUT [34] and de novo methods including PULCHRA [37] and REMO [38] have been developed, to try and resolve the Cα- trace problem. In the case of the backbone atoms, the geometric algorithm introduced by Purisima and Scheraga [39], or some variant thereof, is commonly utilized in these approaches. For the side chain atoms, most approaches to the Cα trace problem rely either on a statistical or on a conformer rotamer library in combination with steric constraints, complemented by an analysis which is based on diverse scoring functions. For the final fine-tuning of the model, all-atom molecular dynamics simulations can also be utilised.

In the present article we introduce and develop new generation visualisation techniques that we hope will become a beneficial complement to existing methods for protein structure analysis, refinement and validation. We use the Cα Frenet frames [40,41] to visualise the side chain. The output we aim at, is a 3D “what-you-see-is-what-you-have” type visual map of the statistically preferred all-atom model, calculable in terms of the Cα coordinates. As such, our approach should have value for example during the construction and validation of the initial backbone and all-atom models of a crystallographic protein structure.

Our approach is based on developments in three dimensional visualisation and virtual reality, that have taken place after the Ramachandran map was introduced. In lieu of the backbone dihedral angles that appear as coordinates in the Ramachandran map and correspond to a toroidal topology, we employ the geometry of virtual spheres that surround each heavy atom. We visually describe all the higher level heavy backbone and side chain atoms on the surface of a sphere, level-by-level along the backbone and side chains, exactly in the manner how they are seen by an imaginary, geometrically determined and Cα based miniature observer who roller-coasts along the backbone and climbs up the side chains, proceeding from one Cα atom to the next. At the location of each Cα our virtual observer orients herself consistently according to the purely geometrically determined Cα based discrete Frenet frames [40,41]. Thus the visualisation depends only on the Cα coordinates, and there is no reference to the other atoms in the initialisation of the construction. The other atoms - including subsequent Cα atoms along the backbone chain - are all mapped on the surface of a sphere that surrounds the observer, as if these atoms were stars in the sky.

At each Cα atom, the construction proceeds along the ensuing side chain, until the position of all heavy atoms have been determined. As such our maps provide a purely geometric and equitable, direct visual information on the statistically expected all- atom structure in a given protein.

The method we describe in this article, can form a basis for the future development of a novel approach to the Cα trace problem. As a complement to the existing approaches such as MAXSPROUT [34], PULCHRA [37] and REMO [38], the method we envision accounts for the secondary structure dependence in the heavy atom positions, which we here reveal. A secondary-structure dependent method to resolve the Cα trace problem should lead to an improved accuracy in the heavy atom positions, in terms of the Cα coordinates. In particular, since rotameric states do display clear secondary structure dependence, a fact that is sometimes overlooked in the development of rotamer libraries. The present article serves as a proof-of-concept.


Structural Biochemistry/Proteins

A protein is a functional biological molecule that is made up of one or more polypeptides that are folded/coiled into a specific structure [1] . Proteins are important macromolecules that serve as structural elements, transportation channels, signal receptors and transmitters, and enzymes. Proteins are linear polymer that are built up of the monomer units called amino acids. There are 20 different amino acids and they are connected by a peptide bond between the carboxyl group and the amino group in a linear chain called a polypeptide. Each protein has different side chains or the "R" groups. Proteins have many different active functional groups attached to them to help define their properties and functions. Proteins cover a wide range of functions, ranging from very rigid structural elements to transmitting information between cells. Each person has several hundred thousands of different proteins in their body. Proteins fold into secondary, tertiary, and quaternary structures based on intra-molecular bonding between functional groups or intermolecular bonding (quaternary only) and can obtain on a variety of three-dimensional shapes depending on the amino acid sequence. All proteins have primary, secondary and tertiary structures but quaternary structures only arise when a protein is made up of two or more polypeptide chains [1] . The folding of proteins is also driven and reinforced by the formation of many bonds between different parts of the chain. The formation of these bonds depends on the amino acid sequence. The study of their structures is important because proteins are essential for every activity in the human body as well as they are the key components of biological materials. Primary structure is when amino acids are linked together by peptide bonds to form polypeptide chains. Secondary structure is when the polypeptide chains fold into regular structures like the beta sheets, alpha helix, turns, or loops. A functional protein is much more than just a polypeptide, it is one or more polypeptides that have been precisely folded into a molecule with a very specific, unique shape which is critical to its function [1] .


Proteins are usually portrayed in 3D structures and categorized into four different characteristics and levels:

Primary: The primary structure of a protein is the level of protein structure which refers to the specific sequence of amino acids [1] . When two amino acids are in such a position that the carboxyl groups of each amino acid are adjacent to each other, they can be combined by undergoing a dehydration reaction which results in the formation of a peptide bond [1] . Amino acids in a polypeptide (protein) are linked by peptide bonds that begin with the N-terminal with a free amino group and ends at C-terminal with a free carboxyl group. rts . The peptide bond is planar and cannot rotate freely due to a partial double bond character. While there is a restricted rotation about peptide bond, there are two free rotations on (N-C) bond and (C-C) bond, which are called torsion angles, or more specifically the phi and psi angles. The freedoms of rotation of these two bonds are also limited due to steric hindrance. Genes carry the information to make polypeptides with a defined amino acid sequence. An average polypeptide is about 300 amino acids in length, and some genes encode polypeptides that are a few thousand amino acids long. It's important to know the primary structure of the protein because the primary structure encodes motifs that are of functional importance in their biological function structure and function are correlated at all levels of biological organization [1] .

Secondary: The amino acid sequence of a polypeptide, together with the laws of chemistry and physics, cause a polypeptide to fold into a more compact structure. Amino acids can rotate around bonds within a protein. This is the reason proteins are flexible and can fold into a variety of shapes. Folding can be irregular or certain regions can have a repeating folding pattern. The coils and folds that result from the hydrogen bonds between the repeating segments of the polypeptide backbone are called secondary structures [1] . Although the individual hydrogen bonds are weak, they are able to support a specific shape for that part of the protein due to the fact that they are repeated many times over a long part of the chain [1] . Secondary structures of a protein are proposed by Pauling and Corey. Its structures are formed by amino acids that are located within short distances of each other. Because of the planar nature of the peptide bonds, only certain types of secondary structure exist. The three important secondary structures are α-helix, β-sheets, and β-turns. Also, the beta sheets can be parallel, antiparallel, or mixed. Antiparallel beta sheets are more stable because the hydrogen bonds are at a ninety degree angles. The a-helix is a coiled structure stabilized by intrachain hydrogen bonds.

Characteristics of the Secondary Structures:

1. α-helix: In an α-helix, the polypeptide backbone forms a repeating helical structure that is stabilized by hydrogen bonds between a carbonyl oxygen and an amine hydrogen. These hydrogen bonds occur at regular intervals of one hydrogen bond every fourth amino acid and cause the polypeptide backbone to form a helix [1] . The most common helical structure is a right-handed helix with its hydrogen bonds parallel to its axis. The hydrogen bonds are formed between carbonyl oxygen and amine hydrogen groups of four amino acid residues away. Each amino acid advances the helix, along its axis, by 1.5 Å. Each turn of the helix is composed of 3.6 amino acids therefore the pitch of the helix is 5.4 Å. There is an average of ten amino acid residues per helix with its side chains orientated outside of the helix. Different amino acids have different propensities for forming x-helix, however proline is a helix breaker because proline does not have a free amino group. Amino acids that prefer to adopt helical conformations in proteins include methionine, alanine, leucine, glutamate and lysine (malek).

2. β-sheet: ß-sheets are stabilized by hydrogen bonding between peptide strands. In a β-sheet, regions of the polypeptide backbone come to lie parallel to each other and are connected by hydrogen bonds [1] . The hydrogen bonds are formed between the carbonyl oxygen and the amine hydrogen of amino acid in adjacent strands in a polypeptide, which means that the hydrogen bonds are inter-stand. β-sheet regions are more extended than an α-helix, and the distance between adjacent amino acids is 3.5 Å. Hydrogen bonding in β-strand can occur as parallel, anti- parallel, or a mixture. Amino acid residues in β- parallel configuration runs in the same orientation. Pleated sheets makes up the core of many globular proteins and also are dominant in some fibrous proteins such as a spiders web [1] . The large aromatics such as: tryptophan, tyrosine and phenylalanine, and beta-branched amino acids like: isoleucine, valine, and threonine prefer to adopt β-strand conformations.This orientation is energetically less favorable because of its slanted, non-vertical hydrogen bonds. Trytophan, tyrosine, and phenylalanine are hydrophobics while the other amino acids are hydrophilics.

3. β-turns: Poly peptide chains can change direction by making reverse turns and loops. Loop regions that connect two anti-parallel β-strands are known as reverse turns or β-turns. These loop regions have irregular lengths and shapes and are usually found on the surface of the protein. The turn is stabilized by hydrogen bond between the backbone of carbonyl oxygen and amine hydrogen. The CO group of the residue, in many reverse turns, which is bonded to the NH group of residue i + 3 . The interaction stabilizes abrupt changes in direction of the polypeptide chain. Unlike the alpha-helices and ß-strands, loops do not have regular periodic structures. However, they are usually rigid and well defined. Since they loops lie on the surface of the proteins, they are able to participate in interactions between proteins and other molecules. Ramachandran plot is a plot that shows the available torsion angles of where proteins can be found. However, in the plot, if there are many dots that locate all over the place, it means that there exists a loop.

Tertiary: As the secondary structure becomes established due to the primary structure, a polypeptide folds and refolds upon itself to assume a complex three-dimensional shape called the protein tertiary structure. Tertiary structure is the overall shape of a polypeptide. [1] Tertiary structure results from the interactions between the side chains (R groups) of the various amino acids [1] . This three dimensional structure is due to intramolecular interactions between the side groups along the polypeptide chain. Its domain typically contains 300 – 400 amino acids, and it adopts a stable tertiary structure when it is isolated from their parent protein. As a polypeptide folds into its functional shape, amino acids that have hydrophobic side chains tend to end up clustered at the core of the protein so that they are out of contact with water [2] . Covalent bonds called disulfide bridges can also affect the shape of a protein [1] . Disulfide Bridges form where two amino acids containing sulfhydryl groups on their side chains are brought close together by how the protein is folding [1] . For some proteins, such as ribonuclease, the tertiary structure is the final structure of a functional protein. Other proteins are composed of two or more polypeptides and adopt a quaternary structure.

Quaternary: While all proteins contain primary, secondary and tertiary structures, quaternary structures are reserved for proteins composed of two or more polypeptide chains [1] . Proteins that have quaternary structures contain more than one polypeptide and each adopt a tertiary structure and then assemble with each other via intermolecular interactions. The quaternary structure of a protein is the overall structure that is the result of the addition of these polypeptide subunits [1] . The individual polypeptides are called protein subunits, which means different polypeptides folded separately. Subunits may be identical polypeptides or they may be different. When proteins consist of more than one polypeptide chain, they are said to have quaternary structure and are also known as multimeric proteins, meaning proteins consisting of many parts. Quaternary structures can also defined as when more than one protein come together to create either a dimer, trimer, tetramer, etc. [2] . Hemoglobin is an example of a quaternary structure that is composed of two alpha subunits and two beta subunits.

Fibrous proteins: Fibrous proteins also known as Schleroprotein are long protein chains shaped liked rodwires. Unlike Globular Protein, they do not denature as easily, and contain many repeats of secondary structures. They are mostly structural proteins that are responsible for organisms in support and protection such as forming connective tissue, muscle fibers, bones, and tendons . The two examples of fibrous proteins are:

1. α –keratin: α –keratin (essential in hair, hooves, horn, fingernails, and etc.) is a coiled-coil protein composed of two intertwining α-helices. Coiled-coil structures are found in other structural proteins, for example, the myosin of skeletal muscle it has heptads repeats correspond to 3.5 amino acids per turn. Residues in the position of a, d, a’ and d’ in the helices of these proteins are usually hydrophobic. The two strands in a coiled-coil are held together by hydrophobic interaction as well as ionic interactions and disulfide bonds.

2. Collagen: Collagen (of tendon, cartilage, blood vessel walls) is the most abundant protein in human’s body. Collagen is a triple helix that is unlike α-helix, it has 3.3 amino acids and 10 Å per turn. Collagen is stabilized by hydrogen bonds, which is formed between the carbonyl oxygen and the amine hydrogen of amino acids situated on neighboring chains and is perpendicular to the fiber axis. It is abundant in proline, and contains hydroxyproline and hydroxlysine. However, due to the abundance of proline, there are no intrachain of hydrogen bonds, and the hydroxylation of proline and lysine requires Vitamin C. Vitamin C deficiency causes scurvy. One third of amino acids of collagen are glycine because of overcrowding only glycines are found in the center of the collagen molecules. Collagen molecules can be cross-linked by covalent bonds to from larger fibers and sheets.

Globular Protein: Globular proteins are folded to bury the hydrophobic side chains. All globular proteins have an inside where the hydrophobic core is arranged. It has an outside toward which the hydrophilic groups are directed. The uncharged polar amino acid residues are usually found on the protein surfaces but it can also occur in the interior. In the latter case, it will hydrogen bonded to other groups, i.e. ser, thr, tyr are all polar, uncharged.

Several factors determine the way that polypeptides adopt their secondary, tertiary and quaternary structures. The amino acid sequences of polypeptides are the defining features that distinguish the structure of one protein from another. As polypeptides are synthesized in a cell, they fold into secondary and tertiary structures, which assemble into quaternary structures for most proteins. As mentioned, the laws of chemistry and physics, together with amino acid sequence, govern this process. Five factors are critical for protein folding and stability:

1. Hydrogen bonds: Hydrogen bonds are formed between a hydrogen bond donor and hydrogen bond acceptor. For amino acids, hydrogen bonding would occur between the backbone of the amine group and the oxygen of the carbonyl group.

2. Ionic bonds: Electrostatic interactions occur between two oppositely charged molecules. Ionic interactions are weaker in water than in vacuum, this is due to a different dielectric constant faced in water between opposing charges within the protein's structure.

3. Hydrophobic effect: The hydrophobic interaction originates from the tendency of non-polar molecules to minimize their interactions with water. When non-polar molecules interact with water, these molecules tend to cluster together in the center to form a micelle.

4. Van der waals forces: Van der waals forces exist between non-polar molecules at close range. Of the three van der waals interactions, interactions between permanent dipoles is the strongest, dipole-induced dipole interactions are weaker than permanent dipole and the London dispersion forces are the weakest. While van der waals forces between individual atoms are weak, the sum of van der waals forces resulting from interactions between many atoms in large macromolecules can be substantial. The strength of van der waals interactions varies with the distance between the atoms and is maximal at the van der waals contact distance.

5. Disulfide bridges: A disulfide bond can be form between two cysteines through oxidation. These are also the strongest covalent bonds within a protein's tertiary structure.

Protein denaturation: is the loss of native conformations of tertiary structure. Denaturing proteins experience either the destruction of disruption of internal tertiary or secondary structure. Denaturation however, does not break the peptide bond between adjacent amino acids, thus not affecting the primary structure of the protein. Denaturation however, will interfere the normal alpha-helix and beta sheets in a protein which ultimately distort its 3D shape.

Denaturation causes the disruption of hydrogen bonding between close proximity amino acids, thus interfering a protein's secondary and tertiary structure. In tertiary structure there are four types of bonding interactions between "side chains" including: hydrogen bonding, ionic bridges, disulfide bonds, and hydrophobic intermolecular interactions. In other words, there are several different conditions to denature the conformation of a protein.

Conditions that denature proteins:

1. Extreme pH (pH < 4 or pH > 9) : alters H-bonding

2. Heat (temp >70oC): thermal effect, disrupts weak forces of non-covalent bonds

3. Detergents or organic solvents : disrupts hydrophobic interaction

4. Chaotropic agents (high concentrations) : e.g., urea and guanidinium chloride

As scientists started discovering more aspects of chemistry, they've actually found the magnitude of complexity in cell chemistry/biology. Although scientists found out that protein had an imperative role in the body, they've also discovered that the proteins assemble themselves at a specific site in the cell, being activated only when necessary. Using the GFP- tagged proteins (fluorescence) in animate cells, the positioning and repositioning of proteins were observed in response to the specific signals. When extracellular signal molecules bind to the receptor proteins, it reels in different proteins towards the inner area of the plasma membrane to create protein apparatus that will pass on the signal.

Humans have 10 PKC enzymes that differ both in their regulation and in their functions. When the PKC gets activated it will move from the cytoplasm to various intracellular locations and will eventually form specific complexes with other proteins thus allowing them to phosphorylate different protein substrates. Various ligases express this kind of behavior such as SCF ubiquitin ligases. These mechanisms involve the collaboration between protein phosphorylation and scaffold proteins that link specific activating, inhibiting, adaptor, and substrate proteins to a discrete part of a cell.

This occurrence is called induced proximity, which describes the reason why minute different forms of the enzymes with the same reaction sites can have different functions. This can be done by covalently modifying the protein’s location in various ways. These alterations construct binding sites on the proteins so that it would bind to scaffold proteins, making them cluster together so that different reactions can take place within a specific location of a cell. Scaffolds therefore allow the cells to group reactions without the need of membranes.

Scaffold proteins were thought to hold the proteins in specific locations relative to each other but in reality, unstructured regions of polypeptide chains connect the proteins that are interacting. This allows the proteins to frequently clash with each other in random orientations, some leading to successful reactions. The tethering of the proteins allows faster reaction rates to occur. Scaffold proteins therefore provide flexible methods of controlling the Cell Chemistry.

DEAD box proteins consist of RNA helicases, they are involved in RNA metabolism processes, and they are conserved in nine domains found in bacteria and viruses to humans. They are 350 amino acids in length. DEAD box proteins are involved in pre-mRNA processing, splicesosome formation, and rearranging of ribonucleoprotein (RNP) complexes. DEAD box proteins are required in the pre-mRNA splicing and the in vivo splicing process. During the pre-mRNA processing, the DEAD box proteins unwind to provide energy to rearrange the five snRNPs (U1, U2, U4, U5, and U6) required in pre-mRNA splicing. In the in vivo splicing, three DEAD box proteins, Sub2, Prp28, and Prp5, are needed. Prp5 helps rearrange the conformation of U2, which allows the U2 sequence bind to the branch point sequence. Prp28 helps the recognition of the 5’ splicing location.

The first DEAD box protein, the ElF4A translation initiation factor, are dependent on RNA ATPase activity. This protein helps unwind the secondary structure, which stops the scanning


Conclusions

It has been shown that protein protomers that form hetero-oligomeric complexes tend to have structures more similar to each other than proteins that do not form this type of supramolecular assemblies. A series of different approaches have contributed to this observation: distances on the proteomic Ramachandran plots, protein structure superpositions, and comparisons based on two domain structure databases (CATH and SCOP).

In agreement with previous studies, it is reasonable to suppose that this surprising similarity between protomers of hetero-oligomeric complexes is due to the evolutionary relationship between hetero-oligomers and earlier homo-oligomers, though gene duplication and paralogs evolution (Archibald et al. 1999 Ispolatov et al. 2005 Lukatsky et al. 2007 Lukatsky et al. 2006 Pereira-Leal et al. 2007). However, in my opinion, further studies are necessary to evaluate the relative importance of evolutionary and physico-chemical restraints on protein structure and dynamics.


Studying backbone torsional dynamics of intrinsically disordered proteins using fluorescence depolarization kinetics

Intrinsically disordered proteins (IDPs) do not autonomously adopt a stable unique 3D structure and exist as an ensemble of rapidly interconverting structures. They are characterized by significant conformational plasticity and are associated with several biological functions and dysfunctions. The rapid conformational fluctuation is governed by the backbone segmental dynamics arising due to the dihedral angle fluctuation on the Ramachandran ϕ–ψ conformational space. We discovered that the intrinsic backbone torsional mobility can be monitored by a sensitive fluorescence readout, namely fluorescence depolarization kinetics, of tryptophan in an archetypal IDP such as α-synuclein. This methodology allows us to map the site-specific torsional mobility in the dihedral space within picosecond-nanosecond time range at a low protein concentration under the native condition. The characteristic timescale of

1.4 ns, independent of residue position, represents collective torsional dynamics of dihedral angles (ϕ and ψ) of several residues from tryptophan and is independent of overall global tumbling of the protein. We believe that fluorescence depolarization kinetics methodology will find broad application to study both short-range and long-range correlated motions, internal friction, binding-induced folding, disorder-to-order transition, misfolding and aggregation of IDPs.

This is a preview of subscription content, access via your institution.


Background

Biofuels are a clean and renewable source of energy, rising as an alternative to fossil fuels, such as those derived from petroleum [1, 2]. They are produced from agricultural materials, for example, sugarcane, corn, soil, seaweed, and so on [3]. Second-generation biofuel production occurs in several steps, such as pre-processing, saccharification, and fermentation. The saccharification step occurs by the synergistic action of three types of enzymes: endoglucanases (E.C. 3.2.1.4), exoglucanases, also called cellobiohydrolases (E.C. 3.2.1.91), and β-glucosidases (E.C. 3.2.1.21) [4, 5]. Endoglucanases act in the cellulose structure, releasing oligosaccharides of different lengths. Cellobiohydrolases hydrolyzes the terminal of these oligosaccharides, releasing mainly cellobiose molecules. Then, β-glucosidases hydrolyzes the cellobiose glycosidic bond, releasing two glucose molecules [4,5,6,7]. However, most β-glucosidases are strongly inhibited by high glucose concentrations [8,9,10]. Thus, these enzymes have been considered by several studies as targets to improve high glucose concentrations tolerance by site-direct mutagenesis or the design of new enzymes [8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42]. Also, many reviews have reported the importance of glucose tolerance for improving the saccharification process [4, 7, 43].

Recently, Salgado et al. [43] proposed a β-glucosidase classification system divided into four groups: (i) β-glucosidases strongly inhibited by glucose (most of them) (ii) β-glucosidases tolerant to glucose (iii) β-glucosidases stimulated by low glucose concentrations but inhibited in high concentrations and (iv) β-glucosidases not inhibited by high glucose concentrations. To the best of our knowledge, the groups ii, iii, and iv are composed of few enzymes. Therefore, many studies aimed to transfer their characteristics to other non-efficient enzymes for biomass hydrolysis. For example, Yang et al. [9] evaluated the importance of a set of amino acid positions through site-direct mutagenesis. They reported that H228T and N301Q/V302F mutations could lead a marine non-resistant β-glucosidase to glucose tolerance. Also, Giuseppe et al. [10] reported that shape and the presence of hydrophobic residues in the middle of the substrate channel could be related to the structural basis of glucose tolerance. Furthermore, mutations in the positions 174, 404, and 441 of a β-glucosidase extracted from the Turpan Depression metagenome, have been reported as necessary for increasing the optimal temperature and reduce the optimal pH [12]. The study of Cao et al. [12] demonstrated that the β-glucosidase of the Turpan Depression metagenome could be classified as glucose-tolerant. However, the wild enzyme presented a low Kcat/Km value when using cellobiose as substrate. Also, the half-life of the wild enzyme at 50 °C was only 1 h. Therefore, this could hinder the employment of this enzyme in cellulose hydrolysis. The combination of three beneficial mutations (W174C/A404V/L441F) was essential to extending the half-life to 48 h, keeping the IC50 and, consequently, the glucose tolerance. The use of the mutant enzyme allowed an improvement of the sugarcane bagasse conversion by 14–35%, which demonstrated that multiple aspects should be considered to propose mutations that improve the activity of β-glucosidases.

Computational approaches have also been used in the search for crucial amino acids to convert non-tolerant to tolerant β-glucosidases. For instance, a set of 15 mutations have been proposed to improve the activity of a non-tolerant β-glucosidase from a marine metagenome [44]. From these 15 proposed mutations, a previous study has provided experimental evidence of enhancing β-glucosidase activity even in high glucose concentrations for three of them: H228C, H228T, and H228V [9]. The residues mutated V302F, N301Q/V302F, F172I, V227M, G246S, T299S, and H228T were also the target of other computational studies that used classic and accelerated molecular dynamics simulation to highlight their role in glucose releasing [45, 46]. Despite all these efforts, the rational design of more efficient β-glucosidases is still a challenge.

Previously, a database containing structures of glucose-tolerant β-glucosidases, called Betagdb, has been proposed [4]. Betagdb database was developed based on papers that reported glucose-tolerant β-glucosidases with experimental validations and structural data from public databases (only 23 occurrences were found at that moment). With the rising and popularization of next-generation sequencing platforms, thousands of β-glucosidase from several organisms were stored in sequence databases, such as UniProt. These data could be better explored to bring new insights into β-glucosidase mechanisms. In this paper, we propose a database of β-glucosidases enzymes called Glutantβase. Our database includes 3842 sequences collected from UniProt of β-glucosidases from the GH1 family (Glycoside Hydrolase Family 1), the most promising family for second-generation biofuel production. For all sequences, we performed comparative modeling, predicted their secondary structure, detected the residues involved in coevolution networks, detailed the conserved residues, the catalytic glutamates, and the residues present in the substrate channel that guides to the active site. Also, we hypothesized that mutations described in the literature as beneficial for improving β-glucosidase activity could be extrapolated to other β-glucosidases. To verify this, we modeled 5607 mutant proteins based on analogous positions of six beneficial mutations described in the literature: H228T [9], V174C [12], A404V [12], L441F [12], H184F [27], and E96K [47]. We performed molecular docking of glucose and cellobiose in the wild and mutant proteins to verify the affinity score variation. Our results show that only mutations in analogous positions of H228T impact in the interactions of glucose and cellobiose, which agree with previous computational and experimental studies [9, 44, 45]. We hope Glutantβase might help engineering tolerant β-glucosidase enzymes to bring improvements in second-generation biofuel production.


Methods

Chemical-Shift-Based Prediction of Secondary Structure Propensities.

In the first step of the CHESHIRE procedure, chemical shifts are used to predict the secondary structure of the protein. The method that we developed, termed 3PRED, uses Bayesian inference to predict the secondary structure of amino acids from the known chemical shifts in combination with the intrinsic secondary structure propensity of amino acids triplets The probability distributions P δ measure the likelihood for individual amino acids of forming specific secondary structures S given a set of experimentally measured chemical shifts (δH α , …, δC β ). The second set of probability distributions P 3 take into account the intrinsic propensities of fragments of three consecutive amino acids (Q 1, Q 2, Q 3) to form given secondary structures (S 1, S 2, S 3). The P 3 distributions act as smoothing potentials to increase the accuracy of the assignments derived from chemical shifts alone through the P δ distributions.

The propensities P 3 were computed by considering all of the structures in the ASTRAL SCOP database (35) having <25% sequence identity according to the secondary structure classification provided by the program STRIDE (36). For the calculations of the probabilities P δ, chemical shifts were calculated by applying SHIFTX (17) to the same set of structures to obtain an extensive database (3PRED-DB), which consisted of 939,639 calculated chemical shifts for each atom type.

Once the probabilities P 3 and P δ are known, for computational convenience they can be recast into pseudoenergies as Thus, the pseudoenergy E of a secondary structure assignment S for a protein of sequence Q and chemical shifts Δ can be approximated as The most likely secondary structure S and the single propensities (P H, P B, P C) are then computed by averaging the assignments with the pseudoenergy function E. We used a Monte Carlo scheme in which E is minimized by a search in the space of the N-dimensional vectors S in which at each move the secondary structure assignment of a single amino acid is changed. Predictions were obtained by considering 10 6 such steps at a pseudotemperature T = 1.

Chemical-Shift-Based Prediction of Dihedral Restrains: TOPOS.

In the second step of the CHESHIRE procedure, the secondary structure propensities computed by 3PRED are used as input in TOPOS, an algorithm based on an approach similar to that of TALOS (2), to predict the backbone torsion angles that are most compatible with the experimental chemical shifts. In TOPOS, for each protein segment of three residues centered at position i in the sequence (the target), the similarity to a triplet centered at position j in a sequence in the ASTRAL SCOP database (the source) is evaluated by computing the similarity function σ(i, j) where Δδ is the secondary chemical shift of a given atom of the source and target protein segment the parameters k h and k s were both set to 0.2, and the values of the remaining parameters and of the amino acid similarity matrix ΔResType were taken from Cornilescu et al. (2). The first terms in Eq. 3 are similar to the TALOS scoring function, the only substantial difference being that we do not consider H N chemical shifts. By contrast, the term k s log P n+j(S n+j) is the secondary structure bias present in TOPOS but not in TALOS. To avoid overfitting problems due to the use of a limited database, TOPOS uses the same extensive database of 3PRED.

The fragments with the highest σ scores, typically 200–500, are then clustered together according to the distance of the backbone torsion angles of the central amino acid. Finally, the average dihedral Φ and Ψ angles for the three best-scoring clusters are reported as prediction.

Prediction of the Structures of Fragments.

The CHESHIRE method is based on the molecular fragment replacement approach, which has been shown to be successful for the determination of protein structures with RDC (27) and in ab initio structure determination (37). In the present method, two types of fragments, of three and nine amino acids, respectively, are selected from the ASTRAL SCOP PDB database. The scoring function takes into account three contributions: (i) the score E shifts between the experimental chemical shifts of the fragment of the protein considered and the chemical shifts of the structure in the database, (ii) the score E restr for the compatibility with the dihedral angle restraints obtained with TOPOS, and (iii) the score E secstr for the match between the predicted secondary structure and the secondary structure of the fragment where the weights are set as

Chemical-shift score.

The chemical-shift score used in the fragment selection is similar to the score used by TOPOS, the only differences are that (i) the ΔResType is not included and (ii) the effect of residues i − 1 and i + 1 on residue i are not taken into account. where E shift(i, j) is given by

Dihedral angle restraint score.

The term E restr penalizes fragments that have torsion angles that are incompatible with the predictions of TOPOS. A fragment is compatible if its distance, on the Ramachandran plot, with at least one of the predicted values is <60°.

Secondary structure score.

The secondary structure score penalizes database segments with secondary structures that differ from those predicted by 3PRED: where P(S j, i) is the probability to have the secondary structure assignment S j at position i.

This step of the CHESHIRE procedure provides at each position along the sequence ten fragments of length three and five fragments of length nine. These fragments are used to generate the low-resolution structures, as described below.

Generation of Low-Resolution Structures.

Molecular representation.

In the initial low-resolution structure generation, a coarse-grained representation of the protein chain was used in which only backbone atoms are explicitly modeled (H, N, C α , C′, O) side chains are represented by a single C β atom. Bond lengths and angles, and the ω backbone torsion angle are kept fixed, while the Φ and Ψ torsion are given the freedom to move.

Energy function.

The energy function used for the low-resolution structure generation is a linear combination of terms that model different features of folded proteins: In the following text, we illustrate the meaning of these energy terms.

Pairwise interactions.

E vdw, E elec, and E EEF1 model van der Waals, electrostatic, and solvation, respectively. The first two were adapted from the CHARMM PARAM19 (38) and the third from ref. 39. The pairwise potential of mean force E PMF was implemented by using all known PDB structures in the ASTRAL SCOP database following Zhou and Zhou (40).

Secondary structure packing.

To model correctly the packing of secondary structure elements, the potential of Baker and coworkers (41) (E SS, E SH, and E HH) was implemented.

Cooperative hydrogen bonding.

This term (E CHB) was implemented according to ref. 42 to favor the formation of β-sheets by β-strands distant in sequence.

Structure generation protocol.

Low-resolution structures were generated by using a Monte Carlo algorithm carried out in an extended configuration space Γ given by the Cartesian product of the protein chain coordinates and a “virtual secondary structure” string where N and M are, respectively, the numbers of atoms and amino acids in the protein chain. These M additional discrete degrees of freedom are used to switch on and off energy terms that depend on the secondary structure of the protein.

Starting from a fully extended chain, conformations are generated by 20,000 Monte Carlo moves using a simulated annealing protocol. Two kinds of moves are applied. In the first (fragment substitution), the torsion angles and the secondary structure string in a randomly selected three- or nine-residue window of the protein chain are replaced with those from a fragment of known structure. In the second, local backbone moves, the torsion angles, but not the secondary structure, of a window of four amino acids are randomly perturbed. The score of the new conformation is calculated, and the move is accepted according to the Metropolis criterion. For each of the proteins studied here, 10,000 trial structures were generated in this way.

Refinement.

Molecular representation.

In the third stage of the CHESHIRE procedure, all atoms, including polar hydrogen atoms, are represented explicitly from the trial structures generated from the previous low-resolution stage. In a first phase, bond lengths, angles, and the ω backbone torsion angles are kept fixed, while the Φ, Ψ, and side chain torsion angles are let free to move. Structures are then optimized by using the energy function described below. Finally, the best-scoring structures are further refined by repeated minimizations and side chain optimizations using the Dunbrack and Cohen rotamers library (43).

Initial structures were obtained by adding the missing atoms to the low-resolution structures according to the following protocol. (i) A fully extended all-atom protein chain is generated by using ideal geometries. (ii) Target Φ and Ψ angles are set to those of the source chain. (iii) An energy minimization of 10,000 steps is performed to remove steric clashes. (iv) An additional energy minimization of 10,000 steps is performed by restraining interbackbone distances to the original ones. (v) A final energy minimization of 10,000 steps is performed without any restraint.

Structure screening.

All structures containing steric clashes as well as those with a radius of gyration larger than R max = 2.83 × M 0.34 , where M is the number of amino acids in the protein (44), were discarded.

Energy function.

The CHESHIRE energy function is a combination of a physicochemical term (E FF) and of a term that describes the correlation (C) between experimental and predicted chemical shifts: where E FF is a background force field given by and log(1 + C)capp is given by where Here, corrX is the correlation between the experimental and the back-calculated chemical shifts for atoms of type X, k ha = 18, and k n = k ca = k cb = 1. The term C is capped at 3.5 to avoid correlations between experimental and back-calculated chemical shift exceeding the error of SHIFTX. With this choice of values, the correlations are biased until they reach a threshold of ≈0.8 for H α atoms and 0.9 for N, C α , and C β atoms.

Force field.

All terms in E FF except E hb are the same defined in Eq. 10 the E hb term models backbone hydrogen bond following Kortemme et al. (45).

Chemical-shift correlation capping.

The chemical-shift correlation term C is capped at 3.5 to avoid correlations between experimental and back-calculated chemical shift that are better that the error of SHIFTX. With this choice of values, the correlations are biased until they reach a threshold of ≈0.8 for H α atoms and 0.9 for N, C α , and C β atoms.

Structure generation protocol.

After addition of the side chain atoms, the E scores of all structures were computed, and the best 500 structures were selected for refinement. The refinement consisted of a simulated annealing Monte Carlo run of 10,000 steps. The use of a Monte Carlo strategy enables us to use a bias on the chemical shifts without requiring the derivatives of the cost function as would be necessary in a molecular dynamics scheme. After refinement, structures were ranked according to their scores, and the best-scoring one was selected as the final result.