Inferring repeat-protein energetics from evolutionary information

Natural protein sequences contain a record of their history. A common constraint in a given protein family is the ability to fold to specific structures, and it has been shown possible to infer the main native ensemble by analyzing covariations in extant sequences. Still, many natural proteins that fold into the same structural topology show different stabilization energies, and these are often related to their physiological behavior. We propose a description for the energetic variation given by sequence modifications in repeat proteins, systems for which the overall problem is simplified by their inherent symmetry. We explicitly account for single amino acid and pair-wise interactions and treat higher order correlations with a single term. We show that the resulting evolutionary field can be interpreted with structural detail. We trace the variations in the energetic scores of natural proteins and relate them to their experimental characterization. The resulting energetic evolutionary field allows the prediction of the folding free energy change for several mutants, and can be used to generate synthetic sequences that are statistically indistinguishable from the natural counterparts.


Introduction
Repeat proteins are composed of tandem repetitions of similar structural motifs of about 20 to 40 amino acids.Under appropriate conditions, these polymers fold into elongated, non-globular structures (Fig 1).It is apparent that the overall architecture is stabilized mainly by short range interactions, in contrast to most globular protein domains that usually adopt very intricate topologies [1].In their natural context, repeat proteins are frequently found mediating protein-protein interactions, with a specificity rivaling that of antibodies [2][3][4].Given their structural simplicity and potential technological applications, repeat-proteins are a prime target for protein design, with very successful examples for a variety of topologies [5][6][7].Most of the current design strategies target the creation of rigid native structures with desired folds that, although beautiful, often lose biological functionality [8].It is becoming clear that the population of 'excited states' is crucial for protein function [9], and thus tackling energetic inhomogeneities in protein structures may be crucial for understanding how biological activities emerge [10].The challenge thus relies in finding an appropriate description for the 'energy' of each system, a daunting task for large molecular objects such as natural proteins.
In principle, the natural variations observed for proteins of the same family must contain information about the sequence-structure mapping.A simple model that just takes into account the frequency of each amino acid in each position is insufficient to capture collective effects, yet, for some architectures it is surprisingly good for the synthesis of non-natural repeat-proteins by 'consensus' design [11][12][13][14].It is apparent that in the case of repeat proteins the local signals play inordinately large roles in the energy distribution, just as expected from their topology [15] and hence, small heterogeneities can be propagated from the local repeat units to higher orders affecting the overall structure and dynamics [16,17].Thus, collective effects may be approximated as small perturbations to local potentials, simplifying the energetic description of complex natural systems [18].
In the last years new methods to analyze correlated mutations across a family of proteins have arisen (mfDCA [19], plmDCA [20,21], Gremlin [22] to name a few).The main hypothesis behind these methods is that biochemical changes produced by a point mutation should be compensated by other mutations (along evolutionary timescales) to maintain protein viability or function.These methods can also be used to disentangle relevant direct correlations from indirect ones.They are very successful at predicting spatial contacts and interactions for many protein topologies [23][24][25][26][27]. Nevertheless, these methods do not take into account the chemical nature of the amino acids, which can be codifying inhomogeneities in the energetic distribution that are crucial for the activity of repeat-proteins [28,29].On this basis, different approaches have been proposed recently to include chemical details in the correlation analyses [30], trying to predict folding stability [31], conformational heterogeneity [23,32,33], mutational effect in the interaction in two-component signaling proteins [27] or the global effect on antibiotic resistance from sequences of β-lactamases [34,35].As many other tools, these were optimized to perform well on globular proteins, and their application to repeat proteins is not straightforward.Besides the point-mutation mechanism, repeat proteins are believed to evolve via duplication and rearrangement of repeats [36], resulting in an inherent symmetry which usually confounds sequence analyses [17].Making use of this symmetry, we have previously proposed a specific version of mfDCA and plmDCA for repeat proteins [37].In this work we develop an alternative 'evolutionary field', able to disentangle biases generated by the repetitive nature of these proteins and which explicitly includes the information of the amino acids that compose a protein.We show that it is able to reflect biochemical properties of the analyzed proteins.We take advantage of the elongated and repetitive structure of these proteins (Fig 1) to extract as much information as possible from the data, and apply the general ideas on three specific families, ankyrin repeats (ANK), leucine-rich repeats (LRR) and tetratricopeptide-like repeats (TPR).

Evolutionary energy for repeat proteins
To study the co-occurrence of mutations in a sequence alignment of a particular protein family, [39] proposed a Hamiltonian or energy expression which resembles a Potts model: where the set of {h i (a i )} parameters, one for each amino acid in each position, accounts for a local propensity of having a specific residue on a particular site of the protein, and the set of {J ij (a i , b j )} indicates the strength of the 'evolutionary' interaction between each possible amino acid in every pair of positions along the protein.There are q = 21 possible values of a i and b j , one for each amino acid and one for the gaps included on the multiple sequence alignments.This expression is evaluated on a particular sequence on an alignment, and the summations go over the L columns of the alignment.A sequence is more favorable or more energetic if it gets lower values of EðsÞ.It can be expected that the population of sequences follows a Boltzmann distribution PðsÞ ¼ 1 Z e À EðsÞ [40].The parameters are thus fitted to reproduce the frequencies of occurrence of each amino acid in each position (f i (a i )) and the joint frequencies of amino acids (f ij (a i , b j )) in an alignment of natural sequences used as input: Nevertheless, for repeat proteins there is another feature we want to capture with an evolutionary energy: the high identity of amino acids constituting consecutive repeats, arisen by the repetitiveness of these families and probably a signature of their evolutionary mechanisms (Fig 2) . Therefore, we propose the following model for repeat proteins: This expression is designed to be applied in sequences constituted by two repeats.λ Id is a parameter that aims at reproducing the probabilities of the percentage of identity (%Id) between consecutive repeats in natural proteins (p id ).Basically, it accounts for higher order correlations not captured by the pairwise terms.For a given sequence we calculate the %Id of the adjacent repeats and sum the parameter λ Id corresponding to that %Id value.When the correct parameters are obtained, this equation can be used to produce an ensemble of sequences consistent with the constraints (f i (a i ), f ij (a i , b j ) and p id ).We work with pairs of repeats as it is the minimum unit that includes the interaction between repeats and the possibility of measure sequence identity between consecutive repeats.In the following section we will show the convergence of the method and the relevant information that can be obtained from it.For further details about the procedure to assign values to the parameters, please refer to Methods section.

Evolutionary energy reproduces ensembles of sequences with natural frequencies and repeat protein characteristics
We construct an alignment of pairs of repeats for each family: ANK (PFAM id PF00023, and final alignment of 20513 sequences of L = 66 residues each), TPR (PFAM id PF00515, and final alignment of 10020 sequences of L = 68 residues each) and LRR (PFAM id PF13516, and final alignment of 18839 sequences of L = 48 residues each).See Methods for further details of construction.We measure f i (a i ), f ij (a i , b j ) and p id .Using a gradient descent procedure we obtain a set of parameters in eq 4 which are able to reproduce f i (a i ), f ij (a i , b j ) and p id .In principle, the number of parameters is large: Lq h i parameters, ðLqÞ 2 À Lq 2 J ij parameters and L 2 þ 1 λ Id .For example, for pairs of ANK repeats this means 1386 h i , 959805 J ij and 34 λ Id .To reduce the The proposed model fits the frequencies of amino acids and natural repeat identities p id .On A, we compare marginal frequencies f i (a i ) (red) and joint frequencies f ij (a i , b j ) (black) on the natural ensemble of sequences (x-axis) and on the set of sequences generated by the model (y-axis).On B, we calculate the distribution of identity between repeats p id for consecutive repeats (solid line), and for natural repeats which are not consecutive, i.e. they are not next to each other in the primary structure (dot lines).Consecutive repeats present a population with high identity between repeats that any pairs of repeats do not show.We compare the distribution produced by the model p model id (dots).https://doi.org/10.1371/journal.pcbi.1005584.g002number of free parameters to fit we use a L 1 -regularization which fixes to zero those parameters which do not contribute significantly to fit the frequencies.This regularization allows us to set to exactly zero between 85 and 91% of the J ij parameters which, when they are free to vary, only reach small values (S3 Fig) .We bound the maximum error permitted in the frequency estimations to 0.02.Refer to Methods for more details.
In the three families studied, the parameters obtained allow us to generate ensembles of sequences which reproduce natural f i (a i ), f ij (a i , b j ) and p id (Fig 2A).Notice that most frequencies are fitted with an error of an order of magnitude lower than the maximum bound imposed (S2 Fig) .The p id distributions are also very well reproduced (Fig 2B).Not only the general shape, but also the populated long tail for highly similar repeats.It is not possible to obtain the same distribution only by fitting amino acid frequencies f i (a i ) and f ij (a i , b j ), it is mandatory to explicitly include the p id by including the parameters λ Id (S1 Fig), suggesting that higher order correlations must be accounted for describing these systems.

Evolutionary energy distinguishes between proteins on a given family and other polypeptides
Once the set of parameters {h i (a i ), J ij (a i , b j ), λ Id } is obtained, it can be used to score any sequence of L amino acids via eq 4. In this section we test if this measure is capable of distinguishing polypeptides that fold in a three dimensional structure similar to members of the repeat protein family from those that do not.
We calculate the distribution of energies of different sets of sequences ( For a positive control we evaluate designed proteins which have been experimentally synthesized.For the ANK family, we consider the library of repeat sequences built by Plu ¨ckthun's laboratory [13] (green lines, Fig 3A).This library was constructed by fixing on each repeat 26 positions out of 33 to the most frequent residue in the multiple sequence alignment.This resulted in a set of sequences that have small variations with respect to the ANK consensus (the sequence with the most frequent amino acid in each position).In our expression, they score a very low energy distribution, overlapping with the most negative tail of the distribution of natural sequences.It is notable that consensus designed ANK have been shown experimentally to be extremely stable.For the TPR family, consensus designed was done by Regan's laboratory [11,12].All pairs of repeats synthesized have the same amino acid sequence, and it's energy score is indicated by a green full square in Fig 3B .Again, the designed sequence matches values at the most left side of the energy distribution of natural sequences, and coincidentally reports high folding stability.From it, other variants with few point mutations to improve binding to a specific ligand have been synthesized.As shown in empty green squares [41] and diamonds [42] in Fig 3B, they have higher energy, but still in the left most side of natural sequences distribution.Recently, a different design strategy was done [43].Based on a non-repetitive protein, but similar to TPR fold, they put togheter various repetitions of the fold, using TPR loops to link them.They obtained a three-repeats protein whose pair of repeats energy are represented on triangles on Fig 3B.This time, they match natural sequences distribution in higher values.
Finally, for the LRR family we contrast with the library of proteins designed by Plu ¨ckthun's group based on the consensus sequence [14].The repetition they considered has 57 amino acids, which includes two types of repeats, one of 28 residues and the other one of 29.As the repeat we are using for LRR is 24 residues long, we aligned both definitions and evaluated the library removing the amino acids not matching our definition.Again, their scores form a narrow distribution, but this time it is not placed on the most favorable side of the natural sequences distribution (Fig 3C).Coincidentally, selected species studied do not show such a high folding stability as the ANK library did.
With these parameters, we are able to generate an ensemble of sequences which are in agreement with the constraints used, via a Monte Carlo simulation (see Methods).The distribution of energies of these simulated sequences matches the natural sequences energies distribution with remarkable accuracy.Moreover, we randomly choose 100 sequences from the natural ensemble and 100 sequences from the simulated one, perform a Smith-Waterman pairwise alignment all against all, calculate the pair similarity using BLOSUM62 matrix and used it as a distance method to plot a dendogram of the sequences (S4 Fig) .Both species appear interspersed, showing that it is not possible to distinguish a natural sequence from a constructed one.Also, we tested familiarity to the ANK family as defined in [44] and found overlapping distributions for both species (S5 Fig).Therefore, simulated sequences represent possible variants to natural repeats.The wide distribution of natural proteins suggests that it should be possible to engineer sequences with more variable repeats, more dissimilar among neighbors and to the consensus than the ones published up to date.

Low evolutionary energy sequences have similar repeats
Are there any invariant properties shared by low energy sequences?Given that repeatproteins may evolve by other mechanisms besides point substitutions, we analyze if low energy sequences are constituted by highly similar repeats and if they are close to consensus sequences.
On Fig 4A we show the relation between the %Id between the repeats and the energy of the sequence.It is evident that low energy sequences are constructed by pairs of highly similar repeats.This could be a transitive effect: if low energy sequences are very similar to the consensus sequence, and the consensus sequence is formed by two identical repeats, we would be seeing that more similarity between repeats causes lower energies.We can see that it is not the case (Fig 4B).We plot the %Id to the consensus against the energy of each sequence.The consensus was calculated with the most frequent amino acid in each position on sequences used as input.We can see that there is no evident correlation between the energy and the similarity to the consensus.Thus, low energy sequences that differ from the consensus one may be constructed.Also, there are no sequences which get a high %Id to the consensus.We conclude that there are different repeats which have low energies within a protein family, and not only the consensus sequence.

Evolutionary energy and folding stability change upon point mutations
Consensus designed ANK proteins are very stable upon chemical and thermal denaturation [13], and, as shown in Fig 3 also score a very low evolutionary energy according to eq 4. Can we quantify the relationship between the stability and the evolutionary energy?
A potential test can be performed by comparing to experiments in which the effect of point mutations was evaluated.These incorporate one, two or three point mutations in natural proteins, and characterize the unfolding free energy ΔG of the wildtype and the mutated variant.A higher ΔG reports a more stable protein.We compare the change in the ΔG between the mutated and the wildtype protein (ΔΔG), and the difference of energy for their sequences according to eq 4. the identity between the repeats that constitute the sequence.Even though the deviation is large, most stable sequences tend to have more similar repeats.On B, we plot the energies of simulated sequences vs the identity to the consensus of the family.In all cases, the identity to the consensus is small and uncorrelated to the energy, indicating that sequences which differ significantly from the consensus can be stable variants of the family.https://doi.org/10.1371/journal.pcbi.1005584.g004 Although the energy expression is learned for pairs of repeats, we can easily extend it to an array of repeats making use of the elongated structure of repeat proteins in which only adjacent repeats interact.From our expression we have parameters assigned to intra-repeat positions (h i with i ¼ 1 . . .L 2 and J ij with i; j ¼ 1 . . .L 2 ), and inter-repeat interactions (J ij with i ¼ 1 . . .L 2 and j ¼ L 2 þ 1 . . .L, and λ Id ).Then for each repeat we can assign an internal energy Þ and a interaction energy Id , which of course depends on the amino acids constituting each repeat.
On Fig 5A, we show the comparison between ΔΔG and the evolutionary energy calculated using Eq 4, done for three different ANK proteins: IκBα [45,46], Notch [47] and p16 [48].It should be noted that different experimental techniques return different values for ΔG for the same protein, non overlapping within experimental error, pointing that other factors contribute to the experimental quantification of ΔΔG.A linear fit returns R 2 % 0.61.From 152 mutations we analyzed, 124 (82%) are predicted favorable when the mutation stabilized the folding of the structure, and unfavorable when they have also been measured to destabilize.The predictions that deviated the most are mutations in Notch from Serine to Proline, which is a structural disruptor, and were not considered in the linear fit.A comparison against FoldX [49] predictions can be found on S6 Fig.
On Fig 5B, we show reported mutations on pp32 [50], a protein belonging to LRR family.Again, measurements with different methods report different values of ΔΔG.The linear fit returns a poor R 2 % 0.21, but 30 (75%) mutations are both predicted and reported unstabilizing.
A similar comparison was performed by [31] for small globular proteins with an expression related to Eq 1.To reduce the number of interaction parameters J ij (a i , b j ) they explicitly used structural information and set to zero all interactions between positions which are not in contact in the native structure.In contrast, we use a L 1 -regularization to fix to zero those parameters which do not contribute significantly to the fitting process and obtain J ij (a i , b j ) = 0 and J ij (a i , b j ) 6 ¼ 0 in all pairs of positions, regardless they are supposed to be in contact or not in the 3D structure.

Interaction parameters are related to the structure and the sequence symmetry
Are the obtained parameters related to structural properties of these proteins?Local fields, h i (a i ), should account for the local propensity of each amino acid in each position, and therefore are expected to be related to f i (a i ).Fig 6A shows that the inferred h i (a i ) parameters are different from the initial condition ln(f i (a i )) for the ANK family; that is, the values obtained for the parameters that account for higher order correlations are relevant.In red we highlight the points related to the consensus amino acid in each position.All of these residues have a strong local field associated to them, justifying why the construction of sequences with these amino acids results in foldable proteins.We also show a contact map of two ANK repeats (PDB id: 1N0R) on Fig 6B : gray background indicates that the two positions given by x and y axis are in contact in the native structure, and white that they are not.On the upper triangle of the figure and in blue crosses, we mark the positions involved in the highest J ij parameters, i.e. those which imply higher coupling.A darker blue indicates that there are more J ij (more combinations of amino acids) between those positions.Most of the highest J ij match a pair of positions in contact in the 3D structure, or two which correspond to the same residue in the adjacent repeat patterns, i.e. i-th position in the first repeat and position j = i+33 in the second repeat.In red crosses we show the lowest J ij , that mark a negative constraint.Again, a darker red means that there are more J ij with low values between those positions.It is apparent that these also involve mostly residues in contact, but shows that other regions are responsible for negative design.

Discussion
We propose a statistical model to account for fine details of the energy distribution in families of repeat proteins using only the sequences of amino acids.The model consists of a specialization of a Potts model to account for the local and pair-wise interactions and an extra term that includes higher order correlations, accounting for the similarity between consecutive repeats.The model is constrained by evolutionary characteristics of the families of proteins: we measure the frequencies of amino acids, co-occurrence of amino acids and the identity between repeats in extant natural proteins.To statistically define these quantities it is necessary to have a large set of sequences, which we showed are currently available for several repeat-protein families [37].No information about the native folded conformation is required.The computation of the evolutionary energy field is computationally demanding, mostly due to long times spent in rigorous Monte Carlo simulations, but once the fitting is done the parameters can be used to score individual sequences fast and easily.
We studied three popular repeat protein families: ANK, TPR and LRR.After pre-processing of the alignments, we had enough sequences (% 20500, 10000 and 18800 respectively) to fit the model to pairs of repeats of each family.We scored the evolutionary energy of all natural sequences in PFAM, and it allowed us to clearly distinguish between natural proteins and random sequences of amino acids: the first have energy values < -50 and show a large spread while all random sequences have energy values % 0. We evaluated designed repeat proteins which have been shown to fold and found that they score within the natural sequences distribution of energies.For the ANK and TPR families, these designed proteins have been shown to be highly stable upon thermal and chemical denaturation and, coincidentally, they are located at the most favorable side of the energy distribution of natural proteins, suggesting that the evolutionary energy score can be related to folding stability.
The energetic model can be used in Monte Carlo simulations to generate sequences that agree with the natural constraints of a given protein family.This ensemble of simulated sequences matches the amino acid frequencies, the identity between repeats and also the energy distribution of natural proteins.We found this set of simulated sequences is statistically indistinguishable from natural counterparts.Thus, the proposed model can be used as a tool to design repeat-protein sequences that have all the natural characteristics evaluated to date.Repeat proteins bind to other polypeptides and are candidates for specific binder scaffolds.Designed repeat proteins have been successfully synthesized and adapted to biomedical applications.Nevertheless, consensus design limits the possible variants as only a small proportion of residues are free to vary.Furthermore, they are extremely stable.Including coupling information can wide the possible sequences that can be studied, and could lead to more malleability of the designed molecules.Moreover, the stability change upon single point mutation can be well predicted by the model using just sequence information.For the ANK family, evolutionary energy variations correlate with the experimental values with an R 2 %0.6.This improves FoldX [49] performance, which additionally requires a reference structure.Moreover, from the 152 experiments analyzed, the 82% predicts the direction of the stability change upon a point mutation.For the LRR family, the correlation is considerably lower, but 75% of the mutations are both predicted and reported in the available bibliography as destabilizing.For both the simulated sequences and for natural counterparts, we found that the similarity between consecutive repeats correlates with lower energy values, and that these are not necessarily similar to the consensus sequence of the family, pointing out that duplication of stretches of sequences may well be an important factor in the evolution of these systems [51].
The existence of a simple and reliable energy function to score the 'evolutionary energy' of repeat-proteins can be used to trace the biological forces that acted upon their history, and to explore to which extent these conflict with the physical necessities of the systems [52].Mapping the energy inhomogeneities along the repeat-arrays may allow us to infer the population of excited states in these proteins, many of which have been related to their physiological mechanisms.

Sequence alignments
Multiple sequence alignments of repeats were obtained from PFAM 27.0 [53].The aligned sequences usually have misdetected initial and final residues.The amino acids at the ends of the repeat-detection do occur in the polypeptide chains (they are not actual deletions) and incorporating them improves the statistics of the real sequences.We completed these positions with the amino acids present on the actual proteins using the provided headers on the alignment and crossing information with UniProt database [54].This leads to a reduction on the number of gaps in our alignments, which usually derives into noisy predictions in correlation analyses [31].After, we created the alignment of pairs of repeats, joining sequences of repeats which are consecutive in a natural protein.Finally, we removed insertions from the alignments by deleting positions which have gaps in more than 80% of the sequences in the alignment.

Frequency calculations
Our model fits the occurrence of amino acids in every position, which we call the marginal frequency of residue a i at position i of the alignment and denote f i (a i ), and the joint occurrence of two amino acids a i and b j simultaneously at two different positions of the alignment, f ij (a i , b j ).To avoid biases by the overrepresentation of some proteins in the database, we used CD-HIT [55] to cluster sequences at 90% of identity and chose a representative sequence from each cluster.Finally, we computed by counting the f i (a i ) and f ij (a i , b j ), and divided by the total number of sequences.

p id calculations
From the same alignment explained in Frequencies calculations, for a sequence which has L residues constituting two consecutive repeats, the %Id between the repeats is the number of amino acids in positions i and i þ L 2 , for i ¼ 1 . . .L 2 which are exactly the same.Gaps are treated as an amino acid.Once we have the values for all sequences in an alignment, we define p id as the proportion of sequences within the alignment with the same %Id between repeats.

Construction of an ensemble of sequences in agreement with a energy equation
Given a set of parameters h i , J ij , λ Id and Eq 4, we use a Monte Carlo procedure and the Metropolis criterion to generate an ensemble of N sequences of length L each.We initiate with a random string of L residues.At each step, we produce a point mutation in any position.If this mutation is favorable, i.e. the energy is lower than that of the original sequence, we accept the mutation.If not, we accept the mutation with a probability of e −ΔE , where ΔE is the difference of energy between the original and the mutated sequence.When accepted, the mutated sequence is used as the original one for next step.We add one sequence to our final ensemble every t steps (we used t = 1000).

Learning the parameters for the model
Our model is proposed to reproduce f i (a i ), f ij (a i , b j ) and p id from the alignment of natural sequences.To learn the set of parameters h i , J ij , λ Id which reproduce them, we used a gradient descent procedure.In each step, an ensemble of N = 80000 sample sequences was produced via Monte Carlo using as energy the expression 4 and the trial parameters.We measured its marginal, joint frequencies and p id and updated the local parameters according to: As the number of parameters for coupling is large (= 21 2 L 2 ), we used a regularization L 1 to force to 0 those parameters which are not contributing significantly to the modeled frequencies.Then, we update these parameters by: Finally, the parameters λ Id are updated according to: We iterated until the maximum difference between the predicted frequencies and the natural sequences was below 0.02.This value was chosen according to the robustness of the frequencies estimations on the available data.We calculated the frequencies on half of the available sequences and compared the results to the frequencies counts on all the available sequences.The largest differences were slightly below 0.02.We believe that this maximum error thus reflects the actual error in the data and it is not reasonable to ask the model for more accuracy than that of the data itself.The code was written in C++ and is available at GitHub: https://github.com/proteinphysiologylab/2017_Espadaetal.h i (a i ) and the marginal frequencies f i (a i ).At center, contact map (grey indicates position in contact, white not in contact on the native structure).On blue, pairs of positions involved in highest J ij (a i , b j ), red lowest J ij (a i , b j ).At right, pairs of amino acids involved in the highest J ij (a i , b j ) parameters (on blue) and in the lowest J ij (a i , b j ) (on red).ANK at the top.TPR at the center.LRR at the bottom.(PDF)

Fig 2 .
Fig 2.The proposed model fits the frequencies of amino acids and natural repeat identities p id .On A, we compare marginal frequencies f i (a i ) (red) and joint frequencies f ij (a i , b j ) (black) on the natural ensemble of sequences (x-axis) and on the set of sequences generated by the model (y-axis).On B, we calculate the distribution of identity between repeats p id for consecutive repeats (solid line), and for natural repeats which are not consecutive, i.e. they are not next to each other in the primary structure (dot lines).Consecutive repeats present a population with high identity between repeats that any pairs of repeats do not show.We compare the distribution produced by the model p model

Fig 3 )
. The ensembles of natural sequences of each protein family used to learn the parameters have a unimodal distribution of energies centered around -100 (Fig 3, red lines).These distributions are clearly differentiated from the energies of random chains of residues (Fig 3, yellow lines), which constitute a basic negative control for our model.

Fig 3 .
Fig 3. Energy score distribution for different ensembles of sequences for ANK (A panel), TPR (B) and LRR (C).Red lines, natural sequences used to train the model on eq 4. Blue lines, sequences simulated by Monte Carlo under expression 4. In the three families, it overlaps with natural sequences, suggesting that simulated sequences imitate the natural ensemble.Yellow lines, strings of random amino acids used as negative control.They show that the energy distinguishes between polypeptides belonging to a protein family and other strings of amino acids.Green lines, squares, diamonds and triangles, energies for designed proteins.https://doi.org/10.1371/journal.pcbi.1005584.g003

Fig 4 .
Fig 4. Most favorable simulated sequences have very similar repeats, yet they are different to the consensus repeat.On A, we plot the energy vs.the identity between the repeats that constitute the sequence.Even though the deviation is large, most stable sequences tend to have more similar repeats.On B, we plot the energies of simulated sequences vs the identity to the consensus of the family.In all cases, the identity to the consensus is small and uncorrelated to the energy, indicating that sequences which differ significantly from the consensus can be stable variants of the family.

Fig 5 .
Fig 5. Variation of energy score as a predictor of the folding stability upon point mutations.We compare difference in unfolding ΔG between a wildtype protein and a mutated variant (x-axis) and the change in energy according to Eq 4. Error bars indicate the experimental standard deviation.On A, for proteins belonging to ANK family, and on B for LRR.https://doi.org/10.1371/journal.pcbi.1005584.g005

Fig 6 .
Fig 6.For the ANK family, on panel A we compare the parameters h i (a i ) to the marginal frequencies.The site-independent model (and initial condition) states that h i (a i ) = ln(f i (a i )).For the final model, this relation is tuned by the higher order correlations.On red, the parameters associated with the most common amino acid in each position are highlighted.On panel B, we compare the contact map of a pair of repeats of 1N0R (gray shadow) and the highest (blue) and lowest (red) J ij (a i , b j ) parameters.The color scale indicates how many parameters involves the two positions (due to different sets of amino acids).Most extreme values fall into residues in contact or in the equivalent position of a repeat.https://doi.org/10.1371/journal.pcbi.1005584.g006

S4Fig.
Dendrogram based on pair-wise similarity.Natural and simulated sequences are indistinguishable from pairwise similarity.(PDF) S5 Fig. A) Logos for the MSA of natural pairs of ANK repeats (top), of simulated pairs of ANK repeats (center) and low energy pairs of simulated ANK repeats (bottom).B) Distribution of familiarity, as defined in Turjanski et al (2016) for the same sets of sequences.Simulated sequences reproduce the distribution of natural proteins and are indistinguishable.(PNG) S6 Fig.For the ANK family, we compared the results of our models on folding stability with FoldX, a popular tool used to predict stability changes upon mutation.As FoldX requires a reference structure, some of the constructions tested in the protein Notch cannot be analyzed, so we excluded them from our predictions for a fair comparison.It can be seen that, overall, our model (right panel) is a better predictor than foldX (left panel) of stability changes upon mutation, and has the advantage of using just the sequence information.(PNG) S7 Fig.At left, comparison between the local field parameters