Size and structure of the sequence space of repeat proteins

TL;DR

This paper estimates the size and explores the complex, rugged structure of the sequence space of repeat proteins, revealing hierarchical subtypes and correlations that influence diversity and have implications for protein design.

Contribution

It introduces a maximum entropy modeling approach to quantify the sequence space of repeat proteins and uncovers the hierarchical, rugged landscape structure of their sequence diversity.

Findings

A significant impact of amino acid correlations on sequence diversity

Identification of a rugged, hierarchical landscape with multiple local minima

Presence of subtypes within protein families based on sequence clustering

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.

Tables

Table 1. Table 1: Entropies (in bits, i.e. units of ln ⁡ ( 2 ) 2 \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 ( S rand = 2 ​ L ​ ln ⁡ ( 21 ) subscript 𝑆 rand 2 𝐿 21 S_{\rm rand}=2L\ln(21) , divided by ln ⁡ ( 2 ) 2 \ln(2) when expressed in bits); independent-site model ( S 1 subscript 𝑆 1 S_{1} ), pairwise interaction model ( S 2 subscript 𝑆 2 S_{2} ); pairwise interaction model with constraints due to repeat similarity λ ID subscript 𝜆 ID \lambda_{\rm ID} ( S full subscript 𝑆 full S_{\rm full} ); pairwise interaction model of two non-interacting repeats learned without ( S ir subscript 𝑆 ir S_{\rm ir} ) and with ( S ir , λ subscript 𝑆 ir 𝜆 S_{\rm ir,\lambda} ) constraints on repeat similarity. Fig.  2 shows graphically some of the information contained in this table.

family	$2L$	$S_{\mathrm{rand}}$	$S_1$	$S_2$	$S_{\mathrm{full}}$	$S_{\mathrm{ir}}$	$S_{\mathrm{ir},\lambda }$
ANK	66	290	181 $\pm$ 0.05	169.7 $\pm$ 0.6	167.2 $\pm$ 0.3	176.7 $\pm$ 0.1	172 $\pm$ 0.4
LRR	48	211	130 $\pm$ 0.05	114 $\pm$ 0.4	113.2 $\pm$ 0.3	123.1 $\pm$ 0.1	118.8 $\pm$ 0.1
TPR	68	299	169 $\pm$ 0.1	145.4 $\pm$ 0.7	141.4 $\pm$ 0.3	157.6 $\pm$ 0.1	146.9 $\pm$ 0.4

Equations

$$P(\mbox{\boldmath$\sigma$})=(1/Z)e^{-E(\mbox{\boldmath$\sigma$})},$$

$$E(\mbox{\boldmath$\sigma$})=-\sum_{i=1}^{2L}h_{i}(\sigma_{i})-\sum_{i,j=1}^{2L}J_{ij}(\sigma_{i},\sigma_{j})+\lambda_{{\rm ID}}(\mbox{\boldmath$\sigma$})\,,$$

$$P(\mbox{\boldmath$\sigma$})=(1/Z)e^{2NF(\mbox{\boldmath$\sigma$})}.$$

$$P(\mbox{\boldmath$\sigma$})=(1/Z)e^{2Nf(\Delta G(\mbox{\boldmath$\sigma$}))}.$$

$$E(\mbox{\boldmath$\sigma$})=-2Nf(\Delta G(\mbox{\boldmath$\sigma$})).$$

$$P_{\rm sel}(\mbox{\boldmath$\sigma$})=(1/Z)e^{\Delta G(\mbox{\boldmath$\sigma$})/T_{\rm sel}},$$

$$⟨\Delta G{⟩}_{T_{\mathrm{sel}}}=\Delta G_{\mathrm{sel}}.$$

$$E(\mbox{\boldmath$\sigma$})=E_{0}-\Delta G(\mbox{\boldmath$\sigma$})/T_{\rm sel},$$

$$S=\ln\mathcal{N}(\mbox{\boldmath$\sigma$}:\Delta G(\mbox{\boldmath$\sigma$})>\Delta G_{\rm sel}).$$

$$S=-\sum_{\mbox{\boldmath$\sigma$}}P(\mbox{\boldmath$\sigma$})\ln P(\mbox{\boldmath$\sigma$}),$$

$$h_i({\sigma }_i{)}^{t+1}\leftarrow h_i({\sigma }_i{)}^t+{ϵ}_m[f_i({\sigma }_i)-f_i^{\mathrm{model}}({\sigma }_i)],$$

$$\lambda_{\rm ID}(\mbox{\boldmath$\sigma$})^{t+1}\leftarrow\lambda_{\rm ID}(\mbox{\boldmath$\sigma$})^{t}-\epsilon_{\rm ID}[P(\rm ID(\mbox{\boldmath$\sigma$}))-P(\rm ID(\mbox{\boldmath$\sigma$}))^{\rm model}].$$

$$J_{ij}({\sigma }_i,{\sigma }_j{)}^{t+1}\leftarrow 0.$$

$$\begin{array}{rl}{J}_{ij}({\sigma }_i,{\sigma }_j{)}^{t+1}\leftarrow {ϵ}_j[& f_{ij}({\sigma }_i,{\sigma }_j)-f_{ij}^{\mathrm{model}}({\sigma }_i,{\sigma }_j)-\\ & \gamma \mathrm{sign}(f_{ij}({\sigma }_i,{\sigma }_j)-f_{ij}^{\mathrm{model}}({\sigma }_i,{\sigma }_j))].\end{array}$$

$$\begin{array}{rl}{J}_{ij}({\sigma }_i,{\sigma }_j{)}^{t+1}\leftarrow J_{ij}({\sigma }_i,{\sigma }_j{)}^t+{ϵ}_j[& f_{ij}({\sigma }_i,{\sigma }_j)-f_{ij}^{\mathrm{model}}({\sigma }_i,{\sigma }_j)-\\ & \gamma \mathrm{sign}(J_{ij}({\sigma }_i,{\sigma }_j{)}^t)].\end{array}$$

$$J_{ij}({\sigma }_i,{\sigma }_j{)}^{t+1}\leftarrow 0.$$

$$f^{\prime}_{i}(\sigma_{i})=f^{\prime}_{i+L}(\sigma_{i})=\frac{1}{2}\Big{(}f_{i}(\sigma_{i})+f_{i+L}(\sigma_{i})\Big{)},$$

$$f_{ij}^{\prime }({\sigma }_i,{\sigma }_j)$$

$$S$$

$$\displaystyle=\sum_{\mbox{\boldmath$\sigma$}}p(\mbox{\boldmath$\sigma$})E(p(\mbox{\boldmath$\sigma$}))+\log{Z}$$

$$\displaystyle=-\sum_{\mbox{\boldmath$\sigma$}}p(\mbox{\boldmath$\sigma$})\log{p(\mbox{\boldmath$\sigma$})}\,,$$

$$E_{\alpha}(\mbox{\boldmath$\sigma$})=-\sum_{i}h_{i}(\sigma_{i})+\alpha\Big{[}-\sum_{ij}J_{ij}(\sigma_{i},\sigma_{j})+\lambda_{{\rm ID}}\Big{]}\,,$$

$$F(1)=F(0)+{\int }_0^1\mathrm{d}\alpha {⟨-\sum _{ij}{J}_{ij}({\sigma }_i,{\sigma }_j)+{\lambda }_{\mathrm{ID}}⟩}_{\alpha },$$

$$F(0)=-\sum _i\log \sum _{{\sigma }_i}{\mathrm{e}}^{h_i({\sigma }_i)},$$

$$S=⟨E⟩+\sum _i\log \sum _{{\sigma }_i}{\mathrm{e}}^{h_i({\sigma }_i)}-{\int }_0^1\mathrm{d}\alpha {⟨-\sum _{ij}{J}_{ij}({\sigma }_i,{\sigma }_j)+{\lambda }_{\mathrm{ID}}⟩}_{\alpha }.$$

$$\begin{split}{\rm Var}(S_{indep})=&\frac{1}{N_{s}}\Big{[}\sum_{i}\sum_{\sigma_{i}}p_{i}(\sigma_{i})\log{p_{i}(\sigma_{i})}^{2}-S_{indep}^{2}\Big{]}\\ &+O(\frac{1}{N_{s}^{2}})\,.\end{split}$$

$$D_{\mathrm{KL}}(ANK||TPR)=⟨E_{\mathrm{TPR}}-E_{\mathrm{ANK}}{⟩}_{\mathrm{ANK}}+F_{\mathrm{ANK}}-F_{\mathrm{TPR}},$$

$$D_{\mathrm{KL}}(TPR||ANK)=⟨E_{\mathrm{ANK}}-E_{\mathrm{TPR}}{⟩}_{\mathrm{TPR}}+F_{\mathrm{TPR}}-F_{\mathrm{ANK}},$$

Full text

I 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 Dryden2008 , even more so when specializing to a given structure Shakhnovich1998 . Clearly, not all variants are equally likely to survive Salisbury1969 ; Mandecki1998 ; Dokholyan2002 . 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 bateman2004pfam ; finn2013pfam ). More interesting information can be extracted from covariation of amino acid usages at pairs of positions neher1994frequent ; pmid24573474 ; Szurmant2018 or using machine-learning techniques Tubiana2018 . Models of protein sequences based of pairwise covariations have been shown to successfully predict pair-wise amino-acid contacts in three dimensional structures weigt2009pnas ; morcos2011pnas ; aurell2013pre ; hopf2012three ; espada2015BMC ; Espada2018 , aid protein folding algorithms Schug2009 ; marks2011protein , and predict the effect of point mutations Espada2018 ; Contini ; Haldane2016 ; weigt2015mbe . 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  MaynardSmith1970 ; Weinreich2006 .

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-$\kappa$B, p16, Notch pmid17176038 ). Usually they are composed of tandem repetitions of $\sim 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.

A recent study Tian2017 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 $\Delta 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 Barton2016a . 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.

II Results

II.1 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(\mbox{\boldmath\sigma}) to find in the family of interest a particular sequence \mbox{\boldmath\sigma}=(\sigma_{1},\ldots,\sigma_{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(\mbox{\boldmath\sigma}) is obtained as the distribution of maximum entropy Jaynes 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 {\rm ID}(\mbox{\boldmath\sigma})=\sum_{i=1}^{L}\delta_{\sigma_{i},\sigma_{i+L}} between consecutive repeats. Using the technique of Lagrange multipliers, the distribution can be shown to take the form Espada2018 :

[TABLE]

with

[TABLE]

where $h_{i}(\sigma)$, $J_{ij}(\sigma_{i},\sigma_{j})$, and $\{\lambda_{\rm ID}\}$, ${\rm ID}=0,1,\ldots,L$, are adjustable Lagrange multipliers that are fit to the data to reproduce the experimentally observed site-dependent amino-acid frequencies $f_{i}(\sigma_{i})$, joint probabilities between two positions, $f_{ij}(\sigma_{i},\sigma_{j})$, and the distribution of Hamming distances between consecutive repeats P({\rm ID}(\mbox{\boldmath\sigma})), 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 $f_{i}(\sigma_{i})-f_{i}^{\rm model}(\sigma_{i})$, $f_{ij}(\sigma_{i},\sigma_{j})-f_{ij}^{\rm model}(\sigma_{i},\sigma_{j})$ and P({\rm ID}(\mbox{\boldmath\sigma}))-P({\rm ID}(\mbox{\boldmath\sigma}))^{\rm 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. Espada2018 . See Sec. IV.2 for more details. We tested the convergence of the model learning by synthetically generating datasets and relearning the model (see Sec. IV.5).

By analogy with Boltmzan’s law, we call E(\mbox{\boldmath\sigma}) 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 E_{\rm full}(\mbox{\boldmath\sigma}). 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 $\lambda_{{\rm ID}}$ term in Eq. 2. We call the corresponding energy function $E_{\rm 2}$. 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 $h_{i}=h_{i+L}$ ($E_{\rm ir}$), or only linked through phylogenic conservation through $\lambda_{\rm ID}$ ($E_{\rm ir,\lambda}$). Finally one can remove all interaction constraints to make all positions independent of each other ($E_{\rm 1}$), or even remove all constraints ($E_{\rm rand}\equiv 0$).

II.2 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(\mbox{\boldmath\sigma}). 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(\mbox{\boldmath\sigma}) to predict the effect of point mutations on the folding stability measured by the free energy difference between the folded and unfolded states, $\Delta G$, called the unfolding energy Espada2018 ; Contini . Synthetic sequences with low E(\mbox{\boldmath\sigma}) have also been shown to reproduce the fold and function of natural sequences Tian2018 . Here, extending an argument already developed in previous work Shakhnovich1993 ; Shakhnovich1993b ; Dokholyan2001 ; Morcos2014 , 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 $\sigma$, i.e. the probability of finding it in a population, is related to its fitness F(\mbox{\boldmath\sigma}). In the limit where mutations affecting the protein are rare compared to the time it takes for mutations to spread through the population, Kimura Kimura1962 showed that the probability of a mutation giving a fitness advantage (or disadvantage depending on the sign) $\Delta F$ over its ancestor will fix in the population with probability $2\Delta F/(1-e^{-2N\Delta F})$, where $N$ is the effective population size. The dynamics of successful substitution satisfies detailed balance Berg2004 , with the steady state probability

[TABLE]

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 $\Delta G$, F(\mbox{\boldmath\sigma})=f(\Delta G(\mbox{\boldmath\sigma})), then the distribution of genotypes we expect to observe in a population is

[TABLE]

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(\mbox{\boldmath\sigma}) and F(\mbox{\boldmath\sigma}) 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 $\Delta G$:

[TABLE]

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

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(\Delta G)$ is a threshold function, $f(\Delta G)=0$ for $\Delta G>\Delta G_{\rm sel}$ and $-\infty$ otherwise (lethal). In that case the probability distribution is P(\mbox{\boldmath\sigma})=(1/Z)\Theta(\Delta G-\Delta G_{\rm sel}), where $\Theta(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 $\Delta G_{\rm sel}$, and is equivalent to a “canonical” description Shakhnovich1993 ; Shakhnovich1993b ; Morcos2014 :

[TABLE]

where the “temperature” $T_{\rm sel}$ is set to match the mean $\Delta G$ between the two descriptions:

[TABLE]

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):

[TABLE]

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 $\langle E\rangle$.

II.3 Equivalence between two definitions of entropies

There are several ways to define the diversity of a protein family. The most intuitive one, followed by Tian2017 , is to count the total number of amino acid sequences that have an unfolding free energy $\Delta G_{\rm sel}$ above a threshold $\Delta G_{\rm sel}$ Shakhnovich1998 . This number naturally defines a Boltzmann entropy,

[TABLE]

Alternatively, starting from a statistical model P(\mbox{\boldmath\sigma}), one can calculate its Shannon entropy, defined as

[TABLE]

as was done in Ref. Barton2016a . 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 $E_{\rm sel}$ in terms of statistical energy. This threshold is given by $E_{\rm sel}=E_{0}-\Delta G_{\rm sel}/T_{\rm sel}$, where $E_{0}$ and $T_{\rm sel}$ can be obtained directly by fitting (Eq. 8) to single-mutant experiments. $E_{\rm sel}$ can also be obtained as a discrimination threshold separating sequences that are known to fold properly versus sequences that do not Tian2017 . In that case, assuming that the linear relationship (Eq. 8) was evaluated empirically using single mutants, this relationship can be inverted to get $\Delta G_{\rm sel}$ 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. $E_{\rm sel}=\langle E\rangle$, where the average is performed using the distribution defined in Eqs. 1-2.

II.4 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: $E_{\rm rand}$, $E_{1}$, $E_{\rm ir}$, $E_{\rm ir,\lambda}$, $E_{2}$, $E_{\rm full}$. 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 Barton2016a , for both parameter fitting and entropy estimation. Here we estimated the entropies using thermodynamic integration of Monte-Carlo simulations, as detailed in Sec. IV.4. 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 ($\sim 40-50\%$) compared to random polypeptide string models of the same length $2L$ (of entropy $S_{\rm rand}=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: $(S_{1}-S_{\rm full})/S_{\rm full}=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 $e^{\rm S_{1}}\sim 3\cdot 10^{54}$ to $e^{S_{\rm full}}\sim 2\cdot 10^{50}$ for ANK, from $10^{39}$ to $10^{34}$ for LRR, and from $7\cdot 10^{50}$ to $4\cdot 10^{42}$ for TPR.

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 $E_{1}$ to $E_{\rm ir}$ to $E_{\rm ir,\lambda}$ to $E_{\rm full}$. This division allows us to quantify the relative importance of phylogenic similarity between consecutive repeats ($\lambda_{\rm ID}$) relative to the impact of functional interactions ($J_{ij}$), 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 $S_{\rm ir}-S_{\rm ir,\lambda}=4.5$ bits for ANK, 4.3 bits for LRR, and 10.7 bits for TPR. The contribution of repeat-repeat interactions ($S_{\rm ir,\lambda}-S_{\rm full}\sim 5$ bits for all three families) is comparable or of the same order of magnitude as that of within-repeat interactions ($S_{1}-S_{\rm ir}=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 $S_{\rm ir}-S_{\rm ir,\lambda}-S_{2}+S_{\rm full}>0$ (Fig. 2, orange bars). This redundancy comes from the fact that, in absence of an explicit constaint on $P({\rm ID})$ in $E_{2}$, the interaction couplings $J_{i,i+L}(\sigma,\sigma)$ between homologous positions in the two repeats is expected to favor pairs of identical residues to mimic the effect of $\lambda_{\rm 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.

II.5 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 $J_{ij}$ was allowed to be non-zero only within a certain interaction range $d(i,j)\leq 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 Fig. S1 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 Fig. S2 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.

II.6 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 $\sigma_{i}$ can only take two values, have been extensively studied in the context of spin glasses Mezard . In these systems, the interaction terms $-J_{ij}(\sigma_{i},\sigma_{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 Bryngelson7524 .

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 $E_{\rm full}$ (learned with all constraints including on $P(\rm ID)$, but taken with $\lambda_{\rm 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 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 ($E_{\rm full}$), and obtained very similar (see Fig. S6 for ANK). When ranked from largest to smallest, the distribution of basin sizes follows a power law (Fig. 5A for ANK and Fig. S3A and Fig. S4A 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.

The partition of sequences into basins allows for the definition of a new kind of entropy $S_{\rm conf}=-\sum_{b}P(b)\ln P(b)$ called configurational entropy, based on the distribution of basin sizes, P(b)=\sum_{\mbox{\boldmath\sigma}\in b}P(\mbox{\boldmath\sigma}), where \mbox{\boldmath\sigma}\in b means that energy minimization starting with sequence $\sigma$ leads to basin $b$. This configurational entropy measures the effective diversity of basins, and is thus much lower than the sequence entropy $S_{\rm full}$, while the difference $S_{\rm full}-S_{\rm conf}$ measures the average diversity of sequences within each basin. We find $S_{\rm conf}=$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 Dokholyan2001 , 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 Figs. 5, S3, and S4), suggesting a hierarchical structure of basins, as is common in spin glasses Mezard .

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, $E_{\rm ir}$. 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 Figs. 5, S3, and S4. 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 $E_{1}$ defines a convex energy landscape with a single local minimum — the consensus sequence — as all constraints $h_{i}$ 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. 3 A). We find that the more interactions we add, the more local minima we find (Fig. S5A and B 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 in entropy observed in 3 A.

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.

II.7 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(\mbox{\boldmath\sigma}) (given by Eqs. 1-2) of two different families, after aligning them together in a single multiple sequence alignment (see Sec. IV.7).

We find essentially no similarity between ANK and TPR, despite them having similar lengths: $D_{\rm KL}({\rm ANK}||{\rm TPR})=227.6$ bits, and $D_{\rm KL}({\rm TPR}||{\rm ANK})=214.1$ bits. These values are larger than the Kullback-Leibler divergence between the full models for these families and a random polypetide, $D_{\rm KL}({\rm ANK}||{\rm rand})=122.8$ bits, and $D_{\rm KL}({\rm TPR}||{\rm 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.

III 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 $10^{36}$ to $10^{50}$ sequences) are of course huge compared to the number of sequences in our initial samples ($\sim 20,500$ for ANK, $\sim 18,800$ for LRR, and $\sim 10,000$ for TPR), but comparable to the total number of proteins having been explored over the whole span of evolution, estimated to be $10^{43}$ in Ref. Dryden2008 .

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 ($E_{1}$) 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 Dokholyan2001 .

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  pmid1749933 . 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 Facco2019 , 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 Figs. 5, S3, and S4), 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 pmid21715155 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 pmid18396879 ; Barrick2008 . 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.

IV Methods

IV.1 Data curation

We use a previously curated alignment of pairs of repeats for each family  Espada2018 : 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 bateman2004pfam ; finn2013pfam . In order to improve the data obtained from the PFAM database, we used original full protein sequences available in UniProt database TheUniProtConsortium01012014 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 $f_{i}(\sigma_{i})$ of an amino acid $\sigma_{i}$ at position $i$ and the joint frequency $f_{ij}(\sigma_{i},\sigma_{j})$ of two amino acids $\sigma_{i}$ and $\sigma_{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 cdhit2002 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({\rm ID}(\mbox{\boldmath\sigma})), with {\rm ID}(\mbox{\boldmath\sigma})=\sum_{i=1}^{L}\delta_{\sigma_{i},\sigma_{i+L}}.

IV.2 Model fitting

In order to obtain a model that reproduces the experimentally observed site-dependent amino-acid frequencies, $f_{i}(\sigma_{i})$, correlations between two positions, $f_{ij}(\sigma_{i},\sigma_{j})$, and the distribution of Hamming distances between consecutive repeats, P({\rm ID}(\mbox{\boldmath\sigma})), we apply a likelihood gradient ascent procedure, starting from an initial guess of the $h_{i}(\sigma_{i})$, $J_{ij}(\sigma_{i},\sigma_{j})$ and \lambda_{\rm ID}(\mbox{\boldmath\sigma}) 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^{-\Delta E}$, where $\Delta 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 $f^{\rm model}_{i}(\sigma_{i})$ and $f^{\rm model}_{ij}(\sigma_{i},\sigma_{j})$, as well as P^{\rm model}({\rm ID}(\mbox{\boldmath\sigma})), and update the parameters of Eq. 2 following the gradient of the likelihood. The local field and \lambda_{\rm ID}(\mbox{\boldmath\sigma}) are updated along the gradient of the per-sequence log-likelihood, equal to the difference between model and data averages:

[TABLE]

[TABLE]

As the number of parameters for the interaction terms $J_{ij}$ is large ($=21^{2}L^{2}$), we force to 0 those that are not contributing significantly to the model frequencies through a $L_{1}$ regularisation $\gamma\sum_{ij,\sigma,\tau}|J_{ij}(\sigma,\tau)|$ added to the likelihood. This leads to the following rules of maximization: If $J_{ij}(\sigma_{i},\sigma_{j})^{t}=0$ and $|f_{ij}(\sigma_{i},\sigma_{j})-f_{ij}^{\rm model}(\sigma_{i},\sigma_{j})|<\gamma$

[TABLE]

If $J_{ij}(\sigma_{i},\sigma_{j})^{t}=0$ and $|f_{ij}(\sigma_{i},\sigma_{j})-f_{ij}^{\rm model}(\sigma_{i},\sigma_{j})|>\gamma$

[TABLE]

If \Big{[}J_{ij}(\sigma_{i},\sigma_{j})^{t}+\epsilon_{j}[f_{ij}(\sigma_{i},\sigma_{j})-f_{ij}^{\rm model}(\sigma_{i},\sigma_{j})-\gamma{\rm sign}(J_{ij}(\sigma_{i},\sigma_{j})^{t})]\Big{]}J_{ij}(\sigma_{i},\sigma_{j})^{t}\geq 0

[TABLE]

If \Big{[}J_{ij}(\sigma_{i},\sigma_{j})^{t}+\epsilon_{j}[f_{ij}(\sigma_{i},\sigma_{j})-f_{ij}^{\rm model}(\sigma_{i},\sigma_{j})-\gamma{\rm sign}(J_{ij}(\sigma_{i},\sigma_{j})^{t})]\Big{]}J_{ij}(\sigma_{i},\sigma_{j})^{t}<0

[TABLE]

The optimization parameters were set to: $\epsilon_{m}=0.1$, $\epsilon_{j}=0.05$, $\epsilon_{\rm ID}=10$, and $\gamma=0.001$.

To estimate the model error, we compute $f_{i}(\sigma_{i})-f_{i}^{\rm model}(\sigma_{i})$ and $f_{ij}(\sigma_{i},\sigma_{j})-f_{ij}^{\rm model}(\sigma_{i},\sigma_{j})$. 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 $\lambda_{\rm ID}$: 5(P(\rm ID(\mbox{\boldmath\sigma}))-P(\rm ID(\mbox{\boldmath\sigma}))^{\rm model}). We repeat the procedure above until the maximum of all errors, $|f_{i}(\sigma_{i})-f_{i}^{\rm model}(\sigma_{i})|$, $|f_{ij}(\sigma_{i},\sigma_{j})-f_{ij}^{\rm model}(\sigma_{i},\sigma_{j})|$ and 5|P(\rm ID(\mbox{\boldmath\sigma}))-P(\rm ID(\mbox{\boldmath\sigma}))^{\rm model}|, goes below $0.02$, as in Ref. Espada2018 .

IV.3 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. 3 A. We start from the independent model $h_{i}(\sigma_{i})=\log{f_{i}(\sigma_{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|\leq 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 $\lambda_{\rm 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 peng2002pnas ; pluckthun2003jmb , 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 $\lambda_{\rm ID}$ ($E_{2}$) we remove $\lambda_{\rm 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 $\lambda$ ($E_{\rm ir}$ and $E_{\rm ir,\lambda}$, 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 $J_{ij}$ terms between different repeats. We also impose that the $h_{i}$ fields and intra-repeats $J_{ij}$ 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 $f_{i}(\sigma_{i})$ and $f_{ij}(\sigma_{i},\sigma_{j})$, such that

[TABLE]

and

[TABLE]

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.

IV.4 Entropy estimation

In practice to calculate the entropy $S$ of the protein families we relate it to the internal energy E=-\log{p(\mbox{\boldmath\sigma})} and the free energy $F=-\log{Z}$:

[TABLE]

We generate sequences according to the energy function in Eq. 2 and use them to numerically compute $\langle E\rangle$. To calculate the free energy we use the auxilliary energy function:

[TABLE]

where the interaction strength across different sites can be tuned through a parameter $\alpha$ that is changed from [math] to $1$. We generate protein sequence ensembles with different values of $\alpha$ and use them to calculate $F$ as a function of $\alpha$, $F(1)=F(0)+\int_{0}^{1}{\rm d}\alpha\frac{{\rm d}F}{{\rm d}\alpha}$:

[TABLE]

where the average over $\alpha$ is taken over the sequences generated with a certain value of $\alpha$, characterized by the ensemble with probability p_{\alpha}(\mbox{\boldmath\sigma})=(1/Z_{\alpha})e^{-E_{\alpha}(\mbox{\boldmath\sigma})}. $F(0)$ is the free energy for an independent sites model:

[TABLE]

where the first sum is taken over protein sites and the second over all possible amino-acids at a given site. Eq. 22 and Eq. IV.4 result in the thermodynamic sampling approximation for calculating the entropy Frankel2007 :

[TABLE]

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

IV.5 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 ${\rm Var}(f_{i}(\sigma_{i}))={\left(p_{i}(\sigma_{i})(1-p_{i}(\sigma_{i}))\right)}/{N_{s}}$ and the covariance between the frequencies of different amino-acids $\sigma$ and $\sigma^{\prime}$ at the same site $i$ is ${\rm Cov}(f_{i}(\sigma_{i}),f_{i}(\sigma^{\prime}_{i}))=-{\left(p_{i}(\sigma_{i})p_{i}(\sigma^{\prime}_{i})\right)}/{N_{s}}$ where $N_{s}$ is the sample size and $p_{i}(\sigma_{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}/{N_{s}}$:

[TABLE]

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 ${1}/{\sqrt{N_{s}}}$ with $N_{s}=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 in Section IV.2 and the initial condition that we chose to start learning from. Fig. S1 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. Fig. S1 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 Fig. 3,Fig. 2, Table 1, Fig. S2.

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. Fig. S1 shows the mean entropy of ANK as in Fig. 3 A 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 IV.2). 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 in the Section IV.2, 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 such as LRR and TPR reaching similar conclusions (Fig. S2).

IV.6 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 $E_{\rm full}$ yields very similar results (see Fig. S6 for ANK)

We perform this energy landscape exploration learning the parameters of the Hamiltonian in Eq. 2 (refer to Section IV.2 for the learning procedure), and then set $\lambda_{\rm 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. 3 A. That is Eq. 2, in which $J_{ij}$ 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 IV.2. 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 Fig. S3 and Fig. S4 and lead to similar conclusions as for the ANK family.

IV.7 Kullback-Leibler divergence

The Kullback-Leibler divergence between two families ${\rm A}$ and ${\rm B}$ is defined as D_{\rm KL}({\rm A}||{\rm B})=\sum_{\mbox{\boldmath\sigma}}p_{\rm A}(\mbox{\boldmath\sigma})\log_{2}{p_{\rm B}(\mbox{\boldmath\sigma})/p_{\rm A}(\mbox{\boldmath\sigma})}. We can substitute the sequence ensembles for ANK and TPR in the definition of the probabilities obtaining:

[TABLE]

[TABLE]

where the notation $\langle\rangle_{\rm ANK}$ means that the average is calculated over sequences drawn from the ANK ensemble: P(\mbox{\boldmath\sigma})_{\rm ANK}=(1/Z_{\rm ANK})e^{-E(\mbox{\boldmath\sigma})_{\rm ANK}}. Therefore $\langle E_{\rm TPR}\rangle_{\rm 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 $\langle E_{\rm ANK}\rangle_{\rm TPR}$). The terms $F_{\rm ANK}$ and $F_{\rm TPR}$ are calculated in the same way as when estimating the entropy through Eqs. (21),(22), as explained in Section IV.4.

For the control against a random polypeptide of length $L$ we use $D_{\rm KL}({\rm FAM}||{\rm rand})=\log{\Lambda}-S({\rm FAM})$, where $\Lambda=21^{L}$ is the total number of possible sequences of length $L$.

**Acknowledgements. ** This work was partly supported by ERC CoG 724208 and ECOS Sud - MINCyT no. A14E04