Selection for Folding Stability Predicts Observed Covariation Between Protein Positions in the PDB

Abstract

b002Statistical couplings between protein sites have been enormously useful for protein structure prediction, but there is still debate on the selective forces that generate them. Here we explicitly predict the couplings considering selection for protein folding stability, both against unfolding and misfolding, which we model through the Stability constrained model of protein evolution (SCPE). The SCPE is based on contact interactions, and it predicts stability against misfolding through the Random Energy Model in the space of contact matrices. We adopt the minimum selection principle, assuming that the multivariate amino acid distribution has minimum Kullback-Leibler divergence with respect to a site-unspecific distribution that is not influenced by protein folding stability. We predict the couplings at first order, assuming that they are small (small coupling approximation, SCA), and obtain an explicit formula that relates them with the contact interaction matrix and with the structural properties of the pairs of sites: presence of native contact, distance along the sequence and buriedness of the sites. We test these predictions on a representative set of the Protein Data Bank. For pairs of sites in contact, both predicted and observed couplings are inversely related with the contact energies of the amino acid pairs, while for pairs not in contact this relation is positive, i.e. negative design destabilizes pairs of amino acids for short-range pairs of sites that can form wrong contacts with high probability. Our results suggest that the strongest couplings correspond to native contacts, but an even larger number of couplings are influenced by negative design. Intriguingly, the most informative pairs are formed by surface sites, suggesting a role for protein-protein interactions and functional motions. We interpret the site-unspecific couplings, which are not influenced by selection for protein folding stability, as the result of variations of global mutational or selective forces across the lineages where the proteins evolve. Consistently, the strongest global couplings are self-couplings between equal amino acids. Cys-Cys, which form disulfide bridges, show the strongest global coupling, followed by pairs of metabolically costly amino acids (Trp-Trp, Tyr-Tyr, Trp-Tyr) and positively charged amino acids that tend to interact with nucleic acids (Lys-Lys, Arg-Arg). Ser, Thr, Asn and Gln form a cluster of globally correlated polar amino acids.

Introduction

b009The evolution of protein sequences is commonly modelled mathematically as a substitution process, i.e. a stochastic process that describes the probability that an amino acid at a given site in a protein is substituted by another one as the function of time. This formalism is very powerful for inferring phylogenetic trees and other evolutionary properties (Yang 2006). Substitution processes typically assume that different sites of the protein evolve independently of each other, which allows feasible computation of the likelihood function required for evolutionary inferences. However, it is known that the pairwise correlations that are neglected in this approach are key for predicting structural interactions between protein sites (Göbel et al. 1994; De Juan et al. 2013). The inference of these evolutionary correlations through neural networks lead to spectacular success in protein structure prediction (Jumper et al. 2021). Large part of these evolutionary correlations are thought to arise from the selective pressure for maintaining the functional structure of the protein and its folding stability (Zeldovich et al. 2007; Goldstein 2008; Liberles et al. 2012; Bastolla et al. 2017). Within this framework, one of us and co-workers developed a model that predicts the stationary amino acid distribution of a protein family with known native structures (Bastolla et al. 2006), from which we derived the stability-constrained protein evolution (SCPE) model, a substitution process for phylogenetic inference based on selection on protein folding stability (Arenas et al. 2015), which we explicitly predict based on a statistical physics model of the misfolded state (Minning et al. 2013). We later extended the SCPE model to enforce structure conservation, since several lines of evidence suggest that mutations that change the native structure are selected against more strongly than mutations that maintain the structure and only affect its folding stability. Following (Echave 2008), we predict the effect of mutations on the protein structure as the linear response of the protein represented as an elastic network model. We obtain in this way the structure and stability constrained protein evolution (SSCPE) models (Lorca-Alonso et al. 2025).

b010For computational treatability, the substitution processes derived from the SSCPE models assumed that protein sites evolve independently of each other. The model takes into account structural interactions effectively, through the effects of mutations on folding stability and structural changes. The SSCPE models produce site-specific amino acid distributions that present qualitative differences at sites buried in the protein core, which tend to be more hydrophobic, more conserved and evolve more slowly, with respect to exposed sites that tend to be polar, have higher sequence variation and evolve more rapidly. These differences arise from the different constraints imposed by natural selection for folding stability, and they are consistent with observations from multiple sequence alignments (MSA) (Echave et al. 2016).

b011Adopting the SCPE model, we relaxed the assumption that protein sites evolve independently and assume that the stationary amino acid distribution presents pairwise couplings described by the Potts model (Wu 1982) pioneered by Lapedes et al. (2002); Weigt et al. (2009); Marks et al. (2011) for inferring evolutionary correlations and studied by several groups (for instance, Morcos et al. 2011; Skwark et al. 2013; Ekeberg et al. 2013; Kamisetty et al. 2013).

b012The study of the pairwise stability constrained model in the framework of the "small-couplings" approximation (SCA) described below was one of the subjects of the PhD thesis of Jonas Minning (Minning 2012), directed by Markus Porto, at that time at the University of Darmstadt, with the external advice of Ugo Bastolla. Unfortunately, despite reiterate attempts, we have not been able to get in contact with Markus Porto, which is the reason why he does not appear among the authors although he is clearly one of the main contributors to this work. Should he read these lines, we will be more than happy if he gets in touch with us, and he is included in the authors list, as he deserves.

b013Here we consider the pairwise SCPE model within the SCA approximation, which provides an analytic formula that we test here, leaving the test of the general theory for future work. Despite its simplicity, the SCA constructively supports the hypothesis that site-specific selection for protein folding stability, together with global mutational and selective forces that act at the level of the whole protein, may produce the observed covariation among protein positions.

Multivariate Amino Acid Distribution

b016Following (Weigt et al. 2009), we model the L-dimensional probability distribution that sites \(A_i\), \(i\in [0,L-1]\) in a protein family are occupied by the amino-acids \(a_i\) as the exponential distribution

b017We use the notation \(x_i= \nu i +a_i\), which has \(L\nu\) possible values where L is the number of sites and \(\nu\) is the number of amino acid types (21, since we consider gaps as an additional symbol). In the following, we indicate by sum over \(a_i\) the sum over the \(\nu\) amino acids at position i, and by sum over \(x_i\) the sum over all L sites and their possible amino acids \(a_i\). The single site parameter vector \(h_{x_i}\) has \(L\nu\) components and the pairwise coupling matrix \(g_{x_i,x_j}\) is an \(L\nu \times L\nu\) symmetric matrix whose elements with \(j=i\) vanish.

b018In the physical literature, Eq.(1) is called "Potts model" (Wu 1982), and it has been extensively studied in the context of disordered systems before being recognized as the exponential pairwise distribution that describes a protein family.

b019In general, we can not normalize the Potts distribution either analytically or numerically, because computing the normalization factor \(Z=\sum _{a_0\cdots a_{L-1}} P(a_0,\cdots a_{L-1})\) requires an exponential number of operations. Thus, its parameters are usually inferred adopting approximated methods, such as the message passing algorithm originally proposed by Weigt et al. (2009), the regularized matrix inversion proposed by Jones et al. (2012), or the pseudo-likelihood approximation (Besag 1975), applied to this problem by Aurell and coworkers (Ekeberg et al. 2013) and (Kamisetty et al. 2013), among others.

Propensities

b021For notational convenience, we define the propensity matrix

b022that represents the ratio between the conditional probability of \(a_i\) conditioned to \(a_j\) and the unconditioned probability, minus 1. It is zero if the amino aids are independent, positive if amino acid \(a_j\) at site j favors \(a_i\) at site i, and negative if it affects it negatively. We call this quantity “propensity” between the two variables \(a_i\) and \(a_j\). This name is used in bioinformatics to describe the propensity of a given amino acid to take any secondary structure (Koehl and Levitt 1999). We also call \(\log (1+Q)\) “logarithmic propensity”. At first order in Q, the propensity and the logarithmic propensity coincide: \(\log (1+Q)\approx Q\).

b023The mutual information between sites i and j, \(M_{ij}\), is related to the propensities through

b024In the limit of small propensities the mutual information is approximated by the weighted mean square propensity.

Zero-Mean Gauge

b026The distribution Eq.(1) presents more parameters than necessary, in the sense that different parameter sets produce exactly the same distribution. These extra parameters must be fixed through normalization conditions or "gauge" (Weigt et al. 2009; Ekeberg et al. 2013). We change parameters as \(p_{x_i}\equiv \exp (h_{x_i})\), \(q_{xi,xj}\equiv \exp (g_{x_i,x_j})-1\) and write the multivariate distribution as

b027We adopt the following gauge conditions:

b028The sum over \(a_i\) denotes the sum over all possible \(\nu\) amino acids at position i. We call these conditions the "zero-mean" gauge because it sets to zero the weighted mean value of the couplings over the amino acids at one of the site, and we adopt it throughout this paper.

b029Note that the propensities satisfy the zero mean condition \(\sum _{a_j}p_{x_j}Q_{x_i,x_j}=\sum _{a_i}p_{x_i}Q_{x_i,x_j}=0\).

Global Distribution

b031We model the complete pairwise distribution through the combination of site specific parameters that represent natural selection for folding stability and global parameters \(p^\textrm{glob}_a\), \(q^\textrm{glob}_{ab}\) that are the same for all sites of the protein. These global parameters reflect the genetic code and the mutational process, which may be different in different genomes and different genomic regions. Indeed, considering the codon level instead of the amino acid level improves the inference of the couplings (Jacob et al. 2015). The global parameters also represent selective pressures that are uniform across the protein (see Discussion).

b032We represent the global pairwise distribution as \(P^\textrm{glob}_{a,b} = p^\textrm{glob}_{a} p^\textrm{glob}_{b} \left( 1 +q^\textrm{glob}_{a,b}\right)\) for all pairs of sites. It is easy to see that the global propensities \(q^\textrm{glob}_{a,b}\) obey the zero-mean gauge with respect to \(p^\textrm{glob}_{a}\).

b033We determine the \(\nu +\nu ^2\) global parameters \(\{p^\textrm{glob}_a\},\{q^\textrm{glob}_{a,b}\}\) through a regularized fit, as explained in Methods.

Site-Specific Couplings: Stability Constrains and Minimal Selection

b035The site-specific components of the modelled distribution reflect selection for maintaining the native structure of the protein and its folding stability. We determine them through the minimum selection principle (Arenas et al. 2015). This principle assumes that the complete distribution is as similar as possible to the global distribution, which is not influenced by selection for correct folding, in the sense that it has minimal Kullback-Leibler (KL) (Kullback and Leibler 1951) divergence from the global distribution, while constraining the average protein folding free energy \(\Delta G\).

b036The minimum selection principle generalizes the commonly adopted maximum entropy principle (see for instance Lapedes et al. 2002; Bastolla et al. 2006) if the amino acids are not equivalent. In fact, maximum entropy is equivalent to minimum KL divergence from the uniform distribution. Since the genetic code, the mutational process and the unspecific selective pressures differentiate the amino acids, the minimum KL divergence is preferable.

b037This principle can be justified based on the analogy between molecular evolution and the Boltzmann distribution of statistical physics (Sella and Hirsh 2005; Mustonen and Lässig 2010; Puente-Sánchez et al. 2024). In particular, the evolutionary process, modelled as a stochastic process in genotype space (Fisher 1930), converges to a stationary exponential distribution analogous to the Boltzmann distribution in statistical physics, with the logarithm of the fitness playing the role of minus the energy and the inverse of the effective population size playing the role of temperature, meaning that small populations are more tolerant towards sequences with low fitness.

b038In the present context, the inverse of the selective temperature, or selective strength, is represented by the Lagrange multiplier \(\Lambda\) that enforces the mean value of the protein folding free energy \(\Delta G\). This is the only free parameter of the stability constrained model. We fit \(\Lambda\) by maximizing the log-likelihood of the observed protein sequences with respect to the model.

b039Additionally, we impose the gauge conditions Eq.(5,6,7) through the Lagrange multipliers \(\mu _i\), \(\zeta _{x_i,j}\) and \(\eta _{x_j,i}\). These are not free parameters, since they are determined by imposing the constraints.

b040We determine the site-specific Potts parameters \(\{p_{x_i}\},\{q_{x_i,x_j}\}\) of the stability constrained model by minimizing the objective function \(\Phi\), which is the sum of the KL divergence plus the functions that describe the constraints.

b041The matrix \(G_{x_i,x_j}\) represents the effective free energy of the pair \(x_i\) and \(x_j\) with respect to unfolding and incorrect folding (misfolding). These are not free parameters: We predict them with our model of protein folding stability, which considers the native contact matrix and the statistics of non-native contacts (Minning et al. 2013), as we describe in section 3.4. This computation is approximate since we approximate the free energy of the misfolded ensemble as a sum of pairwise terms, while the REM prediction involves terms that depend on the amino acid at four protein sites.

b042We compute the optimal values of the site-specific parameters \(p_{x_i}\) and \(q_{x_i,x_j}\) by minimizing the score \(\Phi\), where \(\Lambda\), \(\{p_a^\textrm{glob}\},\{q^\textrm{glob}_{a,b}\}\) are pre-set global parameters and the external information comes from the predicted values of the effective free energies \(G_{x_i,x_j}\).

b043We impose the constraints Eq.(5,6,7) through the Lagrange multipliers \(\mu _i\), \(\zeta _{j,x_i}\) and \(\eta _{i,x_j}\), whose values are determined by the corresponding constraints. In this way, we determine the site-specific couplings and the single-site parameters as the function of the selective strength \(\Lambda\).

b044Finally, we determine the optimal value of the parameter \(\Lambda\) by maximizing the likelihood of the empirical pairwise distributions observed in the PDB for all types of pairs of sites.

b045In our previous work (Arenas et al. 2015), we determined the single site parameters \(p_{x_i}\) assuming that all sites are independent, i.e. assuming that \(q_{x_i,x_j}\equiv q^\textrm{glob}_{a,b}\equiv 0\). Here we relax this assumption, as we also did in Minning (2012) within the small couplings approximation (SCA) discussed below.

The Small-Coupling Approximation (SCA)

b047In the numerical tests presented below, we adopt the small-coupling approximation (SCA) that neglects second and higher order terms in the couplings in the formula of the pairwise probability, i.e. we only consider the first order terms in \(q_{x_i,x_j}\). Because of the zero-mean gauge, the resulting values of the single site probabilities are \(P_i(a_i)=p_{x_i}\). The propensities Eq.(2) are

b048i.e. at the first order in q the propensities are equal to the direct couplings, and we neglect the correlations that arise from indirect couplings, although they can be recovered a-posteriori as we do when we predict the propensities of indirect contacts.

b049The simplicity of the SCA approximation stems from the fact that the couplings of the site of pairs ij only depend on its properties of the pair, and not on the other pairs. This approximation is reminiscent of the first approach that was adopted to infer the couplings \(q_{x_i,x_j}\) through the mutual information, which in the small coupling limit is equal to the weighted mean square \(Q_{x_i,x_j}\). The seminal work by Weigt et al. (2009) and its sequels clearly showed that the approximation based on the mutual information is much less accurate for inferring the couplings than the approach based on the Potts model that also considers the indirect couplings. Nevertheless, we consider here the SCA because of two reasons: (1) The SCA provides an analytic formula that is easy to interpret and sheds light on the selective forces that are behind the couplings. (2) The resulting pairwise distributions depend only on the “local” properties of each pair ij, which allows grouping pairs with similar properties and performing the simple and exhaustive numerical tests that we describe below.

b050Using the approximation (9) and computing the derivatives of Eq.(8) with respect to the parameters \(\{p_i\}\) and \(\{q_{ij}\}\), we can explicitly determine the values of the parameters that minimize \(\Phi\) as

b051In the above equations, we denote with \(\textrm{ZM}\left( G_{x_i,x_j}\right)\) the zero-mean of the matrix \(G_{x_i,x_j}\), whose weighted mean value across the amino acids \(a_i\) and \(a_j\) is zero. Eq.(10) is the same that we obtained in our previous mean-field approximation in which we neglected the couplings and considered sites as independent, i.e. \(q_{x_i,x_j}\equiv 0\). The same formula continues to hold also at first order in the couplings, which does not change the single-site probabilities.

Site-Specific Distributions

b054For assessing Eq. 11, we sampled pairs of sites from a representative subset of the protein data bank, which we grouped in various groups according to their properties.

b055Firstly, we distinguished four types of proteins according to sequence length (short, \(<150\) residues, and long) and average contact energy, related with hydrophobicity (hydrophobic, \(\overline{U}<0\) and polar). These groups allow us analyzing how selection for stability depends on the properties of proteins. Alternatively, we lumped together all proteins in the same group.

b056For every group of proteins, we grouped sites in two classes (buried or exposed) depending on their number of native contacts \(n_i=\sum _j C^\textrm{nat}_{ij}\), thus obtaining three types of pairs of sites: exposed-exposed, exposed-buried and buried-buried. Alternatively, we considered only one type of sites.

b057We then classified pairs of sites i, j according to their native contact \(C^\textrm{nat}_{ij}\): either pairs that form a native contact (C1), or pairs that form an indirect contact through another site (C2) or pairs that do not form direct or indirect contacts.

b058Finally, we grouped the distance of the sites along the sequence \(|i-j|\) (contact range). We found that five groups of contact ranges are sufficient to obtain good results, and we show this value in most of the figures.

b059As a result, we obtained 3 pairs of sites times 3 types of contacts times 5 contact ranges, yielding 45 types of pairs for each group of proteins.

b060For each group of proteins and for each group of pair of sites p, we sampled the pairwise amino acid distribution from protein sequences in the PDB and computed the observed propensities \(Q^\textrm{obs}_{p}(a,b)=F_p(a,b)/(F_i(a)F_j(b))-1\) where \(F_p(a,b)\) is the observed frequency of the pair of amino acids a, b at pairs of sites that belong to class p. We obtained the single-site amino acid frequencies \(F_{i}(a), F_j(b)\) by marginalizing the observed pairwise distribution either symmetrically, if the sites are equivalent (core-core and surface-surface) or asymmetrically if they are different (core-surface).

Global Distribution

b062We determine the global parameters \(\{p^\textrm{glob}_a\},\{q^\textrm{glob}_{a,b}\}\) as follows.

b063First, we sampled the pairwise empirical distributions at all sites, without considering site specificity, i.e. we averaged the site-specific frequencies over all classes of pairs that belong to the same type of proteins. We thus get what we call \(F_{a,b}^\textrm{glob,obs}\), from which we obtain the single-site frequencies \(F^{\textrm{glob,obs}}_{a}\) by marginalization and the propensities as \(q^{\textrm{glob,obs}}_{a,b}=F^{\textrm{glob,obs}}_{a,b}/\left( F^{\textrm{glob,obs}}_{a}F^{\textrm{glob,obs}}_{b}\right) -1\).

b064This data is affected by the reduced sample size, which we treat as an overfitting problem. We fit the global parameters \(p^\textrm{glob}_a\) and \(q^\textrm{glob}_{a,b}\) from the observed parameters adopting a variant of the “rescaled ridge regression” method proposed by one of us and a coworker (Dehouck and Bastolla 2017). Precisely, we minimize the weighted mean of the errors of the fits. For propensities, we adopt the weights \(w_{a,b}=F^\textrm{glob,obs}_a F^\textrm{glob,obs}_b\) that satisfies the zero mean condition \(\sum _{a,b} w_{a,b}q^\textrm{glob,obs}_{a,b}=0\). We then add the L2 regularization term with regularization parameter R.

b065In the spirit of the rescaled ridge regression, we impose an additional condition on the average parameter values, \(\sum _{a,b} w_{a,b} q^\textrm{glob}_{a,b}\), for preventing that the scale of the fitted parameters vanish for infinite R. The additional term is formally similar to an L1 regularization with the signed value of the fitted parameters \(q^\textrm{glob}_{a,b}\) instead of the absolute value. For notational convenience, we indicate with c any pair of amino acids a, b. The scoring function that we minimize is

b066We determine the value of the regularization parameter \(\mu\) by minimizing the mean square error of the fit for given R. This yields \(\mu =R\sum _{c} w_{c} q_{c}^\textrm{glob,obs}\equiv R\langle q^\textrm{glob,obs}\rangle\) (see Appendix). In other words, the rescaled ridge regression sets the scale of the fitted parameters that minimizes the error of the fit, and it only imposes the regularization on the relative values of the parameters, leaving the fit unregularized when there is only one parameter.

b067Finally, we apply the specific heat criterion (Dehouck and Bastolla 2017) and we determine the regularization parameter R by maximizing the derivative of the error with respect to R, which represents a specific heat, in the statistical mechanics analogy proposed by Dehouck and Bastolla. This yields \(R=0.5\) (see Appendix).

b068The final result is \(q_c=\frac{2}{3}q^\textrm{obs}_c+\frac{1}{3}\langle q^\textrm{obs}\rangle\). For the propensities, because of the zero mean gauge, it holds \(\langle q^\textrm{obs}\rangle =0\). For the single-site probabilities \(p^\textrm{glob}_a\), we choose weights \(w_a= p^\textrm{obs}_a\), and \(\langle p^\textrm{obs}\rangle =\sum _a (p^\textrm{obs}_a)^2\). We finally obtain

Assessment of the Predictions

b070We then predict the pairwise distribution for each group of pairs of sites adopting Eq.(11). We use the average value of \(\left\langle C_{ij} \right\rangle\) over all pairs in the group and the observed single-site frequencies, since we want to focus on the prediction of the pairwise couplings.

b071We considered three scores that quantify the agreement between the observed and predicted pairwise amino acid distributions. These scores are related with the log-likelihood scores that we used for optimizing the parameter \(\Lambda\). The additional scores allow a more transparent assessment, and we show them in the figures.

b073where \(F_{p}^\textrm{obs}(a,b)\) is the frequency of the amino acid pairs a, b observed at the pair of type p, \(Q_{p}^\textrm{pred}(a,b)\) is the predicted propensity, and \(w_{i,j}\) is the weight of pair p (see below). The log-likelihood is related with the Kullback-Leibler divergence, except for a term that does not depend on \(Q_{p}^\textrm{pred}\), so that optimizing either score yields the same results:

b075This score is quadratic, and it is more sensitive to outliers. Its value allows us to assess the relative error of the predicted logarithmic propensities.

b076weighted Pearson correlation coefficient \(\textrm{wr}\) between the set observed and predicted logarithmic propensities, \(\log (1+Q_{p}^\textrm{obs}(a,b))\) and \(\log (1+Q_{p}^\textrm{pred}(a,b))\), in which each pair of amino acids a, b is again weighted proportionally to their observed frequency \(F_{p}^\textrm{obs}(a,b)\). We show this quantity because it is easy to interpret.

b077In all of these scores, we weight each pair of amino acids proportionally to its observed number of instances \(F_{p}^\textrm{obs}(a,b)\), which down weights infrequent pairs that present large relative error, improving the reliability.

b078Finally, we fitted the selection parameter \(\Lambda\) separately from each group of proteins by maximizing the weighted sum of the log-likelihood of all pairs p, using as weight of each pair its mutual information, \(w_{p}=M_{p}\equiv \sum _{a,b}F_{p}^\textrm{obs}(a,b)\) \(\log (F_{p}^\textrm{obs}(a,b)/F_i^\textrm{obs}(a)F_j^\textrm{obs}(b))\), so that we weight more the more informative pairs.

b079We performed the maximization through an iterative quadratic interpolation algorithm that evaluates the score for three values of the parameter, interpolates the score with a quadratic curve and determines the maximum, then it substitutes the candidate maximum to one of the three previous values and iterates until convergence.

b080After fitting the selection parameter, we compute the total score for each group of proteins using as weight \(w_{p}\) of each group of pair the number of sampled pairs. In this way, the final scores are not affected by the way in which the pairs are grouped, and their comparison is meaningful.

Stability with Respect to Misfolding

b082The strength of the stability constrained model consists in the fact that we predict the folding stability \(\Delta G\) without any free parameter by adopting a simple contact free energy model. In this model, protein structures are represented as binary contact matrices \(C_{ij}\in \{0,1\}\), with residues i and j considered in contact if any of their heavy atoms are closer than 4.5Å. For a given sequence \(a_i\), we estimate the effective free energy of the native state as \(G^\textrm{nat}=\sum _{i<j}C^\textrm{nat}_{ij}U(a_i,a_j)\,,\) where \(C^\textrm{nat}_{ij}\) is the contact matrix of the average native structure, \(a_i, a_j\) are the amino acids at positions i and j in the sequence, and U(a, b) are the effective contact interaction parameters determined by Bastolla et al. (2000).

b083Importantly, the free energy parameters U(a, b) where not determined from the statistics of amino acids in the PDB as in the quasi-Boltzmann approximation (Sippl 1995), which would be circular with the approach outlined here, but they were determined by imposing that the native state of a representative set of PDB structures are maximally stable against non-native decoys.

b084We estimate the effective free energy of the unfolded state as \(G^\textrm{unf}=L S_U\), where \(S_U\) is the average configurational entropy of an unfolded residue. For the ensemble of wrongly folded (misfolded) compact conformations, we adopt the Random Energy Model (Derrida 1981) in the space of compact contact matrices, which approximate the free energy of the ensemble of misfolded conformation with the mean and the variance of their contact free energies. Although corrections to the REM that involve the third moment of the energies are not negligible (Minning et al. 2013), here we do not consider them for simplicity.

b085In the above formula, the brackets \(\left\langle \cdot \right\rangle\) represent the mean value over the set of the possible compact contact matrices of proteins with L residues sampled from native structures in the PDB. Thus, the first term represents the mean contact free energy of misfolded conformations and the second term represents their variance.

b086Because of the contact correlation terms \(U(A_i,A_j)U(A_k,A_l)\), the free energy of the misfolded ensemble depends on sets of four sites, not just pairs of sites. We therefore average these terms in order to obtain an approximated expression that only depends on pairs of sites:

b087To evaluate the above expression, we sample the moments of the distributions of contacts from the PDB as follows.

b088The contact frequency \(\left\langle C_{ij} \right\rangle\) represents the mean value of the contact between residues i and j in proteins of L residues. In order to sample it, we assume that this is the product of a function c(L) that only depends on the protein length and a function \(f(|i-j|)\) that only depends on the distance of the two residues along the chain, which we estimate from contact matrices in the protein data bank (PDB). \(\left\langle C_{ij} \right\rangle\) increases with L (because the surface to volume ratio decreases with L, decreasing the influence of surface residues that have fewer contacts) and it decreases with \(|i-j|\) (because short range contacts form more easily than long range ones). We sample these functions from observed contact matrices in the PDB, accelerating the subsequent computations.

b089In the second term, \(N_C = \sum _{i<j}C_{ij}\) is the total number of contacts. The average number of contacts per residue, \(\left\langle N_C \right\rangle /L\), proportional to c(L), increases with L due to the decrease of the surface to volume ratio. We compute numerically its correlation with the contact frequencies as \(\left\langle C_{ij} N_C \right\rangle \approx Lc(L)^2 f1(|i-j|)\), sampling the factors \(f1(|i-j|)\) prior to the computation.

b090We define \(n_i=\sum _j C_{ij}\) as the number of contacts formed by residue i. Furthermore, we assume that its average value is the same for all residues (disregarding the important but limited end effect: residues at the chain termini form fewer contacts), and we numerically estimate the correlations between residues at distance \(|i-j|\) along the sequence as \(\left\langle n_i n_j \right\rangle - \left\langle n_i \right\rangle \left\langle n_j \right\rangle \approx c(L)^2 f2(|i-j|)\), again sampling the factors \(f2(|i-j|)\) prior to the computation.

b091Finally, in the spirit of the mean-field approximation, we approximate the terms in Eq.() that depend on the sequence with their average value, except the pairwise terms \(U(A_i,A_j)\). The factor \(\overline{U}\) is the average energy of non-native contacts, which we estimate with the global distribution \(P^\textrm{glob}\):

b093is the average energy of residue \(a_i\) sitting at position i (and same for residue \(a_j\) sitting at position j), which we compute at the beginning of the computation based on the site-specific single site distributions \(p_{x_j}\) and the global amino correlations \(q^\textrm{glob}_{a,b}\), but neglecting the site-specific couplings \(q_{x_i,x_j}\).

b094If \(\overline{U}>0\), the free energy of the misfolded ensemble is larger than the free energy of the unfolded ensemble, and we consider the latter for computing the folding free energy \(\Delta G\):

b095As the result of these approximations, we can decompose the folding free energy as the sum of pairwise terms,

Data Sets

b097We analyzed proteins in the PDB grouped at the 70% sequence identity level (list available in the PDB web server), retaining only structures with more than 40 residues studied by X-ray crystallography, which amounts to 37108 protein chains. In order to accelerate the computations, we produced files that contain the sequences and the contact matrices

Influence of the Misfolding Model

b100First, we optimized the temperature parameter T of the Random Energy Model (REM) applied to contact matrices, Eq.(18), which sets the scale of the free energy parameters U(a, b). We plot in Fig. 1 the scores summed over all pairs and all classes of proteins as the function of the parameter T. Fig. 1A and C show the log-likelihood score and Fig. 1B and D show the mean squared error of the logarithmic propensities. The top plots A,B do not consider the global propensities (\(q^\textrm{glob}_{ab}=0)\), which focuses more on folding stability, and the bottom plots C,D are obtained adopting \(q^\textrm{glob}\).

b101In all cases there is an optimum at intermediate T (\(T=0.8\) without \(q^\textrm{glob}\), \(T=0.7\) with \(q^\textrm{glob}\) for log-likelihood, \(T=0.9\) without \(q^\textrm{glob}\), \(T=0.8\) with \(q^\textrm{glob}\) for wMSE, in the arbitrary units of the contact interaction parameters).

b102The continuous lines in Fig. 1 show the scores computed considering only the native free energy and disregarding the average free energy of the misfolded ensemble, \(\sum _{ij}\left\langle C_{ij} \right\rangle U(A_i,A_j)\). The dashed line is the limit of infinite T, which is obtained neglecting the variance of the energy of the misfolded ensemble, U2.

b103The data show that it is important to consider selection against misfolding, which is usually neglected when modelling the folding free energy.

b104Agreement between observed and predicted pairwise propensities as the function of the parameter T of the REM. A, C: log-likelihood score. B, D: Mean Square Error. In plots A, B, the only fitted parameter was the selection strength \(\Lambda\). In plots C,D we also adopt the global propensity matrix \(q^\textrm{glob}\), which we obtain from the data

b105In all the following figures, we show data obtained at \(T=1.0\) and considering \(q^\textrm{glob}\), unless otherwise stated.

The Propensities are Related with the Contact Energy in Opposite Ways for Pairs with Native Contacts and Non-native Contacts

b107From Eq.(21) and (11), we see that the predicted propensities of amino acid pairs at pairs of sites p that form native contacts (\(C^\textrm{nat}_{p}=1\)) are negatively correlated with their contact energy U(a, b), i.e. amino acids that interact more favorably are predicted to be more frequent for pairs of sites in contact. On the contrary, for pairs of sites that do not form native contacts, the propensities are predicted to be positively correlated with U(a, b), i.e. amino acid pairs that repel each other are predicted to be more frequent, particularly so for short contact range for which \(\left\langle C_{p} \right\rangle\) is large (i.e. they can form contacts more easily).

b108From Eq.(21), we predict that the propensities are correlated with the zero-mean value of the energy parameters, defined as

b109We assess this expectation in Figure 2, which shows the scatter plots of the observed logarithmic propensities versus the zero-mean energy. We group together all proteins, and we distinguish pairs with direct, indirect and no native contacts and with short and long contact ranges.

b110Comparing the different types of contacts, we see that the propensities of native contacts present strong negative correlation with the contact energies, as predicted. This means that stabilizing pairs of amino acids are more frequent. The corresponding slope is much stronger for long range native contacts (Fig. 2D) than for short range native contacts (Fig. 2A), as expected since short range contacts are subject to “negative design” that tends to make them less favorable.

b111In agreement with negative design, short range pairs without native contacts (Fig. 2B) and with indirect contacts (Fig. 2C) are positively correlated with the contact energy, as predicted by our formula. This means that destabilizing pairs of amino acids are more frequent. However, this correlation disappears for long range pairs without contact, which are less subject to negative design (Fig. 2E) and it becomes negative for long range indirect contacts, which are more similar to native contacts than to no contacts (Fig. 2F).

b112We can see that the amino acid pair Cysteine-Cysteine (CC) has one of the largest propensities in all plots. This pair, which can form disulfide bridges, has the most stabilizing contact energy. However, its nature of outlier is not due to its effect on stability but it is due to the fact that CC has the largest global propensity, which we attribute to the mutation process and to global selective pressures (see below).

b113Scatter plots of the observed logarithmic propensities of amino acid pairs versus the zero-mean energy of the same amino acids, distinguishing pairs in contact (A, D), pairs without contacts (B, E) and pairs with indirect contacts (C, F). The first line represents short range contacts (A, B, C), the second line represents long range contacts (D, E, F)

b114Figure 3 plots the scatter plots of observed versus predicted propensities. They represent the same type of pairs as Fig. 2.

b115We can see that the SCA approximation allows reasonable predictions of the logarithmic propensities. The quality of the prediction is higher for long range native contacts (Fig. 3D) than for short range native contacts (Fig. 3A) that are subject to the contrasting forces of negative and positive design.

b116For sites without contact, the prediction is much better for long range pairs that are not subject to negative design (Fig. 3E). This is due to the fact that the prediction in this case is exclusively based on the global propensities \(q^\textrm{glob}\), which may suffer from overfitting despite the regularization that we adopted.

b117Indirect contacts show lower agreement (Fig. 3C and F). In all cases, the amino acid pair CC presents one of the highest propensities.

b118Observed versus predicted logarithmic propensities of amino acid pairs for pairs of sites in contact (A, D), without contacts (B, E) and with indirect contacts (C, F). The first line represents short range contacts (A, B, C), the second line represents long range contacts (D, E, F)

Dependence on Contact Range and Buriedness

b120In Fig. 4, we show how the contact range \(|i-j|\) influences the weighted correlation coefficients between the logarithmic propensities of each amino acid pair, \(\log (1+Q^\textrm{obs}(a,b))\), and the zero-mean of the free energy matrix, Eq.(22), on which our theoretical prediction is based. The figure depicts pairs of sites in contact (top row), not in contact (middle row) and in indirect contact (bottom row). Each column represents one of the four classes of proteins (Short hydrophobic, short polar, Long hydrophobic and Long polar) and each symbol represent a type of pairs of sites (surface-surface, surface-core and core-core).

b121The correlations are negative for pairs of sites in contact (top row, Fig. 4A, B, C, D), as expected based on positive design, but they are less negative for short range contacts as expected based on negative design, since short range contacts have large contact frequency \(\langle C_{ij}\rangle\), i.e. they are frequently formed in misfolded conformations. The strength of the correlations and their dependence on contact range is larger for pairs of sites both in the core (triangles) than pairs on the surface (circles). Correlations are relatively weaker for short polar proteins (Fig. 4B).

b122For pairs of sites that do not form native contacts (middle row, Fig. 4E, F, G, H) the correlation between propensities and contact energies is positive for short range contacts, as expected based on negative design, and it becomes negative for long range contacts. However, these negative correlations at long ranges are not attributable to positive design, as for pairs of sites in contact, but to the negative correlation between the contact energies and the global propensities \(q^\textrm{glob}\), which present positive propensities between pairs of hydrophobic amino acids. In fact, \(q^\textrm{glob}\) determines almost completely the propensities of long range pairs without contacts (Fig. 3E). In long proteins (Fig. 4G,H), the strongest correlations for non-contact pairs happen for pairs of site in the core (triangles), as it is the case for contact pairs, while for short proteins (Fig. 4E,F) they happen for surface-surface sites (the core is not so well defined for short proteins, and many of these proteins are forming complexes, but quaternary contacts are not taken into account in our analysis).

b123Similar to no-contact pairs, indirect contact pairs (Fig. 4I,J,K,L) show positive correlations at short range, as expected based on negative design, and negative at long range, where correlations are more negative than for no-contact pairs, i.e. long range pairs that form indirect contacts tend to interact more favourably than pairs that do not form contacts, which suggests that some kind of positive design is taking place for indirect contacts.

b124weighted Pearson correlation coefficients between the propensities of each amino acid pair, \(Q^\textrm{obs}(a,b)\), and the zero-mean of the free energy matrix, Eq.(22), on which our theoretical prediction is based, as a function of the contact range \(|i-j|\), for pairs in contact (first row, plots A,B,C,D), pairs not in contact (middle row, E,F,G,H) and pairs with indirect contact (bottom row,I,J,K,L). The three curves show three different pairs of sites. The four columns show the different classes of proteins

b125We show in Fig. 5 the quality of the predictions, assessed through the weighted Pearson correlation coefficients between observed and predicted propensities, which are based on four global propensities matrices, one for each group of proteins. As in the previous figure, the three rows represent the three types of contacts (native contact, no contact and indirect contacts) and the four columns represent the four groups of proteins.

b126Comparing the contacts, one can see that the predictions are generally worse for indirect contacts (bottom row), which are not well represented in the SCA approximation. Comparing pairs of sites, we see that the predictions tend to be comparatively slightly worse for core-surface pairs (green curves). Comparing the different contact ranges, we see that the predictions are generally worse for short range contacts, which are more strongly affected by negative design, than for long range contacts, for essentially all types of contacts. Finally, comparing the four types of proteins, we can see that the worst predictions are achieved for short polar proteins, which are less affected by selection for negative design, which is an important component of our model. The peculiarity of short polar proteins could also be related with the fact that many of them occur in multichain complexes, so that the distinction between core and surface sites is not so clear-cut.

b127Weighted Pearson correlation coefficients between the observed and predicted propensities of each amino acid pair, \(Q^\textrm{obs}(a,b)\), and \(Q^{\textrm{pred}}(a,b)\) as a function of the contact range \(|i-j|\), for pairs in contact (first row, A, B, C, D), pairs not in contact (middle row, E, F, G, H) and pairs with indirect contact (bottom row, I, J, K, L). The three curves show three different pairs of sites. The four columns show the different classes of proteins. The dashed line is a reference for the eye at \(r=0.5\)

Mutual Information Increases with the Contact Range for Pairs in Contact and Decreases for Pairs Not in Contact

b129As the result of the different selective pressures that act on native contacts, no-contacts and indirect contacts and on short and long contact ranges, the information content varies for different types of pairs.

b130To illustrate this fact, we computed the observed mutual information for every class of pairs, removing the mutual information obtained from the global propensities \(q^\textrm{glob}\) in order to focus only on the site-specific contributions to mutual information.

b131We show the results in Fig. 6, lumping all proteins together and distinguishing contact types and pairs of sites (surface-surface in Fig. 6A, core-surface in Fig. 6B and core-core in Fig. 6C). As expected, the highest mutual information occurs for pairs in contact, followed by indirect contacts, while pairs that do not form contacts present almost vanishing values (i.e., most of their mutual information arises from the global propensities), except at short range where they are influenced by negative design.

b132As expected, the mutual information of sites in contact tends to increase for longer ranges. In fact, we expect that short range contacts present small mutual information because they are subject to the contrasting forces of positive and negative design.

b133The contrary holds for sites without contacts and with indirect contacts, which present the strongest mutual information at short range, as expected based on negative design.

b134Nevertheless, unexpectedly, pairs of surface sites in contact (Fig. 6A circles) present higher mutual information than pairs of core sites in contact (Fig. 6C circles). Selection for protein folding stability cannot explain this difference, since surface sites contribute less to protein stability. It is intriguing to speculate that this result may be due to selection for function and/or for proper protein interfaces.

b135Mutual information of each class of pairs of sites as a function of the contact range \(|i-j|\), for pairs in contact (cicles), pairs not in contact (squares) and pairs with indirect contact (triangles), distinguishing surface-surface (plot A), core-surface (plot B) and core-core sites (plot C)

Long Proteins are Subject to Stronger Selection

b137Figure 4 suggests that long and hydrophobic proteins are subject to stronger selection for protein folding stability than short and polar proteins. To clarify this aspect, we plotted in Fig. 7A the selection parameter \(\Lambda\) fitted for the four different classes of proteins, considering one global propensity matrix \(q^\textrm{glob}\) for each protein class (Qglob4, black bars), only one \(q^\textrm{glob}\) for all protein class (Qglob1, green bars) and no \(q^\textrm{glob}\) at all (Qglob0, orange bars). We also compared models that consider core and surface sites (NC2, solid bars) with models that consider only one type of site, abolishing the distinction between core and surface sites (NC1, dotted bars).

b138As expected, the more detailed a model is, the lower the selection coefficient is, which indicates that it complies well with the mimimum selection principle. The lowest selection coefficients \(\Lambda\) corresponds to the most accurate global model NC2 Qglob4 (two sites and one global propensity matrix for each protein class).

b139In almost all cases, we find larger selection coefficients for long proteins, which suggests that longer proteins are subject to stronger selection for protein folding stability. This is consistent with our previous conclusion that long proteins are more affected by misfolding problems, and they are subject to stronger negative design (Minning et al. 2013).

b140Nevertheless, when we compare polar and hydrophobic proteins, we find results contrary to our expectations. Misfolding problems are expected to be more serious for hydrophobic proteins, which are subject to stronger negative design (Minning et al. 2013), but the value of \(\Lambda\) is slightly larger for long polar proteins than for long hydrophobic proteins.

b141To further investigate this issue, we plot in Fig. 7C the maximum correlation between propensities and contact energies across no-contact pairs. This correlation allows quantifying negative design without any model. It tends to be larger for hydrophobic than for polar proteins and for long than for short proteins, as expected. Thus, we conclude that our qualitative expectation on negative design is confirmed, but the fitted value of \(\Lambda\) may not be accurate enough to show subtle differences.

b142Similarly, we expect that short proteins experience stronger positive design for stabilizing the native state against unfolding, since the number of contacts per residues increases with chain length due to the surface to volume scaling, and they present fewer native interactions for compensating the loss of conformational entropy upon folding (Bastolla and Demetrius 2005). We quantify positive design for stabilizing the native state against unfolding through the maximum negative correlation between propensities and contact energies across all contact pairs (Fig. 7B). We expect that positive design is stronger for polar proteins, but we find that it is weaker both when assessed through the maximum negative correlation and when assessed through the selection coefficient \(\Lambda\) of short proteins (Fig. 7A), contrary to our expectation.

b143This discrepancy may be due to the fact that many proteins in our data set, in particularly short proteins, tend to occur in protein complexes, and they do not need to be stable on their own. To test this hypothesis, we restricted our data set only to monomeric proteins determined by Xray crystallography, and we tested two thresholds of length, 150 amino acids and 90 amino acids that is often considered the threshold for proteins that fold through two-states folding mechanisms without intermediate states (Barrick 2009). Indeed, for these short lengths the negative design score is very reduced (Fig. 7F orange bars), consistent with the absence of metastable intermediates, and the selection coefficient \(\Lambda\) is larger for short polar proteins than for short hydrophobic proteins (Fig. 7D orange bars), as we expected. However, the positive design score is again larger for short hydrophobic proteins than for short polar proteins (Fig. 7E).

b144Selective parameters \(\Lambda\) fitted for the four different classes of proteins (plots A and D) and maximum correlations between propensities and contact energies for native contacts (B and E) and no contacts (C and F) as a function of different global propensities \(q^\textrm{glob}\) and number of types of sites NC1 and NC2. Plots D,E,F refer to monomeric proteins determined by X ray

Global Amino Acid Correlations are Influenced by the Genetic Code and General Selective Pressures

b146By definition, the global propensities are not site-specific. Therefore, they are not influenced by selection for protein folding stability. They mainly reflect variation of the mutational process and of global (not site-specific) selective forces across the lineages where the proteins evolve. Even if here we do not consider multiple sequence alignments (MSA) or phylogenetic trees from which we can derive evolutionary information, we expect that the global propensities are related with the phylogenetic propensities that can be derived from the MSA and do not depend on site-specific structural or functional constraints.

b147We show in Fig. 8 the ranked values of the global propensities \(q^\textrm{glob}_{a,b}\), focusing on the amino acid pairs that present the most positive (Fig. 8A) and most negative (Fig. 8B) propensities.

b148The first observation is that pairs of the same amino acids have the strongest positive propensities. This is expected since, if the mutational preference for amino acid a is higher in some lineages and lower in other lineages, this variation will create positive self-correlations of the amino acid frequency. More in general, the mutational preferences between amino acids a and b are similar if they are coded by similar codons, predicting positive global propensities in this case. Conversely, if the codons that code two amino acids overlap less than expected by chance, their mutational preferences will be anticorrelated, possibly leading to negative propensities.

b149A similar mechanism holds for variations of the selective preference for any given amino acid. These global preferences may include several selective pressures, which act globally at the level of the whole protein and modulate the mutational propensities:

b150The pressure towards disulfide bridges in short proteins that exist in an oxidative environment (Bastolla and Demetrius 2005). Consistently, Cysteine-Cysteine presents the strongest global propensities \(q^\textrm{glob}_{\textrm{CC}}=0.312\) (Fig. 8).

b151The pressure to avoid biosynthetically costly amino acids, which is particularly relevant in bacterial community, which present hints of cross-feeding where only some community members synthesize some specific amino acids and leak them to the environment (Mustonen and Lässig 2010; Puente-Sánchez et al. 2024), and is lower in organisms that are less limited by resources. Consistently, tryptophane, considered the most biosynthetically costly amino acid, present the second-largest value of the global propensity \(q^\textrm{glob}_{\textrm{WW}}=0.152\), tyrosine presents the 7th highest propensity \(q^\textrm{glob}_{\textrm{YY}}=0.064\) and \(q^\textrm{glob}_{\textrm{WY}}=0.042\). Phenylalanine belongs to the same cluster: \(q^\textrm{glob}_{\textrm{FF}}=0.038\), \(q^\textrm{glob}_{\textrm{FW}}=0.018\) and \(q^\textrm{glob}_{\textrm{FY}}=0.027\).

b152The pressure to avoid codons that frequently mutate to stop codons and can interrupt the translated protein, which mainly applies to the same amino acids discussed above (Tryptophane, Tyrosine and Cysteine), whose codons are close to stop codons.

b153The pressure favoring positively charged amino acids in proteins that interact with nucleic acids that are negatively charged. Consistently, Lysine and Arginine have the third and the sixth-largest global propensity, \(q^\textrm{glob}_{\textrm{KK}}=0.100\) and \(q^\textrm{glob}_{\textrm{RR}}=0.064\). However, \(\textrm{RK}\) is disfavored, \(q^\textrm{glob}_{\textrm{RK}}=-0.058\), possibly because there are specific selective pressures that differentiate Arg and Lys besides their common charge.

b154The pressure towards hydrophobic amino acids in membrane proteins. However, pairs of different hydrophobic residues are not favored, except Methionine and Isoleucine (\(q^\textrm{glob}_{\textrm{MI}}=0.015\)) that are coded by very similar codons.

b155Interestingly, the four polar amino acids Serine, Threonine, Asparagine and Glutamine form a cluster, since they are positively correlated with each other (the highest value is \(q^\textrm{glob}_{\textrm{ST}}=0.030\) and the lowest is \(q^\textrm{glob}_{\textrm{TQ}}=0.003\)). These correlations may also be caused by the similarity of the codons that code for them.

b156Global propensities \(q^\textrm{glob}_{ab }\). Plot A shows the 25 most positive propensities and plot B shows the 25 most negative ones

Discussion

b158In this paper, we introduce the stability constrained model of protein evolution (SCPE) in the pairwise approximation, which generalizes our previous SCPE model that assumes that the protein sites are independent (Arenas et al. 2015). We model the site-specific amino acid distribution of a protein family based on minimal selection for folding stability, requiring that the probability distribution of the amino acids has minimal Kullback-Leibler (KL) divergence from the global (not site-specific) amino acid distribution, while constraining the protein folding stability \(\Delta G\) with respect to both the unfolded state and compact but uncorrectly folded misfolded conformations.

b159Minimum KL divergence coincides with the more familiar maximum entropy principle when amino acids are equiprobable, since maximum entropy is the same as minimum KL divergence from the uniform distribution, but it is more suitable for describing real situations in which the amino acids are not equally distributed. The main assumption of the model is that protein fitness can be approximated with the stability of the native state.

b160Importantly, we do not fit stability constraints from the data, but we predict them from the protein sequence and the experimentally known protein structure with a simple contact model that also takes into account the stability of the misfolded state modelled with the Random Energy Model (Minning et al. 2013). In this way, we predict the parameters of a stability constrained Potts model that describes the amino acid distribution of the protein family in the pairwise approximation (Lapedes et al. 2002; Weigt et al. 2009; Marks et al. 2012; Morcos et al. 2011) only based on global amino acid propensities and site-specific stability constraints. Within the REM, the free energy of the misfolded ensemble depends on groups of four sites, while we approximate it here as a sum of pairwise terms. This approximation may explain in part the reduced accuracy of our predictions for short range pairs, which are dominated by selection on misfolding.

Selective Forces Considered by the SCPE Model

b162We adopt here the small-couplings approximation, i.e. we consider only first order terms in the Potts couplings, which is related with the mutual information approximation originally adopted for inferring the statistical couplings. We obtain an explicit formula that predicts the couplings between amino acid \(a_i\) at site i and amino acid \(a_j\) at site j, and depends on:

b163The global (non-site-specific) couplings \(q^\textrm{glob}_{a_i,a_j}\) between the amino acids \(a_i\) and \(a_j\). In our interpretation, global couplings are the result of variations across organismic lineages of the mutation process and the selective forces that favor or disfavor \(a_i\) and \(a_j\). The global couplings are generally small, and they are strongest when \(a_i=a_j\), since in this case the mutational and selective forces that favor them are perfectly correlated. In particular, we find the strongest self-coupling for Cysteine, which is subject to global selective force for the formation of disulfide bridges. The second one is for Triptophane, which is subject to global selective force for sparing this metabolically costly amino acid. The same applies to Tyrosine and, to a limited extent, Phenilanine, and the three form a cluster of positively correlated amino acids. Another cluster is formed by the polar amino acids Serine, Threonine, Glutamine and Asparagine.

b164Site-specific selective forces due to positive design that stabilize the native state with respect to the unfolded state. They mainly act at pairs of sites that form native contacts. The predicted amino acid propensities of native contacts are approximately equal to minus the interaction free energy of the two amino acids, i.e. amino acid pairs that stabilize the native state are predicted to be more frequently observed at these pairs of sites, with contrasting influence from negative design on short range contacts.

b165Site-specific selective forces due to negative design that destabilize misfolded states. They mainly act at pairs of sites that are at short range along the sequence and do not form native contacts, and they predict amino acid propensities that are approximately equal to the interaction free energy of the two amino acids, i.e. amino acid pairs that destabilize misfolded conformations are more frequently observed at these pairs of sites.

b166We test these predictions on a representative subset of 37108 non-redundant protein chains from the PDB, grouping together 45 classes of pairs of sites with similar properties (native contacts versus no contacts and indirect contacts, 5 different distances along the sequence, core or surface sites). We also consider four types of proteins (long or short, hydrophobic or polar and their combinations).

Test of the SCPE Model Within the SCA Approximation

b168We examine the observed propensities of each class of pairs of sites p, \(Q_{p}(a,b)=P_{p}(a,b)/\left( P_{i}(a)P_{j}(b)\right) -1\), which in the small coupling approximation is proportional to the corresponding coupling.

b169The qualitative dependence of the observed couplings on the contact free energy is for each type of pair in line with what our model predicts (see Fig. 4), and it is modulated by native contacts, the buriedness of the sites, their distance along the sequences, and the length of the proteins as the SCPE model predicts. In particular, the correlations between the couplings and the contact energy are negative for pairs that form native contacts, as expected based on positive design, and they are less negative for short range pairs, as expected based on negative design. They are more negative for core-core sites, as expected based on their importance for protein stability, and they are more negative for longer proteins. For pairs that do not form native contacts, the correlation is positive for short range pairs, as expected based on negative design, and they are stronger for longer proteins, which are expected to be subject to stronger selection against misfolding. This also holds for pairs that form indirect contacts through another residue. On the other hand, the couplings of long range pairs that do not form native contacts are essentially determined by the global couplings, with which they correlate very strongly.

b170The agreement with our predictions is reasonably good, with correlation coefficients r larger than 0.5 for most pairs of sites (see Fig. 3, where the dashed line represents \(r=0.5\)).

b171However, one can see that there are some regimes in which our predictions do not perform very well. In particular, predictions are poor for pairs that form indirect contacts (bottom line in Fig. 3), which are neglected in the small coupling approximation and can be computed only a posteriori. Predictions are also poor for short polar proteins (second column in Fig. 3). This is in part due to the fact that most of these chains (59%) appear in multimeric proteins, so that they may be unstable as monomers. Also, 49% of these chains have been determined by NMR and may be different from average crystal structures that dominate for long proteins. Finally, predictions are poor for short range contacts. There are two likely explanations of this fact: one is the breakdown of the small couplings approximation for short range pairs that do not form native contacts, whose couplings are large (see Fig. 6), and the other one is the approximation of the free energy of the misfolded ensemble predicted through the REM as a sum of pairwise terms, Eq.(), which is necessary for applying the pairwise model but may be not sufficiently accurate.

Unexpected Observations

b173These results suggest that a big part of the couplings in protein evolution originate from the variation of global mutational and selective forces, which act on most pairs of sites, plus couplings due to positive selection for enhancing the stability of the native state, which mainly act at pairs of sites that form native contacts, plus couplings due to negative design for destabilizing pairs that are close along the sequence and do not form native contacts.

b174Perhaps unexpectedly, a larger number of couplings are influenced by negative design against misfolding (positive correlations between observed propensities and contact energies) than those influenced by positive design for stabilizing the native state against unfolding.

b175Another interesting and perhaps unexpected result is that pairs of surface-surface sites in contact possess higher mutual information than pairs of core-core suites in contact (compare circles in Fig. 6A and C) despite the fact that the latters contribute more to protein folding stability. This suggests that surface-surface contacts are subject to other selective pressures besides protein folding stability. Strong candidates are functional movements and interface properties.