article: Inverse statistical approaches to determine protein structure and function from Multiple Sequence Alignments (MSA) are emerging as powerful tools in computational biology. However the underlying assumptions of the relationship between the inferred effective Potts Hamiltonian and real protein structure and energetics remain untested so far. Here we use lattice protein model (LP) to benchmark those inverse statistical approaches. We build MSA of highly stable sequences in target LP structures, and infer the effective pairwise Potts Hamiltonians from those MSA. We find that inferred Potts Hamiltonians reproduce many important aspects of ‘true’ LP structures and energetics. Careful analysis reveals that effective pairwise couplings in inferred Potts Hamiltonians depend not only on the energetics of the native structure but also on competing folds; in particular, the coupling values reflect both positive design (stabilization of native conformation) and negative design (destabilization of competing folds). In addition to providing detailed structural information, the inferred Potts models used as protein Hamiltonian for design of new sequences are able to generate with high probability completely new sequences with the desired folds, which is not possible using independent-site models. Those are remarkable results as the effective LP Hamiltonians used to generate MSA are not simple pairwise models due to the competition between the folds. Our findings elucidate the reasons for the success of inverse approaches to the modelling of proteins from sequence data, and their limitations.
  Prediction of protein structure from sequence remains a major goal of computational structural biology with significant practical implications. While great progress was achieved in de novo prediction of structures and even direct folding of smaller globular water soluble proteins [1,2], structure prediction remains challenging for larger proteins, and membrane and other non-globular proteins. For these cases indirect methods such as homology modeling are most promising. An approach to predict structure from the statistics of multiple sequence alignment (MSA) was proposed twenty years ago [3]. The underlying assumption for these statistical approaches is that residues that covary in a MSA are likely to be in close proximity in protein structure [4]. Recently, this approach has been significantly improved and has become a practical tool to extract structural information from sequence data [5,6], and to help in folding proteins [7]. Progress was made possible by the application of the maximum entropy principle of statistical mechanics [8,9] to derive the distribution of sequences in a protein family, i. e. their probability of appearing in an MSA. This approach is similar in spirit to the derivation of Gibbs distribution in Statistical Mechanics, with an effective Hamiltonian constructed to reproduce single-site amino-acid frequencies and pairwise amino-acid correlations. As such, the effective Hamiltonian is a Potts model [10] with parameters (the site-dependent fields and pairwise couplings acting on amino acids) fitted to reproduce the MSA statistics. This approach allows one to disentangle direct (corresponding to couplings) from indirect (mediated by other sites) correlations [6,11]; large couplings are much better predictors of contacts than large correlations. The exploitation of covariation information in proteins extends beyond structural prediction, and is potentially useful for homology detection, for characterizing the fitness changes resulting from mutations, or for designing artificial proteins with ‘natural’ properties [12]. Experiments show in particular that a sizeable fraction of artificial sequences generated to respect the 1- and 2-point amino-acid frequencies calculated from the natural MSA of the WW domain (a short protein domain with ≃30 amino acids) acquire the native fold [13,14]. This fraction vanishes when artificial sequences are generated which respect the pattern of single-site frequencies only. Recent studies have shown that inferred maximum-entropy models are helpful to predict the effects of mutations in various protein sequences [15–17], including HIV virus sequences [18]. Despite those successes fundamental issues remain poorly understood. To what extent do the inferred couplings reflect the physical (energetics) and structural (contact) properties of proteins? How much are the covariation properties of one family influenced by the need to prevent sequences from folding into competing structures? Are pairwise correlations, and hence, couplings generally sufficient to capture the statistical features of protein families and define generative models? In this work we answer those questions using lattice-based protein (LP) models [19–21]. Using LP allows us to evaluate the strengths and limitation of the Potts inference methods by comparing recovered solutions with known exact ones. Protein-like sequences are generated using a well-defined Hamiltonian that employs 20 amino-acid types and the Miyazawa-Jernigan energy function to approximate contact energies [22]. LP share many common properties with real proteins, for instance nontrivial statistical features of the sequences associated with a given fold, with variable conservation along the protein chain [20,23]. Surprisingly, the covariation properties of LP had not been much studied so far [24]. Here, we do so by generating MSA data for various LP folds, and apply to those sequence data the same inverse approaches used for real proteins (Fig 1). In particular, we show that the inferred couplings are excellent predictors of contacts in the protein structure. In addition, while the inferred Potts pairwise couplings mostly reflect the energetics and the contacts in the fold corresponding to the MSA, they also strongly depend on the nature of the other folds competing with the native structure, and have transparent interpretations in terms of positive and negative designs. Furthermore, we show that the pairwise Potts model is generative: it produces with high probability sequences with the right fold, a performance which cannot be achieved by models reproducing single-site frequencies alone. This is a non-trivial result as the log probability of LP sequences constrained to fold in a given structure is not a sum of pairwise contributions only, but includes multi-body interactions to all orders due to contributions from competing folds.
  We show in Fig 2A the positive predictive value (PPV) for contact prediction (Methods) for the four structures SA, SB, SC, SD, based on the ranking of the mutual information (MI) scores [26] and of the inferred Potts couplings, with the mean field (DCA) [6], the pseudo likelihood (PLM) [27], and the adaptive cluster expansion (ACE) [28–30] procedures. A fifth method, called Projection, shown with magenta lines in Fig 2A will be introduced later on. Mean-field DCA is a very fast, approximate method to infer the couplings. PLM is another approximation method, known to be remarkably accurate for contact prediction on real protein data. ACE is slower, but provides a more precise estimate of the coupling parameters (Methods). As in the case of real protein data, Potts-based contact predictions, either with DCA, PLM or ACE, generally outperform MI-based predictions. MI, indeed, does not disentangle large indirect correlations mediated by one or more sites from direct correlations due to contacts. PLM and ACE couplings are more precise than their DCA counterparts and accordingly, give better contact predictions for all four structures (Table 2). As shown later, most of the missed contacts with ACE (3 or 4, depending on the structure in Fig 2A) are carried by the site at the center of the structure. Why this is so and how to improve contact prediction (magenta lines in Fig 2A) will become clear from the detailed interpretation of the Potts couplings below. While the results above were obtained for uniformly sampled MSAs, real protein data suffer from imperfect sampling, e. g. resulting from phylogenetic correlations. To analyze the effects of biases in the sampling on structural predictions in a simple and controlled way we have introduced a sampling procedure that favors sequences with small Hamming distances to a reference, wild-type sequence (Methods). As shown in Fig 2B, modifying the strength of the bias, b, allows us to interpolate smoothly between the uniform sampling above (b = 0) and strongly biased MSAs with reduced sequence diversities (large b values). Results for contact prediction obtained with the PLM procedure are reported in Fig 2B. We observe that performances worsen as the bias in the MSAs increases. Empirical reweighting of the sequence statistics according to their similarity with the other sequences in the MSA (Methods), as usually done for real protein data, efficiently improves the quality of the structural prediction for intermediate bias values. For very strong biases (b = 0. 1 in Fig 2B), the effective number of sequences resulting from reweighting is so low that results become very poor. We now test the ability of the inferred Potts model to be generative, that is, to produce sequences having high probability of folding into the native structure. To do so, we infer the pairwise Potts model, Eq [9], from the MSA, hereafter referred to as ‘natural’, of structure SB with the ACE procedure. We sample the Potts distribution using Monte Carlo simulations, thereby generating a new MSA, hereafter referred to as ‘Potts-ACE’. We then compute the probabilities of folding (into SB) of all the sequences in the Potts-ACE MSA, see Eq [5]. Results are shown in Fig 3A. Strikingly, the ‘majority of the sequences have high folding probabilities. Conversely, sequences generated with an Independent-site Model (IM, see Methods) following the same procedure are unlikely to have the right fold, as shown in Fig 3A. Pairwise interactions are therefore crucial to design new sequences with the right fold. Because the Potts-ACE model is generative, its entropy can be used as a quantitative estimate of the designability of the structure, with results reported in Table 1 and S1 Text, Section I. B. To better characterize the generative properties of Potts-ACE, we show in Fig 3B the scatter plot of the ‘energies’ H P o t t s - A C E [ A; h, J ], Eq [9], with the inferred Potts-ACE model vs. the effective Hamiltonian -β log Pnat (SB|A) for the best folding sequences A in the MSA generated with the Potts-ACE model for structure SB. We observe a roughly monotonic relationship between the two quantities: the lower the energy with the Potts-ACE model, the more likely is the sequence to be a good folder (Fig 3B). However, while the Potts-ACE model rightly predicts that sequences with increasing energies are less likely to be good folders, the value of H P o t t s - A C E is less and less predictive about the precise value of Pnat, as shown by the increasing dispersion in Fig 3B. Note that the low-energy part of this scatter plot, corresponding to the best sequences generated with the Potts-ACE model, is quantitatively similar to its counterpart for the ‘natural’ MSA, see Fig F in S1 Text. We stress that Potts couplings have to be inferred with high precision to generate high-quality sequences. We show in Fig 3A that sequences generated by the Potts-Gaussian model, which makes use of the approximate DCA couplings (details in S1 Text, Section II. B), have very high folding probabilities, but are extremely concentrated around the consensus sequence. The Potts-Gaussian model, contrary to the IM and the Potts-ACE model, fails to reproduce the diversity of sequences observed in the natural MSA. This failure is a direct consequence of the Gaussian approximation (S1 Text, Section IV). The Potts-PLM model reproduces the diversity of natural MSA, but generates sequences, whose folding properties largely vary from very good to very poor (Fig 3A). The poor folding discrimination is mostly due to the large regularization parameter used in the inference (Section II. C and Fig B in S1 Text). In summary, the Potts-ACE model, in contradistinction with IM, DCA and PLM, is capable to generate a large set of diverse sequences that fold with high probability in the target native structure. We now study in more detail the properties of the Potts couplings. We show in Fig 4A the scatter plot of the Potts couplings Jij (a, b) inferred with the ACE method for structure SB vs. the Miyazawa-Jernigan energetics parameters −E (a, b) used to compute the LP energies (Methods); scatter plots for the other structures are shown in Figs H, I, and J in S1 Text. For each pair of sites i, j, we observe a linear dependency, J i j (a, b) ≈ λ i j × - E (a, b), (1) where the slope λij is positive for the pairs in contact in the native structure (red symbols in Fig 4A), and negative, or zero for the pairs not in contact (blue symbols). In the following, we characterize each pair i, j with its slope, computed through λ i j = - ∑ a b J i j (a, b) E (a, b) ∑ a b E (a, b) 2. (2) We interpret λij as a measure of the coevolutionary pressure on the sites i, j, due to the design of the native structure. This interpretation is supported by the following theoretical and approximate expression for the pressure (see Methods and derivation in S1 Text, Section III): λ i j T H ≈ c i j - c ¯ i j log 1 + 1 β N S e - Δ. (3) The numerator, c i j - c ¯ i j, measures the difference between the native contact map, cij (= 1 if i, j are in contact and 0 otherwise), and the average contact map, c ¯ i j, of the structures in competition with the native fold, each weighted by e−Δ, where Δ is the typical gap between the energy of sequences folded into the native structure and their energies in the competing structures (Methods) [31]. The native and average contact maps for structure SB are shown, respectively, in the lower and upper triangles of Fig 4B. The denominator in Eq [3] depends on the inverse sampling temperature, β, in the Monte Carlo procedure used to generate the MSA, the effective number, NS, of competing structures, and the typical energy gap, Δ (Methods and S1 Text, Section III. D). Fig 4C shows the pressures, λij, computed from the inferred couplings, Eq [2], vs. the contact-map difference, δ c i j ≡ c i j - c ¯ i j, across all pairs of sites i, j for the native fold SB. We observe a monotonic dependence of the pressure λij with δcij; in particular, λij has the same sign as δcij, except for a few, small (in absolute value) λij and δcij. Note that pairs i, j such that j − i is even valued can never be in contact in any structure on a cubic lattice: c i j = c ¯ i j = 0, as can be seen in the checkerboard pattern of c ¯ i j in Fig 4B. The associated pressures, λij, are weak and indicate the level of noise in the inference (magenta symbols in Fig 4C). A detailed classification of the pairs of sites is shown in Fig 4C, based on the identification of the closest competitor to SB, structure SF, represented in Fig 4D. The closest competitor is defined as the structure that is closest in energy (smallest gap) to the native structure, see Methods, Eq [11]. Informally speaking, SF is the ‘most’ dangerous structure as sequences folding in SB could fold in SF upon a few undesirable mutations. As the native fold is not in competition with a single structure, the notion of closest competitor is somewhat reductive, but is helpful to understand the mechanisms controlling covariation and the meaning of the inferred couplings. In particular, comparison of the contact maps of SB and SF allows us to define four classes of contacts, and to illustrate concretely the mechanisms for positive and negative designs [24,32]. Pairs with large and positive δcij are in contact in the native fold, but not in the competitor structures. They are under strong covariation pressure to stabilize the native fold and not the competitors (positive design). We see in Fig 4C that the pairs of sites in contact in SB and not in SF (class Unique Native, UN) are associated with large couplings. On the contrary, pairs with small and positive δcij are in contact both in the native and competitor folds, and are subject to weak pressures. As a result, the corresponding couplings are weak, and contact prediction is difficult. Pairs of sites in contact in both SB and its closest competitor, SF, form the Shared Native (SN) class in Fig 4C. For the structures we have studied, those contacts are generally established between the site at the center of the cube and its four neighbors (see fold SB in Fig 2A and its contact map in Fig 4B) as these central contacts are typically shared with closest competitor folds (Fig 4D and Fig K in S1 Text). Sites i, j not in contact in the native fold correspond to negative δcij. Most of those pairs are not in contact in the competitor structures either and correspond to weak pressures and weak inferred couplings. This phenomenon is illustrated by the Non-Native (NN) class, which gathers the pairs of sites that are neither in contact in SB more in SF, see Fig 4C. However, some pairs do correspond to contacts in the competitor folds (large and negative δcij), and are subject to strong and negative pressures. The anti-correlation between the couplings Jij (a, b) and the Miyazawa-Jernigan parameters −E (a, b) is the result of negative design: the native structure is favored by rendering contacts that appear only in competitor folds unfavorable [32]. Negative design is clearly visible for the pairs of sites in the Close Competitor (CC) class in Fig 4C, which are in contact in SF but not in the native fold, SB. As the scores used for contact predictions are based on the squared couplings (Methods), pairs of sites in the CC class, associated with large couplings (anti-correlated with the Miyazawa-Jernigan energy parameters), may have larger scores than true contacts in the SN class, associated with weaker couplings (positively correlated with MJ). To avoid those false positives, we rank pairs of sites according to their projection on the MJ matrix, defined in Eq [2] (Methods). This ranking method, referred to as ‘projection’, improves contact prediction, see PPV in magenta in Fig 2A. It allows us to detect, in particular, the contacts of the central site, which are often shared among the native and competing structures (see contact maps in Fig 4B for structure SB and its competitors), and which correspond to weak, positive pressures. We find, for the four structures studied, that 27 out of the 28 largest projections do correspond to real contacts (Table 2). The evolutionary pressures, averaged over the pairs of sites in the UN, SN and CC classes, give λUN ≈ 2. 8, λSN ≈ 1. 4, λCC ≈ −0. 9. Those results are comparable to, but smaller than, the average theoretical predictions from Eq [3]: λ U N T H = 4. 5, λ S N T H = 2. 7 and λ C C T H = - 1. 8. Larger values, closer to the theoretical estimates are obtained when the regularization parameter, penalizing large couplings, is made smaller (Methods). To further test the theoretical prediction of the pressure, Eq [3], we have varied two features of the Monte Carlo procedure used to generate the MSA: (1) the pool of possible competing structures appearing in the denominator of the folding probability, Pnat, in Eq [5], and (2) the inverse sampling temperature β (Methods). We observe that the pressures λij increase when the pool of competing structures is restricted to structures similar to the native one, or when β is increased, as sequences are constrained to have higher values for Pnat (Fig L and Table A in S1 Text).
  Lattice proteins offer a fully controlled, nontrivial benchmark to characterize the factors that affect the success of inverse statistical approaches in predicting structure from multi sequence alignments (MSA). The maximum entropy Potts model, inferred to reproduce the low-order statistics of an MSA, can be successfully used for structure prediction. Potts-based predictions, in particular with accurate sophisticated inference methods (PLM, ACE), outperform Mutual Information approaches on lattice protein data (Fig 2A), as for real protein data. An important finding in the present work is the approximate linear relationship between the inferred Potts couplings, Jij (a, b), and the Miyazawa-Jernigan energetic parameters, E (a, b), across the pairs of amino acids a, b for a given pair of sites i, j, Eq [1] and Fig 4A. The slope of this linear relationship, λij, is a measure of the covariation ‘pressure’ on the pair i, j due to the structural constraints. The pressure depends on (i) the difference between the contact maps of the native fold (cij) and the ones of the competing structures (c ¯ i j), see Fig 4B; (ii) the proximity (in terms of energy, Eq [11]) of the native fold with its competitor structures; (iii) the strength of constraints in the design, controlled here by the inverse sampling temperature (Eq [3]). A consequence of (i) is that the intensities of the inferred couplings, Jij (a, b), measured through the scores Fij in Eq [10], are directly related to the local properties of the contact map on sites i, j (and not on other sites), which explains the success of inverse approaches to disentangle direct from indirect effects in amino-acid correlations. However, the relationship also involves the average contact map c ¯, originating from the competitor folds. As c ¯, contrary to the native contact map c, is not sparse, neither are the couplings J. Couplings can, in particular, take non-zero values for i, j in contact in the competing folds (c ¯ i j close to 1), but not in the native fold (cij = 0), see Fig 4A and 4B. Those ‘repulsive’ couplings give rise to negative pressures λij, and prevent the sequence from forming contacts that would result in the wrong fold, a clear illustration of negative design [24,32,33]. Reciprocally, sites in contact in the native fold (cij = 1) need not be subject to strong covariation and associated to strong couplings, if they are also in contact in the competing folds (c ¯ i j close to 1). Our results clearly show that strong couplings do not necessarily imply, nor require structural contact (in the native fold). They also explain the success and limitations of contact prediction based on the magnitude of the inferred couplings (Fig 2A). A second major result of our approach, shown in Fig 3A and 3B, is that the inferred Potts model, when accurately inferred, here, with the ACE procedure, generates sequences that fold with high probability in the native state and show amino-acid diversity similar to the one observed in the original MSA. This performance is remarkable. While the energy of a LP sequence in a given structure is a sum of pairwise contributions only, the effective Hamiltonian of the sequences constrained to fold with high probability in this structure (and not in the competing folds) includes high-order multi-body interactions (originating from the competing folds) between amino acids, as can be expected for real proteins. Furthermore, independent-site models are unable to produce sequences with the correct fold (Fig 3A). Our result corroborates the works of Ranganathan and collaborators [13,14], who experimentally showed that sequences built according to a reshuffling of the MSA of the small WW domain respecting site conservation and (approximately) pairwise correlations folded correctly with a good (30%) probability. The generative character of the Potts model, combined with the linear relationship between couplings and energetics parameters, agrees with previous studies showing that designed and real protein-like sequences are generated by Gibbs measure with Hamiltonian reflecting real protein energetics [20,32]. We stress that couplings have to be inferred with great care, going beyond the mean-field DCA approximation, in order for the Potts model to be generative (Fig 3A). The analysis of real protein data could benefit from our analysis in several ways. First, the understanding of the relationship between the couplings and the contact map obtained here indicates that the pressures λij (Eq [2]) are excellent estimators of contacts (Fig 2A), and may outperform the usual scores, F i j (A P C) (Eq [10]). The use of λij allows one to both keep track of the sign of the couplings (and avoid false positives corresponding to pressures originating from negative design, see above) and increase signal-to-noise ratio, by removing noise in the inferred couplings not aligned along the Miyazawa-Jernigan energy matrix. In real proteins, the physico-chemical properties of amino acids suggest the existence of different classes of contacts, whose energetics could be inferred from large-scale analysis of databases and used to define class templates E. Projections of the inferred Potts couplings on those templates could in turn be used for contact predictions. Secondly, the detailed analysis of lattice proteins done here has allowed us to better understand how the existence of positive and negative design influences the covariation statistics of the MSA, and, hence, the properties of the inferred couplings. The relevance of negative design was already put forward in literature, e. g. in the study of the TIM-barrel fold protein [32]. A promising path to further investigate negative design could be the study of coexistence between natively folded and misfolded proteins [34]. The direct estimate of c ¯ in that case would be straightforward, as it would essentially correspond to the contact map of the misfolded competing structure. More generally c ¯ in real proteins can be determined from the ensemble of misfolds obtained either from all-atom folding simulations [2] or from threading approaches. This is the subject of future work. Last of all, our work suggests a practical way to quantify the designability of a protein family, i. e. to measure the number of sequences ‘belonging’ to the family. Current approaches rely on the computation of the maximal eigenvalue of the native contact matrix c [25], though our study suggests that the matrix c ¯, gathering information about the competing protein folds, is also a determinant of the designability. Given the generative character of the inferred Potts model, a direct and more accurate measure of the designability is provided by its entropy computed with the ACE procedure, see Methods and Table 1. This approach can be applied to real protein families, an example is provided in [35]. In future studies, the total control we have over lattice proteins will allow us to study the importance of covariation and the success of inverse approaches in ‘Gedanken’ experiments related to many important issues of interest for real proteins, and going beyond the structural aspects studied here. Among those let us cite the detection of homology. Current protein databases, such as PFAM, classify query sequences into families based on Hidden Markov Model profiles, an extension of independent-site models capable of taking into account deletions and insertions [36]. It is an open problem to understand whether coupling-based models, exploiting covariation, could, contrary to HMM, recognize sequences with low homology in a reliable and computationally tractable way. Another very important issue is the estimation of fitness landscapes, or fitness changes in responses to one or more mutations [37,38]. Covariation-based models have been recently introduced to predict escape paths of pathogens (virus or bacteria) from drugs or vaccines in this context [18]. Lattice proteins offer a unique benchmark to understand deeply and quantitatively and, ultimately, improve those approaches.
  We consider model proteins, whose L = 27 amino acids occupy the sites of a 3 × 3 × 3 cubic lattice [19–21]. Four of the 103,406 possible configurations of the backbone (excluding symmetries), hereafter called folds or structures, are shown in Fig 2A. Unless otherwise said, we restrict ourselves to a representative subset of N = 10,000 folds [39]. Two amino acids are said to be in contact if they are nearest neighbors on the lattice (but not on the backbone). The contact map c (S) of structure S is its 27 × 27 adjacency matrix: c i j (S) = 1 if i, j are in contact, 0 otherwise. Amino acids in contact interact through the Miyazawa-Jernigan (MJ) 20 × 20 energy matrix E [22]. The energy of a sequence A = (a1, …, a27) of amino acids folded into structure S is E (A | S) = ∑ i < j c i j (S) E (a i, a j). (4) The Boltzmann probability that sequence A folds into structure S is given by (for unit temperature): P nat (S | A) = e - E (A, S) ∑ S ′ = 1 N e - E (A, S ′). (5) Good folders are sequences A with large gaps between E (A | S) and the energies E (A | S ′) with the other folds [20]. We generate a multi-sequence alignment (MSA) for the native fold, say, S, through Monte Carlo (MC) simulations with the Metropolis rule [32]. The simulation starts from a sequence A, and attempts to mutate randomly one amino acid at a time, say, a i → a i ′; Let A′ denote the new sequence. If Pnat (S|A′) > Pnat (S|A) the mutation is accepted; otherwise it is accepted with probability [Pnat (S|A′) / Pnat (S|A) ]β (<1). Parameter β plays the role of an inverse algorithmic temperature that sets the stringency of sequence selection in the sampling MC procedure. The corresponding effective Hamiltonian, H LP [ A | S ] = - log P nat (S | A), (6) includes not only the native-structure energy E (A, S) in Eq [4], with pairwise interactions between the amino acids ai, but also contributions coming from all the other folds S′ in Eq [5], with multi-body interactions at any orders ≥ 2. In practice we choose β = 103, which ensures that sequences fold in the native structure S with probability 0. 995 or larger and that thermalization is fast (Fig A in S1 Text). A conservative estimate of the equilibration time of the MC dynamics (defined as the time needed for Pnat to relax to its equilibrium value independently of the initial sequence condition) is < 1,000 MC steps. The MSA is then made of the sequences Aτ, with τ = 1, …, M = 50,000, generated at regular intervals of 1,000 steps. Each sequence Aτ is therefore drawn according to the equilibrium measure equal to Pnat (S|Aτ) β, up to a multiplicative normalisation factor. To mimic the biases in real protein MSAs, e. g. resulting from phylogeny, we may add to the effective Hamiltonian a penalty term proportional to the number of mutations with respect to a wild-type sequence A (wt), H biased LP [ A | S, b ] = - log P nat (S | A) − b ∑ i δ a i, a i (wt), (7) and repeat the MC sampling procedure above, with M = 100,000. High values of the bias strength b produce MSAs with high sequence identity, see Fig 2B. We compute the one- and two-point amino-acid empirical frequencies computed from the sequences Aτ in the MSA of a given native structure, f i (a) = 1 M e f f ∑ τ = 1 M ω τ δ a, a i τ, f i j (a, b) = 1 M e f f ∑ τ = 1 M ω τ δ a, a i τ δ b, a j τ, (8) where δ denotes the Kronecker function. The weight ωτ sequence Aτ is defined as the inverse of the number of sequences in the MSA whose sequence identities with Aτ are larger than 27 × (1 − w), where w is an arbitrary cut-off. The effective number of sequences is M e f f = ∑ τ = 1 M ω τ. Throughout this work no reweighting was imposed (w = 0 and ωτ = 1), except for the biased sampling case, where some results were obtained with cut-off w = 0. 2, see Fig 2B. We consider the maximum entropy (least constrained [8]) distribution over the sequences, that reproduces the one- and two-point empirical frequencies. This maximum entropy distribution is the Gibbs measure of a Potts model, with q = 20-state variables, whose energy is given by H Potts [ A; h, J ] = - ∑ i h i (a i) - ∑ i < j J i j (a i, a j), (9) where the sets of L × q fields, h = {hi (a) }, and of L (L − 1) /2 × q2 couplings, J = {Jij (a, b) }, are chosen so that the frequencies computed from the Potts model distribution reproduce those of the data, Eq [8]. Note that, as the distribution depends on the energies of the sequence configurations through their relative differences the Potts parameters are defined up to an arbitrary gauge (S1 Text, Section II). We resort to three methods to solve this hard computational problem. Within the mean-field Direct Coupling Approximation (DCA) [6], J is approximated as minus the pseudo-inverse of the connected correlation matrix c, with entries cij (a, b) = fij (a, b) − fi (a) fj (b). Data are regularized through the introduction of a large pseudocount (S1 Text, Section II. B). The Pseudo-Likelihood Method (PLM) is a more sophisticated approximation to compute the direct couplings J, exploiting the full MSA (and not only the frequencies). A large regularization L2-penalty, γ = 0. 01, was considered during the inference [27] (S1 Text, Section II. C). The Adaptive Cluster Expansion (ACE) is an accurate but slower procedure, which recursively builds the most relevant clusters of strongly interacting sites [28–30]. Couplings were inferred in the least-probable gauge, in which each interaction parameters corresponding to the least frequent amino acids are set to zero, with a weak L2–penalty γ = 5/M; M is here the number of the sequences in the MSA. See S1 Text, Section II. D for more information about the implementation of the methods, and the ability of the Potts-ACE distribution to reproduce high-order statistics of the MSA. In addition to h, J ACE gives access to the entropy of the inferred Potts model, which is an estimate of the designability of the native fold. For the simpler Independent-site Model (IM), the energy H IM is a sum of field contributions only (J = 0 in Eq [9]), see S1 Text, Section II. A for details. Computation of the Potts couplings J with DCA, PLM or ACE allows us to define the scores [6] (with the Averaged Product Correction [40]) F i j (A P C) = F i j − F i. F j. F. . , where F i j = ∑ a, b J i j (a, b) 2, (10) and Fi. = ∑l Fil, F. . = ∑i Fi. . The values of the scores depend on the choice of the gauge. For the PLM and DCA procedures it is customary to choose the so-called zero-sum gauge, in which the sums of the coupling matrix elements Jij (a, b) over each column or row vanish. This gauge is the one minimizing the Frobenius norm. For the ACE procedure, we have used the couplings inferred in the least-probable gauge, see paragraph above and S1 Text, Section II. D, to compute both the Frobenius (leading to black dots in Fig 2A) and the projection (magenta symbols) scores; we have verified that contact predictions do not change for the four structures when scores are computed from couplings expressed in the zero-sum gauge. Within the mutual information approach, scores are simply given by MIij = ∑a, b fij (a, b) log [fij (a, b) / (fi (a) fj (b) ) ] [26] (and are then subject to the above APC correction), and do not distinguish direct from indirect correlations. Scores are then ranked in decreasing order, and used to predict contacts. The positive predictive value at rank k is the fraction of the top k scores whose corresponding pairs of sites (i, j) are in contact. We define the gap Δ (S|Snat) between structure S and the native fold Snat through e−Δ (S|Snat) =〈 e−[ E (A| S) −E (A |Snat) ] 〉A, (11) where the average is taken over sequences A in a MSA of Snat. The distributions of gaps for Snat = SA, SB, SC, SD are shown in Fig G in S1 Text. The closest competitor to Snat is defined as the structure S minimizing Δ (S|Snat), see S1 Text, Section III. We approximate the number NS of competitor folds and their typical gap Δ with the native structure Snat through N S e − Δ ≡ ∑ S (≠ S n a t) e − Δ (S | S n a t). The average contact map of those competitors is defined as c ¯ i j ≡ 1 N S e - Δ ∑ S (≠ S n a t) e - Δ (S | S n a t) c i j (S).
summary: Inverse statistical approaches, modeling pairwise correlations between amino acids in the sequences of homologous proteins across many different organisms, can successfully extract protein structure (contact) information. Here, we benchmark those statistical approaches on exactly solvable models of proteins, folding on a 3D lattice, to assess the reasons underlying their success and their limitations. We show that the inferred parameters (effective pairwise interactions) of the statistical models have clear and quantitative interpretations in terms of positive (favoring the native fold) and negative (disfavoring competing folds) protein sequence design. New sequences randomly drawn from the statistical models are likely to fold into the native structures when effective pairwise interactions are accurately inferred, a performance which cannot be achieved with independent-site models.
keywords: sequencing techniques
  split-decomposition method
  multiple alignment calculation
  mathematics
  forecasting
  statistics (mathematics)
  protein structure prediction
  protein structure
  molecular biology techniques
  thermodynamics
  research and analysis methods
  sequence analysis
  entropy
  sequence alignment
  proteins
  mathematical and statistical techniques
  statistical methods
  molecular biology
  physics
  protein structure comparison
  biochemistry
  computational techniques
  biology and life sciences
  physical sciences
  protein sequencing
  macromolecular structure analysis
title: Benchmarking Inverse Statistical Approaches for Protein Structure and Design with Exactly Solvable Models