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: (1) 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 . It can be expected that the population of sequences follows a Boltzmann distribution [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: (2) (3) 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).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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 (dots).

https://doi.org/10.1371/journal.pcbi.1005584.g002

Therefore, we propose the following model for repeat proteins: (4) 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, J_ij parameters and λ_Id. For example, for pairs of ANK repeats this means 1386 h_i, 959805 J_ij and 34 λ_Id. To reduce the number of free parameters to fit we use a L₁-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 (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.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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

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 Plü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 Plü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 repeat-proteins 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.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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.

https://doi.org/10.1371/journal.pcbi.1005584.g004

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.

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 and J_ij with ), and inter-repeat interactions (J_ij with and , and λ_Id). Then for each repeat we can assign an internal energy and a interaction energy , 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² ≈ 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.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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

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² ≈ 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₁-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) ≠ 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.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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

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 , for 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: (5) As the number of parameters for coupling is large (= 21²L²), we used a regularization L₁ to force to 0 those parameters which are not contributing significantly to the modeled frequencies. Then, we update these parameters by: (6)

Finally, the parameters λ_Id are updated according to: (7)

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.