Size and structure of the sequence space of repeat proteins

Abstract

The coding space of protein sequences is shaped by evolutionary constraints set by requirements of function and stability. We show that the coding space of a given protein family—the total number of sequences in that family—can be estimated using models of maximum entropy trained on multiple sequence alignments of naturally occuring amino acid sequences. We analyzed and calculated the size of three abundant repeat proteins families, whose members are large proteins made of many repetitions of conserved portions of ∼30 amino acids. While amino acid conservation at each position of the alignment explains most of the reduction of diversity relative to completely random sequences, we found that correlations between amino acid usage at different positions significantly impact that diversity. We quantified the impact of different types of correlations, functional and evolutionary, on sequence diversity. Analysis of the detailed structure of the coding space of the families revealed a rugged landscape, with many local energy minima of varying sizes with a hierarchical structure, reminiscent of fustrated energy landscapes of spin glass in physics. This clustered structure indicates a multiplicity of subtypes within each family, and suggests new strategies for protein design.

About the Authors:

Jacopo Marchi

Roles Formal analysis, Investigation, Methodology, Validation, Visualization, Writing – original draft

Affiliation: Laboratoire de physique de l’École normale supérieure (PSL University), CNRS, Sorbonne Université, and Université de Paris, 75005 Paris, France

Ezequiel A. Galpern

Roles Data curation, Formal analysis, Methodology

Affiliations Protein Physiology Lab, Universidad de Buenos Aires, Facultad de Ciencias Exactas y Naturales, Departamento de Química Biológica, Buenos Aires, Argentina, CONICET - Universidad de Buenos Aires, Instituto de Química Biológica de la Facultad de Ciencias Exactas y Naturales (IQUIBICEN), Buenos Aires, Argentina

ORCID logo http://orcid.org/0000-0001-9516-3985

Rocio Espada

Roles Data curation

Affiliation: Laboratoire Gulliver, Ecole supérieure de physique et chimie industrielles (PSL University) and CNRS, 75005, Paris, France

Diego U. Ferreiro

Roles Conceptualization, Data curation, Funding acquisition, Investigation, Writing – review & editing

Affiliations Protein Physiology Lab, Universidad de Buenos Aires, Facultad de Ciencias Exactas y Naturales, Departamento de Química Biológica, Buenos Aires, Argentina, CONICET - Universidad de Buenos Aires, Instituto de Química Biológica de la Facultad de Ciencias Exactas y Naturales (IQUIBICEN), Buenos Aires, Argentina

Aleksandra M. Walczak

Contributed equally to this work with: Aleksandra M. Walczak, Thierry Mora

Roles Conceptualization, Formal analysis, Funding acquisition, Investigation, Project administration, Supervision, Writing – original draft, Writing – review & editing

* E-mail: [email protected] (AMW); [email protected] (TM)

Affiliation: Laboratoire de physique de l’École normale supérieure (PSL University), CNRS, Sorbonne Université, and Université de Paris, 75005 Paris, France

Thierry Mora

Contributed equally to this work with: Aleksandra M. Walczak, Thierry Mora

Roles Conceptualization, Formal analysis, Investigation, Methodology, Project administration, Supervision, Writing – original draft, Writing – review & editing

* E-mail: [email protected] (AMW); [email protected] (TM)

Affiliation: Laboratoire de physique de l’École normale supérieure (PSL University), CNRS, Sorbonne Université, and Université de Paris, 75005 Paris, France

ORCID logo http://orcid.org/0000-0002-5456-9361

Introduction

Natural proteins contain a record of their evolutionary history, as selective pressure constrains their amino-acid sequences to perform certain functions. However, if we take all proteins found in nature, their sequence appears to be random, without any apparent rules that distinguish their sequences from arbitrary polypeptides. Nonetheless, the volume of sequence space taken up by existing proteins is very small compared to all possible polypeptide strings of a given length [1], even more so when specializing to a given structure [2]. Clearly, not all variants are equally likely to survive [3–5]. To better understand the structure of the space of natural proteins, it is useful to group them into families of proteins with similar fold, function, and sequence, believed to be under a common selective pressure. Assuming that the ensemble of protein families is equilibrated, there should exist a relationship between the conserved features of their amino acid sequences and their function. This relation can be extracted by examining statistics of amino-acid composition, starting with single sites in multiple alignments (as provided by e.g. PFAM [6, 7]). More interesting information can be extracted from covariation of amino acid usages at pairs of positions [8–10] or using machine-learning techniques [11]. Models of protein sequences based of pairwise covariations have been shown to successfully predict pair-wise amino-acid contacts in three dimensional structures [12–17], aid protein folding algorithms [18, 19], and predict the effect of point mutations [17, 20–22]. However, little is known on how these identified amino-acid constraints affect the global size, shape and structure of the sequence space. Accounting for these questions is a first step towards drawing out the possible and the realized evolutionary trajectories of protein sequences [23, 24].

We use tools and concepts from the statistical mechanics of disordered systems to study collective, protein-wide effects and to understand how evolutionary constraints shape the landscape of protein families. We go beyond previous work which focused on local effects—pairwise contacts between residues, effect of single amino-acid mutations—to ask how amino-acid conservation and covariation restrict and shape the landscape of sequences in a family. Specifically, we characterize the size of the ensemble, defined as the effective number of sequences of a familiy, as well as its detailed structure: is it made of one block or divided into clusters of “basins”? These are intrinsically collective properties that can not be assessed locally.

Repeat proteins are excellent systems in which to quantify these collective effects, as they combine both local and global interactions. Repeat proteins are found as domains or subdomains in a very large number of functionally important proteins, in particular signaling proteins (e.g. NF-κB, p16, Notch [25]). Usually they are composed of tandem repetitions of ∼30 amino-acids that fold into elongated architectures. Repeat proteins have been divided into different families based on their structural similarity. Here we consider three abundant repeat protein families: ankyrin repeats (ANK), tetratricopeptide repeats (TPR), leucine-rich repeat (LRR) that fold into repetitive structures (see Fig 1). In addition to interactions between residues within one repeat, repeat protein evolution is constrained by inter-repeat interactions, which lead to the characteristic accordeon-like folds. Through these separable types of constraints, as well as the possibility of intra- and inter-familly comparisons, repeat proteins are perfect candidates to ask questions about the origins and the effects of the constraints that globally shape the sequences.

[Figure omitted. See PDF.]

Fig 1. Repeat proteins fold into characteristic accordeon-like folds.

Example structures of three protein families are shown, ankyrin repeats (ANK), tetratricopeptide repeats (TPR), leucine-rich repeat (LRR), with the repeating unit highlighted in magenta. All show regular folding patterns with defined contacts in and between repeats.

https://doi.org/10.1371/journal.pcbi.1007282.g001

A recent study [26] addressed the question of the total number of sequences within a given protein family, focusing on ten single-domain families. They took a similar thermodynamic approach to the one followed here, but had to estimate experimentally the free energy threshold ΔG below which the sequences would fold properly. Here we overcome this limitation by forgoing this threshold entirely. Instead we determine the sequence entropy directly, which is argued to be equivalent to using a threshold free energy by virtue of the equivalence of ensembles. We precisely quantify the sequence entropy of three repeat-protein families for which detailed evolutionary energetic fields are known [27]. We explore the properties of the evolutionary landscape shaped by the amino-acid frequency constraints and correlations. We ask whether the energy landscape, defined in sequence space of repeat proteins, is made of a single basin, or rather of a multitude of basins connected by ridges and passes, called “metastable states”, as would be expected from spin-glass theory. Using the specific example of repeat proteins makes it possible to analyze the source of the potential landscape ruggedness, and use it to identify which repeat-protein families can be well separated into subfamilies. The rich metastable state structure that we find demonstrates the importance of interactions in shaping the protein family ensemble.

Results

Statistical models of repeat-protein families

We start by building statistical models for the three repeat protein families presented in Fig 1 (ANK, TPR, LRR). These models give the probability P(σ) to find in the family of interest a particular sequence σ = (σ1, …, σ2L) for two consecutive repeats of size L. The model is designed to be as random as possible, while agreeing with key statistics of variation and co-variation in a multiple sequence alignment of the protein family. Specifically, P(σ) is obtained as the distribution of maximum entropy [28] which has the same amino-acid frequencies at each position as in the alignment, as well as the same joint frequencies of amino acid usage in each pair of positions. Additionally, repeat proteins share many amino acids between consecutive repeats, both due to sharing a common ancestor and to evolutionary selection acting on the protein. To account for this special property of repeat proteins, we require that the model reproduces the distribution of overlaps between consecutive repeats. Using the technique of Lagrange multipliers, the distribution can be shown to take the form [17]:(1)with(2)where hi(σ), Jij(σi, σj), and {λID}, ID = 0, 1, …, L, are adjustable Lagrange multipliers that are fit to the data to reproduce the experimentally observed site-dependent amino-acid frequencies fi(σi), joint probabilities between two positions, fij(σi, σj), and the distribution of Hamming distances between consecutive repeats P(ID(σ)), which is equivalent to maximize the likelihood of the data under the model. We fit these parameters using a gradient ascent algorithm: we start from an initial guess of the parameters, then generate sequences via Monte-Carlo simulations and update the parameters proportionally to the difference between the empirical and model generated observables , and P(ID(σ)) − P(ID(σ))model. We repeat the previous steps until the model reproduces the empirical observables defined above, with a target precision motivated according to the finite size of our original dataset, as in Ref. [17]. See Methods for more details. We tested the convergence of the model learning by synthetically generating datasets and relearning the model (see Methods).

By analogy with Boltmzan’s law, we call E(σ) a statistical energy, which is in general distinct from any physical energy. The particular form of the energy (2) resembles that of a disordered Potts model. This mathematical equivalence allows for the possibility to study effects that are characteristic of disordered systems, such as frustration or the existence of an energy landscape with multiple valleys, as we will discuss in the next sections.

Eq 2 is the most constrained form of the model, which we will denote by Efull(σ). One can explore the impact of each constraint on the energy landscape by removing them from the model. For instance, to study the role of inter-repeat sequence similarity due to a common evolutionary origin, one can fit the model without the constraint on repeat overlap ID, i.e. without the λID term in Eq 2. We call the corresponding energy function E2. One can further remove constraints on pairwise positions that are not part of the same repeat, making the two consecutive repeats statistically independent and imposing hi = hi+L (Eir), or only linked through phylogenic conservation through λID (Eir,λ). Finally one can remove all interaction constraints to make all positions independent of each other (E1), or even remove all constraints (Erand ≡ 0).

Statistical energy vs unfolding energy

The evolutionary information contained in multiple sequence alignments of protein families is summarized in our model by the energy function E(σ). Since this information is often much easier to access than structural or functional information, there is great interest in extracting functional or structural properties from multiple sequence alignments, provided that there exists a clear quantitative relationship between statistical energy and physical energy.

Such a relationship was determined experimentally for repeat proteins by using E(σ) to predict the effect of point mutations on the folding stability measured by the free energy difference between the folded and unfolded states, ΔG, called the unfolding energy [17, 20]. Synthetic sequences with low E(σ) have also been shown to reproduce the fold and function of natural sequences [29]. Here, extending an argument already developed in previous work [30–33], we show how this correspondance between statistical likelihood and folding stability arises in a simple model of evolution.

Evolutionary theory predicts that the prevalence of a particular genotype σ, i.e. the probability of finding it in a population, is related to its fitness F(σ). In the limit where mutations affecting the protein are rare compared to the time it takes for mutations to spread through the population, Kimura [34] showed that the probability of a mutation giving a fitness advantage (or disadvantage depending on the sign) ΔF over its ancestor will fix in the population with probability 2ΔF/(1 − e−2NΔF), where N is the effective population size. The dynamics of successful substitution satisfies detailed balance [35], with the steady state probability(3)

Again, one may recognize a formal analogy with Boltzmann’s distribution, where F plays the role of a negative energy, and N an inverse temperature. If we now assume that fitness is determined by the unfolding free energy ΔG, F(σ) = f(ΔG(σ)), then the distribution of genotypes we expect to observe in a population is(4)

Note that a similar relation should hold even if we relax the hypotheses of the evolutionary model. While in more general contexts (e.g. high mutation rate, recombination), the relation between ln P(σ) and F(σ) may not be linear, such nonlinearities could be subsumed into the function f.

Identifying terms in the two expressions (1) and (3), we obtain a relation between the statistical energy E, and the unfolding free energy ΔG:(5)

For instance, if we assume a linear relation between fitness and ΔG, f(ΔG) = A + BΔG, then we get a linear relationship between the statistical energy and ΔG, as was found empirically for repeat proteins [17].

Strikingly, the relationship f does not have to be linear or even smooth for this correspondance to work. Imagine a more stringent selection model, where f(ΔG) is a threshold function, f(ΔG) = 0 for ΔG > ΔGsel and −∞ otherwise (lethal). In that case the probability distribution is P(σ) = (1/Z)Θ(ΔG − ΔGsel), where Θ(x) is Heaviside’s function. Using a saddle-point approximation, one can show that in the thermodynamic limit (long proteins, or large L) the distribution concentrates at the border ΔGsel, and is equivalent to a “canonical” description [30, 31, 33]:(6)where the “temperature” Tsel is set to match the mean ΔG between the two descriptions:(7)

This correspondance is mathematically similar to the equivalence between the micro-canonical and canonical ensembles in statistical mechanics.

Statistical energy and unfolding free energy are linearly related by equating (Eq 1) and (Eq 6):(8)despite f being nonlinear. Eq 8 is in fact very general and should hold for any f in the thermodynamic limit in the vicinity of ⟨E⟩.

Equivalence between two definitions of entropies

There are several ways to define the diversity of a protein family. The most intuitive one, followed by [26], is to count the total number of amino acid sequences that have an unfolding free energy ΔGsel above a threshold ΔGsel [2]. This number naturally defines a Boltzmann entropy,(9)

Alternatively, starting from a statistical model P(σ), one can calculate its Shannon entropy, defined as(10)as was done in Ref. [27]. What is the relation between these two definitions?

By the same saddle-point approximation as in the previous section, the two are identical in the thermodyamic limit (large L), provided that the condition (Eq 7) is satisfied. We can thus reconcile the two definitions of the entropy in that limit.

To calculate the Boltzmann entropy (Eq 9), one needs to first evaluate the threshold Esel in terms of statistical energy. This threshold is given by Esel = E0 − ΔGsel/Tsel, where E0 and Tsel can be obtained directly by fitting (Eq 8) to single-mutant experiments. Esel can also be obtained as a discrimination threshold separating sequences that are known to fold properly versus sequences that do not [26]. In that case, assuming that the linear relationship (Eq 8) was evaluated empirically using single mutants, this relationship can be inverted to get ΔGsel in physical units.

Calculating the Shannon entropy Eq (10), on the other hand, does not require to define any threshold. However, the threshold in the equivalent Boltzmann entropy can be obtained using Eqs 7 and 8, i.e. Esel = ⟨E⟩, where the average is performed using the distribution defined in Eqs 1 and 2.

Entropy of repeat protein families

To compare how the different elements of the energy function affect diversity, we calculate the entropy of ensembles built of two consecutive repeats from a given protein family for the different kinds of models described earlier, from the least constrained to the most constrained: Erand, E1, Eir, Eir,λ, E2, Efull. In the case of models with interactions, calculating the entropy directly from the definition Eq (10) is impossible due to the large sums. A previous study of entropies of protein families used an approximate mean-field algorithm, called the Adaptive Cluster Expansion [27], for both parameter fitting and entropy estimation. Here we estimated the entropies using thermodynamic integration of Monte-Carlo simulations, as detailed in Methods. This method is expected to be asymptotically unbiased and accurate in the limit of large Monte-Carlo samples.

The resulting entropies and their differences are reported in Table 1 and Fig 2. All three considered families (ankyrins (ANK), leucine-rich repeats (LRR), and tetratricopeptides (TPR)) show a large reduction in entropy (∼40 − 50%) compared to random polypeptide string models of the same length 2L (of entropy Srand = 2L ln(21)). Interactions and phylogenic similarity between repeats generally have a noticeable effect on family diversity, although the magnitude of this effect depends on the family: (S1 − Sfull)/Sfull = 7% for ANK, versus, 13% for LRR, and 16% for TPR. Thus, although interactions are essential in correctly predicting the folding properties, they seem to only have a modest effect on constraining the space of accessible proteins compared to that of single amino-acid frequencies. However, when converted to numbers of sequences, this reduction is substantial, from to for ANK, from 1039 to 1034 for LRR, and from 7 ⋅ 1050 to 4 ⋅ 1042 for TPR.

[Figure omitted. See PDF.]

Table 1. Entropies (in bits, i.e. units of ln(2)) of sequences made of two consecutive repeats, for the three protein families shown in Fig 1.

Entropies are calculated for models of different complexity: model of random amino acids (Srand = 2L ln(21), divided by ln(2) when expressed in bits); independent-site model (S1), pairwise interaction model (S2); pairwise interaction model with constraints due to repeat similarity λID (Sfull); pairwise interaction model of two non-interacting repeats learned without (Sir) and with (Sir,λ) constraints on repeat similarity. Fig 2 shows graphically some of the information contained in this table.

https://doi.org/10.1371/journal.pcbi.1007282.t001

[Figure omitted. See PDF.]

Fig 2. Contributions of within-repeat interactions (S1 − Sir green), repeat-repeat interactions (Sir,λ − Sfull, purple), and phylogenic bias between consecutive repeats (Sir − Sir,λ, blue), to the entropy reduction from an independent-site model.

All three contributions are comparable, but with a larger effect of within-repeat interactions and phylogenic bias in TPR. The fourth bar (orange) quantifies the redundancy between two constraints with overlapping scopes: the constraint on consecutive-repeat similary, and the constraint on repeat-repeat correlations. This redundancy is naturally measured within information theory by the difference of impact (i.e. entropy reduction) of a constraint depending on whether or not the other constraint is already enforced.

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

By considering models with more and more constraints, and thus with lower and lower entropy, we can examine more finely the contribution of each type of correlation to the entropy reduction, going from E1 to Eir to Eir,λ to Efull. This division allows us to quantify the relative importance of phylogenic similarity between consecutive repeats (λID) relative to the impact of functional interactions (Jij), as well as the relative weights of repeat-repeat versus within-repeat interactions (Fig 2). We find that phylogenic similarity contributes substantially to the entropy reduction, as measured by Sir − Sir,λ = 4.5 bits for ANK, 4.3 bits for LRR, and 10.7 bits for TPR. The contribution of repeat-repeat interactions (Sir,λ − Sfull ∼ 5 bits for all three families) is comparable or of the same order of magnitude as that of within-repeat interactions (S1 − Sir = 4.3 bits for ANK, 6.9 bits for LRR, and 11.4 bits for TPR). This result emphasizes the importance of physical interactions between neighboring repeats in the whole protein.

On a technical note, we also find that pairwise interactions encode constraints that are largely redundant with the constraint of phylogenic similarity between consecutive repeats, as can be measured by the double difference Sir − Sir,λ − S2 + Sfull > 0 (Fig 2, orange bars). This redundancy comes from the fact that, in absence of an explicit constaint on P(ID) in E2, the interaction couplings Ji,i+L(σ, σ) between homologous positions in the two repeats is expected to favor pairs of identical residues to mimic the effect of λID. This redundancy motivates the need to correct for this phylogenic bias before estimating repeat-repeat interactions.

Comparing the three families, ANK has little phylogenic bias between consecutive repeats, and relatively weak interactions. By contrast, TPR has a strong phylogenic bias and strong within-repeat interactions.

Effect of interaction range

We wondered whether interactions constraining the space of accessible proteins had a characteristic lengthscale. To answer this question, for each protein family in Fig 1, we learn a sequence of models of the form Eq 2, in which Jij was allowed to be non-zero only within a certain interaction range d(i, j) ≤ W, where the distance d(i, j) between sites i and j can be defined in two different ways: either the linear distance |i − j| expressed in number of amino-acid sites, or the three-dimensional distance between the closest heavy atoms in the reference structure of the residues. Details about the learning procedure and error estimation are given in the Methods; see also S1 Fig for an alternative error estimate.

The entropy of all families decreases with interaction range W, both in linear and three-dimensional distance, as more constraints are added to reduce diversity (Fig 3 for ANK, and S2 Fig for LRR and TPR). The initial drop as a function of linear distance (Fig 3A) is explained by the many local interactions between nearby residues in the sequence. The entropy then plateaus until interactions between same-position residues in consecutive repeats are included in the W range, which leads to a sharp entropy drop at W = L. This suggests that long range interactions along the sequence generally do not constrain the protein ensemble diversity, except for interactions at exactly the scale of the repeat. This result suggests that the repeat structure is an important constraint limiting protein sequence exploration. These observations hold for all three repeat protein families. The importance of 3D structure in reducing the entropy can also be appreciated in the entropy decay as a function of physical distance (Fig 3B for ANK) where most of the entropy drop happens within the first 10 angstroms, indicating that above this characteristic distance interactions are not crucial in constraining the space of accessible sequences.

[Figure omitted. See PDF.]

Fig 3. Entropy reduction as a function of the range of interactions between residue sites.

A) Entropy of two consecutive ANK repeats, as a function of the maximum allowed interaction distance W along the linear sequence. The entropy of the model decreases as more interactions are added and they constrain the space of possible sequences. After a sharp initial decrease at short ranges, the entropy plateaus until interactions between complementary sites in neighbouring repeats lead to a secondary sharp decrease at W = L − 1 = 32 (dashed line), due to structural interactions between consecutive repeats. B) Entropy of two consecutive ANK repeats as a function of the maximum allowed three-dimensional interaction range. The entropy decreases rapidly until ∼10 Angstrom, after which decay becomes slower. In both panels entropies are averaged over 10 realizations of fitting the model; See Methods for details of the learning and entropy estimation procedure. Error bars are estimated from fitting errors between the data and the model; see Methods and S1 Fig for error bars calculated as standard deviations over 10 realizations of model fitting.

https://doi.org/10.1371/journal.pcbi.1007282.g003

Multi-basin structure of the energy landscape

The energy function of Eq (2) takes the same mathematical form as a disordered Potts model. These models, in particular in cases where σi can only take two values, have been extensively studied in the context of spin glasses [36]. In these systems, the interaction terms −Jij(σi, σj) imply contradictory energy requirements, meaning that not all of these terms can be minimized at the same time—a phenomenon called frustration. Because of frustration, natural dynamics aimed at minimizing the energy are expected to get stuck into local, non-global energy minima (Fig 4), significantly slowing down thermalization. This phenomenon is similar to what happens in structural glasses in physics, where the energy landscape is “rugged” with many local minima that hinder the dynamics. Incidentally, concepts from glasses and spin glasses have been very important for understanding protein folding dynamics [37].

[Figure omitted. See PDF.]

Fig 4. A rugged energy landscape is characterized by the presence of local minima, where proteins sequences can get stuck during the evolutionary process.

The set of sequences that evolve to a given local minimum defines the basin of attraction of that mimimum.

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

We asked whether the energy landscape of Eq (2) was rugged with multiple minima, and investigated its structure. To find local minima, we performed a local energy minimization of Efull (learned with all constraints including on P(ID), but taken with λID = 0 to focus on functional energy terms). By analogy with glasses, such a minimization is sometimes called a zero-temperature Monte-Carlo simulation or a “quench”. The minimization procedure was started from many initial conditions corresponding to naturally occuring sequences of consecutive repeat pairs. At each step of the minimization, a random beneficial (energy decreasing) single mutation is picked; double mutations are allowed if they correspond to twice the same single mutation on each of the two repeats. Minimization stops when there are no more beneficial mutations. This stopping condition defines a local energy minimum, for which any mutation increases the energy. The set of sequences which, when chosen as initial conditions, lead to a given local minimum defines the basin of attraction of that energy mimimum (Fig 4). The size of a basin corresponds to the number of natural proteins belonging to that basin.

Performing this procedure on natural sequences of consecutive repeat pairs from all three families yielded a large number of local minima (Fig 5). To control for the phylogenetic bias that links natural sequences, we repeated this analysis on sequences synthetically generated from the model (Efull), and obtained very similar results (see S6 Fig for ANK). When ranked from largest to smallest, the distribution of basin sizes follows a power law (Fig 5A for ANK and S3A and S4A Figs for LRR and TPR). The energy of the minimum of each basin generally increases with the rank, meaning that largest basins are also often the lowest. Despite this multiplicity of local minima, the Monte-Carlo dynamics that we used in previous sections for learning the model parameters and for estimating the entropy did not get stuck in these minima, suggesting relatively low energy barriers between them.

[Figure omitted. See PDF.]

Fig 5. Interactions within and between repeats sculpt a rugged energy landscape with many local minima.

Local minima were obtained by performing a zero-temperature Monte-Carlo simulation with the energy function in Eq (2), starting from initial conditions corresponding to naturally occurring sequences of pairs of consecutive ANK repeats. A, bottom) Rank-frequency plot of basin sizes, where basins are defined by the set of sequences falling into a particular minimum. A, top) energy of local minima vs the size-rank of their basin, showing that larger basins often also have lowest energy. Gray line indicates the energy of the consensus sequence, for comparison. B) Pairwise distance between the minima with the largest basins (comprising 90% of natural sequences), organised by hierarchical clustering. The panel right above the matrix shows the size of the basins relative to the minima corresponding to the entries of the distance matrix. A clear block structure emerges, separating different groups of basins with distinct sequences. C-D) Same as A) and B) but for single repeats. Since single repeats are shorter than pairs (length L instead of 2L), they have fewer local energy minima, yet still show a rich multi-basin structure. Equivalent analyses for LRR and TPR are shown in S3 and S4 Figs.

https://doi.org/10.1371/journal.pcbi.1007282.g005

The partition of sequences into basins allows for the definition of a new kind of entropy Sconf = −∑b P(b) ln P(b) called configurational entropy, based on the distribution of basin sizes, P(b) = ∑σ∈b P(σ), where σ ∈ b means that energy minimization starting with sequence σ leads to basin b. This configurational entropy measures the effective diversity of basins, and is thus much lower than the sequence entropy Sfull, while the difference Sfull − Sconf measures the average diversity of sequences within each basin. We find Sconf = 5.1 bits for ANK, 6.0 bits for LRR, and 10.4 bits for TPR. As each basin corresponds to a distinct sub-family within each family [32], this entropy quantifies the effective number of these subgroups.

While basins are very numerous, they are also not independent of each other. An analysis of pairwise distances (measured as the Hamming distance between the local minima) between the largest basins reveals that they can be organised into clusters (panels B of Fig 5, S3 and S4 Figs), suggesting a hierarchical structure of basins, as is common in spin glasses [36].

The impact of repeat-repeat interactions on the multi-basin structure can be assessed by repeating the analysis on the model of non-interacting repeats, Eir. In that model the two repeats are independent, so it suffices to study local energy minima of single repeats—local minima of pairs of repeats follow simply from the combinatorial pairing of local minima in each repeat. The analyses of basin size distributions, energy minima, and pairwise distances in single repeats are shown in panels C and D of Fig 5, S3, and S4 Figs. We still find a substantial number of unrelated energy minima, suggesting again several distinct subfamilies even at the single-repeat level. For comparison, the configurational entropy of pairs of independent repeats is 6.9 bits for ANK, 6.7 for LRR, and 7.6 for TPR. While for ANK and LRR repeat-repeat interactions decrease the configurational entropy, as they do for the conventional entropy, they in fact increase entropy for TPR, making the energy landscape even more frustrated and rugged.

Note that the independent sites model E1 defines a convex energy landscape with a single local minimum—the consensus sequence—as all constraints hi can be optimized independently. To address how the interactions contribute in shaping the sequence space, going from a convex to a rugged landscape, we repeated the analysis with a limited linear interaction range W of 3 and 10 (models of Fig 3A). We find that the more interactions we add, the more local minima we find (S5A and S5B Fig for ANK with W = 3, and C and D for W = 10). The minima cluster into clearer sub-blocks structure as the interaction range is increased, consistent with the entropy reduction observed in Fig 3A.

In summary, the analysis of the energy landscape reveals a rich structure, with many local minima ranging many different scales, and with a hierarchical structure between them.

Distance between repeat families

Lastly, we compared the statistical energy landscapes of different repeat families. Specifically, we calculated the Kullback-Leibler divergence between the probability distributions P(σ) (given by Eqs 1 and 2) of two different families, after aligning them together in a single multiple sequence alignment (see Methods).

We find essentially no similarity between ANK and TPR, despite them having similar lengths: DKL(ANK||TPR) = 227.6 bits, and DKL(TPR||ANK) = 214.1 bits. These values are larger than the Kullback-Leibler divergence between the full models for these families and a random polypetide, DKL(ANK||rand) = 122.8 bits, and DKL(TPR||rand) = 157.6 bits. LRR is not comparable to ANK or TPR as it is much shorter, and a common alignment is impractical. These large divergences between families of repeat proteins show that different families impose quantifiably different constraints, which have forced them to diverge into different troughs of non-overlapping energy landscapes. This lack of overlap makes it impossible to find intermediates between the two families that could evolve into proteins belonging to both families.

Discussion

Our analysis of repeat protein families shows that the constraints between amino acids in the sequences allows for an estimation of the size of the accessible sequence space. The obtained numbers (ranging from 141 bits to 167 bits, corresponding to 1036 to 1050 sequences) are of course huge compared to the number of sequences in our initial samples (∼20, 500 for ANK, ∼18, 800 for LRR, and ∼10, 000 for TPR), but comparable to the total number of proteins having been explored over the whole span of evolution, estimated to be 1043 in Ref. [1].

In particular, we have quantified the reduction of the accessible sequence space with respect to random polypeptides. While most of this reduction is attributable to conservation of residues at each site, interactions between amino acids, both within and between consecutive repeats, significantly constrain the diversity of all repeat families. The break-up of entropy reduction between the three different sources of constraints—within-repeat interactions, between-repeat interactions, and evolutionary conservation between consecutive repeats—is fairly balanced, although TPR stands out as having more within-repeat interactions and more conservation between neighbours, suggesting that it may have had less time to equilibrate.

All studied repeat families have rugged energy landscapes with multiple local energy minima. Note that the emergence of this multi-valley landscape is a consequence of the interactions between amino acids: models of independent positions (E1) only admit a single energy minimum corresponding to the consensus sequence. This multiplicity of minima allow us to collapse multiple sequences to a small number of coarse-grained attractor basins. These basins suggest that mutations between sequences within one coarse-grained basin are much more likely than mutating into sequences in other basins. In general, our results paint a picture of further subdivisions within a family, and define sub-families due to the fine grained interaction structure. Going beyond single families, this analysis suggest a view in which natural proteins all live in a global evolutionary landscape, of which families would be basins, or clusters of basins, with a hierarchical structure [32].

This overall picture of the sequence energy landscape is reminiscent of the hierarchical picture of the structural energy landscape of globular proteins, an overall funneled shape with tiers within tiers [38]. The form of the energy landscape forcibly shapes the accessible evolutionary paths between sequences. The rugged and further subdivided structure shows that the uncovered constraints are global, and not just pairwise between specific residues. Therefore even changing two residues together, as is often done in laboratory experiments, is not enough to recover the evolutionary trajectories. While other approaches have explored local accessible directions of evolution [39], our results suggest more global, non local modes of evolution between clusters.

Interestingly, the sequences that correspond to the energy minima of the landscapes are not found in the natural dataset. This observation can be either due to sampling bias (we have not yet observed the sequence with the minimal energy, although it exists), or this sequence may not have been sampled by nature. Alternatively, there may be additional functional constraint that are not included in our model to avoid these low energy sequences (e.g. a too stable protein may be difficult to degrade).

Even more intriguingly, sequences with minimal energy do not correspond to the consensus sequence of the alignment (whose energy is marked by a gray line in panel A of Fig 5, S3, and S4 Figs), suggesting that the consensus sequence can be improved upon. All three repeat protein families studied here have been shown to be amenable to simple consensus-guided design of synthetic proteins. Synthetic proteins based on the consensus sequences of multiple alignments [40] were found to be foldable and very stable against chemical and thermal denaturation. Mutations towards consensus amino acids in the ANK family members have been experimentally shown to both stabilize the whole repeat-array and they may tune the folding paths towards nucleating folding in the consensus sites [41, 42]. Our results suggest that interactions may play an additional role in stabilizing the sequences, and propose alternative solutions to the consensus sequences in the design of synthetic proteins.

Methods

Data curation

We use a previously curated alignment of pairs of repeats for each family [17]: ANK (PFAM id PF00023 with a final alignment of 20513 sequences of L = 66 residues each), LRR (PFAM id PF13516 with a final alignment of 18839 sequences of L = 48 residues each) and TPR (PFAM id PF00515 with a final alignment of 10020 sequences of L = 68 residues each). Those multiple sequence alignments of repeats were obtained from PFAM 27.0 [6, 7]. In order to improve the data obtained from the PFAM database, we used original full protein sequences available in UniProt database [43] to add available information using the headers of the original alignement. Firstly, to decrease the number of gaps positions, misdetected initial and final amino acids in repeats were completed with residues from full sequences. Secondly, individual repeats which appeared consecutively in natural proteins were joined into pairs. Finally, positions with more than 80% of gaps along the alignment were removed, eliminating in this way insertions.

From the multiple sequence alignement of each family, they were calculated the observables that we use to constrain our statistical model. Particularly, we calculated the marginal frequency fi(σi) of an amino acid σi at position i and the joint frequency fij(σi, σj) of two amino acids σi and σj at two different positions i and j. These quantities were calculated using only sequences selected by clustering at 90% of identity computed with CD-HIT [44] and then normalizing by the amount of sequences. In this way, the occurrences of residues in every position are not biased by overrepresentation of proteins in the database. Furthermore, to take into account the repeated nature of the protein families that we are considering, an additional observable was calculated, the distribution of sequence overlap between two consecutive repeats, P(ID(σ)), with .

Model fitting

In order to obtain a model that reproduces the experimentally observed site-dependent amino-acid frequencies, fi(σi), correlations between two positions, fij(σi, σj), and the distribution of Hamming distances between consecutive repeats, P(ID(σ)), we apply a likelihood gradient ascent procedure, starting from an initial guess of the hi(σi), Jij(σi, σj) and λID(σ) parameters.

At each step, we generate 80000 sequences of length 2L through a Metropolis-Hastings Monte-Carlo sampling procedure. We start from a random amino-acid sequence and we produce many point mutations in any position, one at a time. If a mutation decreases the energy (2) we accept it. If not, we accept the mutation with probability e−ΔE, where ΔE is the difference of energy between the original and the mutated sequence. We add one sequence to our final ensemble every 1000 steps. Once we generated the sequence ensemble, we measure its marginals and , as well as Pmodel(ID(σ)), and update the parameters of Eq 2 following the gradient of the likelihood. The local field and λID(σ) are updated along the gradient of the per-sequence log-likelihood, equal to the difference between model and data averages:(11) (12)

As the number of parameters for the interaction terms Jij is large (= 212 L2), we force to 0 those that are not contributing significantly to the model frequencies through a L1 regularisation γ∑ij, σ, τ|Jij(σ, τ)| added to the likelihood. This leads to the following rules of maximization: If Jij(σi, σj)t = 0 and (13)

If Jij(σi, σj)t = 0 and (14)

If (15)

If (16)

The optimization parameters were set to: ϵm = 0.1, ϵj = 0.05, ϵID = 10, and γ = 0.001.

To estimate the model error, we compute and . We also calculate the difference of generated and natural repeat similarity distribution for all the possible repeats Hamming distances, penalized by a factor 5 to better learn the parameter λID: 5(P(ID(σ)) − P(ID(σ))model). We repeat the procedure above until the maximum of all errors, , and 5|P(ID(σ)) − P(ID(σ))model|, goes below 0.02, as in Ref. [17].

Models with different sets of constraints

Using this procedure we can calculate the model defined in Eq 2 with different interaction ranges used in the entropy estimation in Fig 3A. We start from the independent model hi(σi) = log fi(σi). We first learn the model in Eq 2 with J = 0. We then re-learn models with interactions between sites i, j along the linear sequence such that |i − j| ≤ W, in a seeded way starting from the previous model. The first and last point of Fig 3 correspond to the independent site model with λID and the full model in Eq 2

The entropy in Fig 3B is calculated in the same way as in Fig 3, but now interactions are turned on progressively according to physical distance in the 3D structure rather than the linear sequence distance. In order to obtain the physical distance between residues we use as a reference structure the first two repeats of a consensus designed ankyrin protein 1n0r [45, 46], which have exactly 66 amino-acids. We define the 3D separation between two residues as the minimum distance between their heavy atoms in the reference structure.

To learn the Potts model without λID (E2) we remove λID from Eq 2 and re-learn the Potts field using the full model parameters as initial contition.

To learn the single repeat models with and without λ (Eir and Eir,λ, we take as initial condition the model with interactions below the length of a repeat (W = L − 1, dashed vertical line in Fig 3), and then learn a model removing all the Jij terms between different repeats. We also impose that the hi fields and intra-repeats Jij terms are the same in each repeat, and the experimental amino-acid frequencies to be reproduced by the model are the average over the two repeats of the 1- and 2-points intra-repeats frequencies fi(σi) and fij(σi, σj), such that(17)and(18)if i and j represent sites within the same repeat. In this way we obtain a model for a single repeat that can be extended to both the repeats in the original set of sequences of our dataset.

Entropy estimation

In practice to calculate the entropy S of the protein families we relate it to the internal energy E = −log p(σ) and the free energy F = −log Z:(19)

We generate sequences according to the energy function in Eq 2 and use them to numerically compute ⟨E⟩. To calculate the free energy we use the auxilliary energy function:(20)where the interaction strength across different sites can be tuned through a parameter α that is changed from 0 to 1. We generate protein sequence ensembles with different values of α and use them to calculate F as a function of α, :(21)where the average over α is taken over the sequences generated with a certain value of α, characterized by the ensemble with probability . F(0) is the free energy for an independent sites model:(22)where the first sum is taken over protein sites and the second over all possible amino-acids at a given site. Eqs 22 and 19 result in the thermodynamic sampling approximation for calculating the entropy [47]:(23)

We generate 80000 sequences using Monte Carlo sampling for the energy in Eq 20 with 50 different α values, equally spaced between 0 and 1 at a distance of 0.02, and then numerically compute the integral in Eq 23 using the Simpson rule.

Entropy error

The entropy estimate is subject to three sources of uncertainty: the finite-size of the dataset, convergence of parameter learning, and the noise in the thermodynamic integration. We estimate the contribution of each of these errors using the independent sites model. In the independent sites model each site i is simply described by a multinomial distribution with weights given by the observed amino-acid frequencies in the datasets. The variance in the estimation of the frequencies from a finite size sample is Var(fi(σi)) = (pi(σi)(1 − pi(σi)))/Ns and the covariance between the frequencies of different amino-acids σ and σ′ at the same site i is where Ns is the sample size and pi(σi) are the weights of the true multinomial distribution sampled. Through error propagation from these quantities we calculate the variance in the entropy of the independent sites model, to first order in 1/Ns:(24)

Evaluating this equation using the empirical frequencies p = f assuming they are sampled from an underlying multinomial distribution, gives an estimate of the standard deviation of 0.05. We assume that the interaction terms do not change the order of magnitude of this estimation. Also the standard deviation in the averages in Eq (23) scales as with Ns = 80000.

The parameter inference is affected not only by noise, but also by a systematic bias depending on the parameters of the gradient ascent described above and the initial condition that we chose to start learning from. S1 Fig shows the average entropy of 10 realizations of the learning and thermodynamic integration procedure for the ANK family and its standard deviation as error bars. If we learn the models with an increasing W window progressively we get a different profile than learning each point starting from the independent model, and above L these two profiles are more distant than the magnitude of the standard deviation, signalling a systematic bias. S1 Fig also shows that progressively learning the model results in a better parameters convergence to values that give lower entropy values.

In order to estimate how this bias is reflected in the entropy estimation we take the single-site amino-acid frequencies produced by the inferred energy function in the last Monte-Carlo phase of the learning procedure and calculate the corresponding entropy for this independent-sites model. We compute the absolute value of the difference between this estimate of the entropy and the independent-sites entropy calculated from the dataset. Again in doing this we assume that neglecting the interaction terms does not change the order of magnitude of this error. These procedure results in the errorbars shown in Figs 3 and 2, Table 1, S2 Fig.

We repeat 10 realizations of both the parameter inference procedure and the entropy estimation, and in Fig 3 we show the average entropy of these 10 numerical experiments for the ANK family where error bars are estimated as explained above to sketch the order of magnitude of the error coming from systematic bias in the parameters learning. S1 Fig shows the mean entropy of ANK as in Fig 3A with the standard deviations of the realizations entropy as error bars, to give an idea of the combined noise in the thermodynamic integration and in the gradient descent, starting from the same initial conditions and with the same update parameters (see Section). The combined noise is smaller than the entropy decrease at 33 residues, showing the decrease is real.

To further check the robustness of the entropy estimation procedure, we generate two synthetic ANK datasets, one with an independent sites model, the other with a model of two non-interacting repeats obtained as explained above, and relearn the model from the synthetic datasets. Repeating the learning and entropy estimation procedure on each on the synthetic protein families gives results that are consistent with the model used for the dataset generation. The entropy of the model learned taking an independent sites dataset does not decrease with the interaction range W and the entropy of the model learned taking a non-interacting repeats dataset does not show any drop around the repeat length.

We repeat the procedure described for the LRR and TPR repeat-proteins families reaching similar conclusions (S2 Fig).

Calculating the basins of attraction of the energy landscape

In order to characterize the ruggedness of the inferred energy landscapes and the sequence identity of the local minima, we start from all the sequences in the natural dataset as initial conditions and for each of them we perform a T = 0 quenched Monte-Carlo procedure. Repeating this analysis on sequences synthetically generated from Efull yields very similar results (see S6 Fig for ANK)

We perform this energy landscape exploration learning the parameters of the Hamiltonian in Eq 2 (refer to Section for the learning procedure), and then set λID = 0 in the energy function because we want to investigate the shape of the energy landscape due to selection rather than the phylogenic dependence.

We scan all the possible mutations that decrease the sequence energy and then draw one of them from a uniform random distribution. The possible mutations are all single point mutations. If the same amino-acid is present in the same relative position in the two repeats we allow for double mutations that mutate those two positions to a new amino-acid, that is identical in both repeats, at the same time. We do this so that the phylogenetic biases that are still partially present in the parameters of the model do not result in spurious local minima biasing the quenching results. The Monte-Carlo procedure ends when every proposed move results in a sequence with an increased energy, and the identified sequence is a local minimum of the energy landscape.

To explore how turning on interactions makes the energy landscape more rugged, we perform the same procedure with the Hamiltonian corresponding to two intermediate interaction ranges in Fig 3A. That is Eq 2, in which Jij was allowed to be non-zero only within a certain interaction range W. We picked W = 3 and W = 10.

In order to assess what is the role of the inter-repeat interactions we repeat this T = 0 quenched Monte-Carlo procedure on single repeats, with all the unique repeats in the natural dataset as initial condition. The learning procedure of the Hamiltonian for a single repeat is explained in Section. In this single repeat case the possible mutations are just the single point mutations.

Once we have the local minima of the energy landscape, we obtain the coarse-grained minima using the Python Scipy hierarchical clustering algorithm. In this hierarchical clustering the distance between two clusters is calculated as the average Hamming distance between all the possible pairs of sequences belonging each to one cluster. As a result we plot the clustered distance matrix, the clustering dendogram and the basin size corresponding to the distance matrix entries.

In the end we can repeat the quenching procedure described above for LRR and TPR families. The result are sketched in S3 and S4 Figs and lead to similar conclusions as for the ANK family.

Kullback-Leibler divergence

The Kullback-Leibler divergence between two families A and B is defined as DKL(A||B) = ∑σ pA(σ) log2 pB(σ)/pA(σ). We can substitute the sequence ensembles for ANK and TPR in the definition of the probabilities obtaining:(25) (26)where the notation ⟨⟩ANK means that the average is calculated over sequences drawn from the ANK ensemble: . Therefore ⟨ETPR⟩ANK is the average TPR energy function evaluated, via the structural alignment between the two families, on 80000 sequences generated through a Monte Carlo sampling of the ANK model 2 (and analogously for ⟨EANK⟩TPR). The terms FANK and FTPR are calculated in the same way as when estimating the entropy through Eqs 21 and 22, as explained in Section.

For the control against a random polypeptide of length L we use DKL(FAM||rand) = log Λ − S(FAM), where Λ = 21L is the total number of possible sequences of length L.

Supporting information

S1 Fig. Reproducibility of entropy estimation.

Entropy as a function of the maximum linear interaction range W along the sequence. Green curve: entropy of the ANK family with error bars calculated as standard deviations over 10 model learning realizations, where models are learned by incrementally adding more interaction terms as W is increased, taking the model learned at W − 1 as initial condition. This plot is the same as in Fig 3A but with the different error bar estimates, showing that our results are robust to the details of error estimation. Red curve: entropy obtained after de novo learning for each W, starting from a non-interacting model as initial condition. With those initial conditions the learning gets stuck, leading to systematically overestimating the entropy and missing the second entropy drop at W = L − 1. See Section for details of the learning and entropy estimation procedure.

https://doi.org/10.1371/journal.pcbi.1007282.s001

(TIFF)

S2 Fig. Range dependence of entropy in LRR and TPR families.

Entropy of the LRR (A) and TPR (B) family as a function of the maximum interaction distance W along the sequence. The entropy of the model decreases as a more interactions are added and they constrain the space of possible sequences. As with ANK, the entropy first drops, plateaus, then drops again at the distance corresponding to homologous positions along the two repeats (W = L − 1 = 23 for LRR, and 33 for TPR, dashed line). This second drop indicates that there is a typical distance along the sequence, corresponding to the repeat length, where interactions due to structural properties constrain the sequence ensemble. The error bars are estimated approximately from errors in learning (see Section). Entropies are averaged over 5 realizations of the learning and entropy estimation procedure.

https://doi.org/10.1371/journal.pcbi.1007282.s002

(TIFF)

S3 Fig. Analysis of local energy minima for pairs of consecutive repeats of LRR.

Energy minima were obtained by zero-temperature dynamics. Sequences falling into a given minimum with these dynamics define its basin of attraction. A, bottom) rank-frequency plot of the sizes of the basins of attraction. A, top) energy minimum of each basin. Gray line shows the energy of the consensus sequence B) Pairwise Hamming distances between energy minima, organised by hierarchical clustering. The panel right above the matrix shows the the size of the basins relative to the minima corresponding to the entries of the distance matrix. C and D) Same analysis as A) and B), but for single LRR repeats.

https://doi.org/10.1371/journal.pcbi.1007282.s003

(TIFF)

S4 Fig. Analysis of local energy minima for pairs of consecutive repeats of TPR.

Energy minima were obtained by zero-temperature dynamics. Sequences falling into a given minimum with these dynamics define its basin of attraction. A, bottom) rank-frequency plot of the sizes of the basins of attraction. A, top) energy minimum of each basin. Gray line shows the energy of the consensus sequence B) Pairwise Hamming distances between energy minima, organised by hierarchical clustering. The panel right above the matrix shows the the size of the basins relative to the minima corresponding to the entries of the distance matrix. C and D) Same analysis as A) and B), but for single TPR repeats.

https://doi.org/10.1371/journal.pcbi.1007282.s004

(TIFF)

S5 Fig. Interactions within repeats increase the ruggedness of the energy landscape.

Local minima were obtained by performing a zero-temperature Monte-Carlo simulation with the energy function in Eq (2) with non-zero Jij within linear interaction range W, starting from initial conditions corresponding to naturally occurring sequences of pairs of consecutive ANK repeats, for W = 3 (A and B) and W = 10 (C and D). See Fig 5 for the full model (W = 2L). A and C, bottom: Rank-frequency plot of basin sizes, where basins are defined by the set of sequences falling into a particular minimum. A and C, top: energy of local minima vs the size-rank of their basin. Gray line indicates the energy of the consensus sequence, for comparison. B and D: Pairwise distance between the minima with the largest basins (comprising 90% of natural sequences), organised by hierarchical clustering. The panel right above the matrix shows the size of the basins relative to the minima corresponding to the entries of the distance matrix. The block structure starts emerging as interactions are turned on (D versus B).

https://doi.org/10.1371/journal.pcbi.1007282.s005

(TIFF)

S6 Fig. Analysis of local energy minima from generated pairs of consecutive repeats of ANK.

Energy minima were obtained by zero-temperature dynamics starting from sequences generated in silico from Efull. Sequences falling into a given minimum with these dynamics define its basin of attraction. A, bottom) rank-frequency plot of the sizes of the basins of attraction. A, top) energy minimum of each basin. Gray line shows the energy of the consensus sequence B) Pairwise Hamming distances between energy minima, organised by hierarchical clustering. The panel right above the matrix shows the size of the basins relative to the minima corresponding to the entries of the distance matrix. C and D) Same analysis as A) and B), but for single ANK repeats.

https://doi.org/10.1371/journal.pcbi.1007282.s006

(TIFF)

Citation: Marchi J, Galpern EA, Espada R, Ferreiro DU, Walczak AM, Mora T (2019) Size and structure of the sequence space of repeat proteins. PLoS Comput Biol 15(8): e1007282. https://doi.org/10.1371/journal.pcbi.1007282

References

1. Dryden DTF, Thomson AR, White JH. How much of protein sequence space has been explored by life on Earth? Journal of the Royal Society InterfaceRoyal Society Interface. 2008;5(April):953–956.

2. Shakhnovich EI. Protein design: a perspective from simple tractable models. Current Biology. 1998;3:45–58.

3. Salisbury FB. Natural Selection and the Complexity of the Gene. Nature. 1969;244:342–343.

4. Mandecki W. The game of chess and searches in protein sequence space. Biotopic. 1998;16(May):200–202.

5. Dokholyan NV, Shakhnovich B, Shakhnovich EI. Expanding protein universe and its origin from the biological Big Bang. PNAS. 2002;99:14132–14136. pmid:12384571

6. Bateman A, Coin L, Durbin R, Finn RD, Hollich V, Griffiths-Jones S, et al. The Pfam protein families database. Nucleic acids research. 2004;32(suppl 1):D138–D141. pmid:14681378

7. Finn RD, Bateman A, Clements J, Coggill P, Eberhardt RY, Eddy SR, et al. Pfam: the protein families database. Nucleic acids research. 2013; p. gkt1223.

8. Neher E. How frequent are correlated changes in families of protein sequences? Proceedings of the National Academy of Sciences. 1994;91(1):98–102.

9. Morcos F, Hwa T, Onuchic JN, Weigt M. Direct coupling analysis for protein contact prediction. Methods Mol Biol. 2014;1137:55–70. pmid:24573474

10. Szurmant H, Weigt M. Inter-residue, inter-protein and inter-family coevolution: bridging bridging the scales. Current Opinion in Structural Biology. 2018;50:26–32. pmid:29101847

11. Tubiana J, Cocco S, Monasson R. Learning protein constitutive motifs from sequence data. eLife. 2019;8:e39397:1–61. pmid:30857591

12. Weigt M, White RA, Szurmant H, Hoch JA, Hwa T. Identification of direct residue contacts in protein–protein interaction by message passing. Proceedings of the National Academy of Sciences. 2009;106(1):67–72.

13. Morcos F, Pagnani A, Lunt B, Bertolino A, Marks DS, Sander C, et al. Direct-coupling analysis of residue coevolution captures native contacts across many protein families. Proceedings of the National Academy of Sciences. 2011;108(49):E1293–E1301.

14. Ekeberg M, Lövkvist C, Lan Y, Weigt M, Aurell E. Improved contact prediction in proteins: using pseudolikelihoods to infer Potts models. Physical Review E. 2013;87(1):012707.

15. Hopf TA, Colwell LJ, Sheridan R, Rost B, Sander C, Marks DS. Three-dimensional structures of membrane proteins from genomic sequencing. Cell. 2012;149(7):1607–1621. pmid:22579045

16. Espada R, Parra RG, Mora T, Walczak AM, Ferreiro DU. Capturing coevolutionary signals inrepeat proteins. BMC bioinformatics. 2015;16(1):207. pmid:26134293

17. Espada R, Parra RG, Mora T, Walczak AM, Ferreiro DU. Inferring repeat-protein energetics from evolutionary information. PLoS computational biology. 2017; p. 1–16.

18. Schug A, Weigt M, Onuchic JN, Hwa T, Szurmant H. High-resolution protein complexes from integrating genomic information with molecular simulation. Proceedings of the National Academy of Sciences of the United States of America. 2009;106(52):22124–9. pmid:20018738

19. Marks DS, Colwell LJ, Sheridan R, Hopf TA, Pagnani A, Zecchina R, et al. Protein 3D structure computed from evolutionary sequence variation. PloS one. 2011;6(12):e28766. pmid:22163331

20. Contini A, Tiana G. A many-body term improves the accuracy of effective potentials based on protein coevolutionary data. The Journal of Chemical Physics. 2015;143(2):025103. pmid:26178131

21. Haldane A, Flynn WF, He P, Vijayan RSK, Levy RM. Structural propensities of kinase family proteins from a Potts model of residue co-variation. Protein Science;25(8):1378–1384. pmid:27241634

22. Figliuzzi M, Jacquier H, Schug A, Tenaillon O, Weigt M. Coevolutionary landscape inference and the context-dependence of mutations in beta-lactamase TEM-1. Molecular biology and evolution. 2015;. pmid:26446903

23. Maynard Smith J. Natural Selection and the Concept of a Protein Space. Nature. 1970;225:563–564.

24. Weinreich DM, Delaney NF, Depristo MA, Hartl DL. Darwinian evolution can follow only very few mutational paths to fitter proteins. Science (New York, NY). 2006;312(5770):111–114.

25. Li J, Mahajan A, Tsai Md, Li J, Mahajan A, Tsai Md. Current Topics / Perspectives Ankyrin Repeat: A Unique Motif Mediating Protein − Protein Interactions Ankyrin Repeat: A Unique Motif Mediating Protein—Protein Interactions †. Biochemistry. 2006;45(November):15168–15178.

26. Tian P, Best RB. How Many Protein Sequences Fold to a Given Structure? A Coevolutionary Analysis. Biophysj. 2017;113(8):1719–1730.

27. Barton JP, Chakraborty AK, Cocco S, Jacquin H, Monasson R. On the Entropy of Protein Families. Journal of Statistical Physics. 2016;162(5):1267–1293.

28. Jaynes ET. Information Theory and Statistical Mechanics. Phys Rev. 1957;106:620–630.

29. Tian P, Louis JM, Baber JL, Aniana A, Best RB. Co-Evolutionary Fitness Landscapes for Sequence Design. Angew Chemie—Int Ed. 2018;57(20):5674–5678.

30. Shakhnovich EI, Gutin AM. A new approach to the design of stable proteins. Protein Eng. 1993;6(8):793–800. pmid:8309926

31. Shakhnovich EI, Gutin AM. Engineering of stable and fast-folding sequences of model proteins. Proc Natl Acad Sci. 1993;90(15):7195–7199. pmid:8346235

32. Dokholyan NV, Shakhnovich EI. Understanding Hierarchical Protein Evolution from First Principles. Journal of Molecular Biology. 2001;312:289–307. pmid:11545603

33. Morcos F, Schafer NP, Cheng RR, Onuchic JN, Wolynes PG. Coevolutionary information, protein folding landscapes, and the thermodynamics of natural selection. 2014;111(34):12408–12413.

34. Kimura M. On the Probability of Fixation of Mutant Genes in a Population. Genetics. 1962;47:713–719. pmid:14456043

35. Berg J, Willmann S, Lässig M. Adaptive evolution of transcription factor binding sites. BMC evolutionary biology. 2004;4:42. pmid:15511291

36. Mezard M, Parisi G, Virasoro M. Spin Glass Theory and Beyond. WORLD SCIENTIFIC; 1986. Available from: https://www.worldscientific.com/doi/abs/10.1142/0271.

37. Bryngelson JD, Wolynes PG. Spin glasses and the statistical mechanics of protein folding. Proceedings of the National Academy of Sciences. 1987;84(21):7524–7528.

38. Frauenfelder H, Sligar SG, Wolynes PG. Proteins. Science. 1991;254:1598–1603.

39. Facco E, Pagnani A, Russo ET, Laio A. The intrinsic dimension of protein sequence evolution. PLoS Comput Biol. 2019;15(4):e1006767. pmid:30958823

40. Boersma YL, Plu A. DARPins and other repeat protein scaffolds: advances in engineering and applications. Current opinion in biotechnology. 2011;22:849–857. pmid:21715155

41. Tripp KW, Barrick D. Rerouting the Folding Pathway of the Notch Ankyrin Domain by Reshaping the Energy Landscape. Journal of the American Chemical Society. 2008;(19):5681–5688. pmid:18396879

42. Barrick D, Ferreiro DU, Komives EA. Folding landscapes of ankyrin repeat proteins: experiments meet theory. Current Opinion in structural biology. 2008;18:27–43. pmid:18243686

43. Consortium U, et al. UniProt: the universal protein knowledgebase. Nucleic acids research. 2017;45(D1):D158–D169.

44. Li W, Jaroszewski L, Godzik A. Tolerating some redundancy significantly speeds up clustering of large protein databases. Bioinformatics. 2002;18(1):77–82. pmid:11836214

45. Mosavi LK, Minor DL, Peng Zy. Consensus-derived structural determinants of the ankyrin repeat motif. Proceedings of the National Academy of Sciences. 2002;99(25):16029–16034.

46. Binz HK, Stumpp MT, Forrer P, Amstutz P, Plückthun A. Designing repeat proteins: well-expressed, soluble and stable proteins from combinatorial libraries of consensus ankyrin repeat proteins. Journal of molecular biology. 2003;332(2):489–503. pmid:12948497

47. Frankel D, Smit B. Understanding Molecular Simulation: From Algorithms to Applications. Academic Press; 2007.