Parsimonious evolutionary scenario for the origin of allostery and coevolution patterns in proteins

TL;DR

This paper proposes an evolutionary scenario where proteins develop allostery and coevolution patterns as by-products of selection for discriminative ligand binding, driven by conformational changes and physical constraints.

Contribution

It introduces a parsimonious evolutionary model linking allostery, coevolution, and ligand discrimination, supported by a physical model and testable experimental predictions.

Findings

Proteins' allostery and coevolution patterns can arise as evolutionary spandrels.

Conformational changes are essential for discriminative ligand binding.

Most coevolving residues are related to thermal stability, not direct contact.

Abstract

Proteins display generic properties that are challenging to explain by direct selection, notably allostery, the capacity to be regulated through long-range effects, and evolvability, the capacity to adapt to new selective pressures. An evolutionary scenario is proposed where proteins acquire these two features indirectly as a by-product of their selection for a more fundamental property, exquisite discrimination, the capacity to bind discriminatively very similar ligands. Achieving this task is shown to typically require proteins to undergo a conformational change. We argue that physical and evolutionary constraints impel this change to be controlled by a group of sites extending from the binding site. Proteins can thus acquire a latent potential for allosteric regulation and evolutionary adaptation because of long-range effects that initially arise as evolutionary spandrels. This scenario accounts for the groups of conserved and coevolving residues observed in multiple sequence alignments. However, we propose that most pairs of coevolving and contacting residues inferred from such alignments have a different origin, related to thermal stability. A physical model is presented that illustrates this evolutionary scenario and its implications. The scenario can be implemented in experiments of protein evolution to directly test its predictions.

Figures

A). While conformational proofreading is particularly relevant to proteins involved in search processes [vlaminck ](#bib.bib40), exquisite discrimination may be more relevant to enzymes, which have to bind to a reactant but release its product. Despite differences, the two mechanisms imply a similar notion of fine-tuning, which we propose here to be a sufficient factor for evolving an extended set of coupled sites.

Equations

$$U(x,a,\ell )=-\sum _{⟨i,j⟩}{J}_{ij}(a_i,a_j)x_i{x}_j-\sum _i{h}_i(a_i)x_i-\ell x_b$$

$$\varphi (a)=\min (-\Delta F_r(a),\Delta F_w(a)),$$

$${\omega }_{i,r}(a)=\frac{1}{2}\left|⟨x_i{⟩}_{a,{\ell }_r}-⟨x_i{⟩}_{a,{\ell }_0}\right|,{\bar{\omega }}_r(a)=\frac{1}{N}\sum _{i=1}^N{\omega }_{i,r}(a)$$

$$D_i=\sum _{a_i=1}^q{f}_i(a_i)\ln \frac{f_i(a_i)}{q^{-1}},$$

$$p_k(a)=\frac{e^{-\beta (\Delta F_k(a)-{\mu }_k)}}{1+e^{-\beta (\Delta F_r(a)-{\mu }_r)}+e^{-\beta (\Delta F_w(a)-{\mu }_w)}},$$

$$U(x,-h)=U(-x,h).$$

$$U(x,h)=\frac{1}{2}k(|x|-r{)}^2-hx,x\in \mathbb{R},$$

$$F(h)=-\frac{1}{2}\frac{h^2}{k}-|h|r,\hat{x}(h)=\frac{h}{k}+\mathrm{sign}(h)r.$$

$$U(x,h)=k{\left(\sqrt{x^2+d^2}-r\right)}^2-hx,x\in \mathbb{R},$$

$$\hat{x}(h)\left(1-\frac{r}{\sqrt{\hat{x}(h{)}^2+d^2}}\right)=\frac{h}{2k}.$$

$$U(x,h)=\frac{1}{2}k(x-r{)}^2-hx,x\in \mathbb{R}.$$

$$F(h)=-\frac{h^2}{2k}-hr+\frac{1}{2\beta }\ln \frac{\beta k}{2\pi }.$$

$$U(x,h)=-hx,x\in \{-1,+1\},$$

$$F(h)=-{\beta }^{-1}\ln \left(e^{\beta h}+e^{-\beta h}\right).$$

Full text

Proteins combine a capacity to perform exquisite tasks such as specific binding and catalysis with a capacity to adapt, both on physiological time scales through allostery and on evolutionary time scales through mutations. What principles underly this duality? How does it originate in the course of evolution? For proteins folding into a stable structure, a few generic properties constrain the explanations that we may provide. One is the ubiquity of long-range effects experimentally observed in all studied proteins: perturbations at a distance from the active site, whether in the form of mutations or of binding to another molecule, can affect significantly protein function licata1995long ; Daugherty:2000gz ; Morley:2005kc . Another is the coexistence of two types of patterns of coevolution revealed by statistical analyses of alignments of protein sequences: structurally connected groups of conserved and coevolving residues called sectors Lockless:1999uf ; Halabi:2009jc , and structurally contacting pairs of residues distributed across the structure Morcos:2011jg . Finally, another feature reported to be essential to protein machineries is their capacity to undergo conformational changes fersht1985enzyme .

Here, we propose that a key to understand the link between high-performance and adaptability in proteins is to understand the physical and evolutionary implications of one of their most fundamental requirements: discriminative binding. Discriminative binding, the capacity to bind to a particular ligand but not to other similar ones, is central to the function of many if not most proteins, from signaling proteins that respond to particular inputs and avoid cross-talk pawson2000protein to enzymes that bind to the transition state of a reaction but release its product fersht1985enzyme , transcription factors that recognize specific promoters among a profusion of similar motifs todeschini2014transcription , or antibodies that are highly specific to particular antigens VanRegenmortel:1998uv .

As binding takes place at a particular location on a protein surface, it is not a priori expected to be sensitive to distant perturbations, whether in the form of mutations or interactions with other molecules. But achieving exquisite discrimination imposes particular constraints. Physically, we shall show how it typically involves a conformation change, or even a switch between two states that pre-exist any interaction with a ligand. Evolutionarily, only few sequences can achieve discrimination which is shown to require a finely tuned binding site. This tuning is generally difficult to accomplish based only on the few amino acids directly interacting with the ligand. Additional tuning knobs are effectively provided by sites coupled to the binding site. Since only few residues interact directly with it, this generally implies recruiting distant sites. In this scenario, binding specificity is thus controlled through a conformational change by a group of sites extending beyond the binding site. This sensitivity of a functional phenotype to multiple and possibly distant sites endow proteins with a latent capacity to evolve allosteric regulation and adapt to new selective pressures. When integrating an additional constraint, thermal stability, we shall show that it also accounts for the different patterns of intra-protein coevolution that multiple sequence alignments of natural proteins report.

We demonstrate this evolutionary scenario by means of simple physical models. These models are not meant to describe any particular protein but to capture the main constraints relevant to a discussion of long-range effects within stable protein folds: the fact that physical interactions are short-range. We analyze a variety of such physical models to ensure that our conclusions do not depend on the nature of the interactions: in all cases, we find that the physical implications of exquisite discrimination is a finely-tuned conformational switch between the ligand-free and the ligand-bound states. To examine the additional implications of evolutionary constraints, we then focus on one of the simplest models, a spin model where each site can be in just two states. This model can for instance be interpreted as describing side-chain fluctuations between a few discrete rotameric states in the context of a rigid backbone, one of the basic mechanisms by which long-range effects arise in proteins dubay2011long . We show that the evolutionary constraints cause with high probability a large part of the system to be involved in the discrimination. Sensitivity to perturbations distant from the active site thus arise as a by-product of a selection for discrimination, i.e., as an evolutionary spandrel spandrel . This side effect is shown to favor adaptation to new selective pressures when considering perturbations caused by mutations, and to enable the evolution of allostery when considering perturbations caused by interactions with other molecules. It is also shown to cause a generic coevolution pattern found in protein sequences, the extended groups of conserved and coevolving amino acids called sectors Halabi:2009jc . Another coevolution pattern is also generically found in protein sequences, which consists in isolated pairs of contacting amino acids distributed across the structure Morcos:2011jg . This pattern is shown to be captured by the model when introducing a constraint on thermal stability. The model thus provides an explanation for conformational switches, the origin of allostery, evolvability and the different patterns of coevolution observed in multiple sequence alignments as a consequence of only two selective constraints: discriminative binding and thermal stability. While it does not exclude other evolutionary scenarios that may lead to the same or different physical effects, it describes an interplay between evolutionary and physical constraints that is consistent with observations in many proteins and directly testable in experiments of protein evolution.

I From discriminative binding to conformational changes

Formally, protein evolution involves three types of variables: (i) physical degrees of freedom $x$ describing the conformation of the molecule, (ii) evolutionary degrees of freedom $a$ defined by the protein sequence, and (iii) environmental degrees of freedom $\ell$ defined by the surrounding medium, here restricted to either the solvent or a ligand interacting at a particular binding site. These three variables are linked in a potential $U(x,a,\ell)$ that dictates through Boltzmann’s law how the different physical conformations are sampled at thermal equilibrium: the probability of conformation $x$ is $P(x|a,\ell)=\exp[-\beta(U(x,a,\ell)-F(a,\ell))]$, where $F(a,\ell)=-\beta^{-1}\ln\int dx\exp[-\beta U(x,a,\ell)]$ represents the free energy of the system and $\beta$ the inverse temperature.

A problem of discrimination arises when two ligands $\ell_{r}$ and $\ell_{w}$ can potentially be substituted for the solvent $\ell_{0}$ at the binding site, but only $\ell_{r}$ is desirable. Finding a sequence $a$ that solves this problem generally involves minimizing $\Delta F_{r}(a)=F(a,\ell_{r})-F(a,\ell_{0})$, the binding free energy to the right ligand $\ell_{r}$, while maximizing $\Delta F_{w}(a)=F(a,\ell_{w})-F(a,\ell_{0})$, the binding free energy to the wrong ligand $\ell_{w}$. When $\ell_{r}$ and $\ell_{w}$ are similar, this may lead to a trade-off. A precise formulation of this trade-off depends on the concentrations of the two ligands as well as on the relative cost and benefit that binding to them entail. Here, we consider a strong form of discrimination and impose $F(a,\ell_{r})<F(a,\ell_{0})<F(a,\ell_{w})$. Under these conditions, binding with $\ell_{r}$ is favorable, $\Delta F_{r}(a)<0$, but binding with $\ell_{w}$ is unfavorable, $\Delta F_{w}(a)>0$. Different formulations of the same trade-off can also lead to the same conclusions (Appendix).

Satisfying $F(a,\ell_{r})<F(a,\ell_{0})<F(a,\ell_{w})$ when $\ell_{r}$ and $\ell_{w}$ are much more similar to each other than they are to $\ell_{0}$ is generally either impossible, or possible only for a restricted set of sequences $a$. This is best illustrated with a few elementary models. Consider first elastic networks, a mechanical framework commonly used as a coarse-grained description of proteins Sanejouand:2013kj . The simplest conceivable elastic network is a single mass attached to a fixed point by a spring. We assume here that it is additionally subject to a constant force $h=\ell-a$ controlled by the evolutionary and environmental degrees of freedom $a$ and $\ell$ (Fig. 1A). The potential is therefore $U(x,a,\ell)=k(|x|-r)^{2}/2-(\ell-a)x$ where $x$ is the position of the mass along a dimension to which it is confined, $k$ the stiffness of the spring and $r$ its equilibrium length. In this model, $\ell$ and $a$ may take arbitrary real values. When $\ell_{0}<\ell_{w}<\ell_{r}$, it is readily shown that satisfying $F(a,\ell_{r})<F(a,\ell_{0})<F(a,\ell_{w})$ requires $a$ to satisfy $(\ell_{0}+\ell_{w})/2<a<(\ell_{0}+\ell_{r})/2$ (Fig. 1A and Appendix). This corresponds to the mean conformation $\langle x\rangle_{a,\ell}=\int dxP(x|a,\ell)x$ switching from a negative to a positive value when the solvent $\ell_{0}$ is replaced by the ligand $\ell_{r}$. Importantly, this conformational switch is required only when $\ell_{w}$ and $\ell_{r}$ are sufficiently similar: if $\ell_{w}<\ell_{0}<\ell_{r}$, the solutions do not involve a change of sign (Appendix).

In this first model, achieving exquisite discrimination involves a switch between two states that are local minima of the potential $U(x,a,\ell_{0})$ (Fig. 1A). This, however, need not be the case as seen by considering an harmonic potential $U(x,a,\ell)=k(x-r)^{2}/2-(\ell-a)x$, which also requires $a$ to verify $(\ell_{0}+\ell_{w})/2<a<(\ell_{0}+\ell_{r})/2$ but cannot sustain multiple states (Appendix). In this model, a conformational change nevertheless occurs upon binding with amplitude $\Delta x=|\ell_{r}-\ell_{0}|/k$. A model with the same phenomenology can also be constructed from two springs (Fig. 1B). Varying one parameter in this model, we can continuously interpolate between a one-state and a two-state model (Fig. S1).

Exquisite discrimination can in principle be achieved without a conformational change, by a rigid lock-and-key mechanism, but fine-tuning the evolutionary parameters is in any case necessary (see shape space model in Appendix). A conformational change is, however, expected as soon as a minimal form of flexibility is present. A limiting case is a spin model where $x$ takes only two values, $x=\pm 1$. With $U(x,a,\ell)=(a-\ell)x$, the discrimination problem takes again the same form: when $\ell_{0}<\ell_{w}<\ell_{r}$, $a$ must be tuned to satisfy $(\ell_{0}+\ell_{w})/2<a<(\ell_{0}+\ell_{r})/2$ for the condition $F(a,\ell_{r})<F(a,\ell_{0})<F(a,\ell_{w})$ to be fulfilled, in which case binding induces a conformational switch from $\langle x\rangle_{a,\ell_{0}}<0$ to $\langle x\rangle_{a,\ell_{r}}>0$ (Fig. 1C and Appendix). Graphically, the conformational change is again linked to the need for $\ell_{0}$ and $\ell_{r}$ to be associated with different “branches” of the free energy ($h<0$ and $h>0$ for $F(h=a-\ell)$ in Fig. 1).

In summary, achieving exquisite discrimination requires a flexible system to be evolutionarily tuned and to physically change conformation upon binding. In proteins, however, the evolutionary degrees of freedom are not controlled by continuous variables but by a limited set of twenty amino acids. At the binding site, relatively few values of $a$ are thus available to tune $F(a,\ell)$. As shown below, a generic solution is to enlarge the evolutionary space by coupling the binding site to other sites. Interestingly, having multiple tuning knobs not only favors adaptation to a particular selective constraint but also to alternative constraints. Besides, if distant sites are thus involved, they are likely to be allosteric: a perturbation at those sites, such as an interaction with another molecule, will alter the binding properties. These implications of evolutionary and physical constraints are independent of the underlying mechanisms and may be demonstrated for a wide class of potentials Zadorin19 . We therefore illustrate them using a physical model with the simplest form of flexibility, a spin model.

II Two-dimensional spin model

As a simple model with a non-trivial geometry, we consider a spin model defined on a two-dimensional lattice with periodic boundary conditions along one dimension (Fig. 2A). This model is similar to the model introduced in Ref. Hemery:2015ei with two differences. First, for simplicity, we consider binary variables (spins) rather than continuous variables. Second and most importantly, the evolutionary scenario that we consider is totally different: instead of selecting explicitly for a long-range effect by varying the boundary conditions at the two extreme sides of the lattice, we select for binding specificity by varying the boundary conditions at a single site of the lattice. The point is indeed to show that long-range effects can evolve spontaneously with significant probability from a selection for exquisite discrimination, even though this selection is very localized.

Each node $i$ is associated with two variables, a physical variable $x_{i}$ that can take two values $x_{i}=\pm 1$, and an evolutionary variable $a_{i}$ that can take $q$ values $a_{i}=1,\dots,q$. The evolutionary variables determine the fields $h_{i}(a_{i})$ and couplings $J_{ij}(a_{i},a_{j})$ to which the physical variables are subject (Fig. 2B). Finally, a site $b$ is chosen to represent the binding site, which is subject to an additional external field $\ell$ representing an interaction with a ligand or with the solvent (red node in Fig. 2A). The potential is

[TABLE]

where $x\in\{-1,+1\}^{N}$, $a\in\{1,\dots,q\}^{N}$ and $\ell\in\mathbb{R}$, $N$ being the total number of sites. The couplings $J_{ij}(a_{i},a_{j})$ are possibly non-zero only between nearest neighbors in the two-dimensional lattice shown in Fig. 2A, so as to reflect the short-range nature of physical interactions. The free energy given a sequence $a$ and a ligand $\ell$ is as usual $F(a,\ell)=-\beta^{-1}\ln\sum_{x}\exp(-\beta U(x,a,\ell))$ with $\beta$ the inverse temperature.

As a score of the extent to which a sequence $a$ discriminates the right ligand $\ell_{r}$ from the wrong ligand $\ell_{w}$, we consider

[TABLE]

where $\Delta F_{r}(a)=F(a,\ell_{r})-F(a,\ell_{0})$ and $\Delta F_{w}(a)=F(a,\ell_{w})-F(a,\ell_{0})$ are the binding free energies to the two possible ligands. With this formulation, $F(a,\ell_{r})<F(a,\ell_{0})<F(a,\ell_{w})$ is equivalent to $\phi(a)>0$. Other fitness functions implementing the same trade-off between $\Delta F_{r}(a)$ and $\Delta F_{w}(a)$ may also be used and lead to similar conclusions (Fig. S2). We sample sequences $a$ with probability $P(a)\propto e^{\gamma\phi(a)}$ using a Metropolis Monte Carlo algorithm metropolis53 , which corresponds to an evolutionary dynamics in the origin-fixation limit, with $\phi(a)$ interpreted as a fitness function and $\gamma$ as an effective population size McCandlish:2014us (see Appendix for details).

For illustration, we consider a cylinder with $L=10$ layers of couplings and $W=5$ nodes per layer (Fig. 1A) for a total of $N=W(L+1)=55$ sites, in the range of smallest folding protein domains. With this geometry, the free energy $F(a,\ell)$ and other thermodynamical quantities can be computed exactly by transfer matrices baxter82 . We take $q=5$, smaller than the number of natural amino acids but sufficient to generate a large evolutionary space of size $q^{N}\simeq 10^{38}$. To define $h_{i}(a_{i})$ and $J_{ij}(a_{i},a_{j})$, we draw their values at random from normal distributions, $h_{i}(a_{i})\sim\mathcal{N}(0,\sigma_{h}^{2})$ and $J_{ij}(a_{i},a_{j})\sim\mathcal{N}(0,\sigma_{J}^{2})$, independently for each $i,j$ and each $a_{i},a_{j}$. Fixing the energy scale by setting $\beta=1$, we take $\sigma_{h}=\sigma_{J}=3$ so that some fields and couplings may be large relative to $\beta^{-1}$, as in natural proteins where the strength of some physical interactions, e.g., covalent bounds and steric constraints, can significantly exceed the scale of thermal fluctuations. We take $\ell_{0}=0$, $\ell_{w}=2$ and $\ell_{r}=3$, a choice of parameters for which $\phi(a)>0$ is achievable with a single spin, provided the adequate evolutionary diversity is available. Lastly, we take $\gamma=100$ to sample near-optimal values of $\phi(a)$.

With these parameters, most evolutionary trajectories lead to $\phi(a)>0$ after $1000$ iterations ($>98\%$, see inset of Fig. 3A). To visualize and quantify conformational changes induced by substituting $\ell_{r}$ for $\ell_{0}$, we introduce

[TABLE]

where $\langle x_{i}\rangle_{a,\ell}=\int dxP(x|a,\ell)x_{i}$ stands for the mean value of $x_{i}$ given $a,\ell$. $\omega_{i,r}(a)$ quantifies the extent to which site $i$ undergoes a conformational change upon binding to the right ligand: $\omega_{i,r}(a)=0$ indicates no change, while $\omega_{i,r}(a)=1$ indicates a maximal change between two polarized states. The site-averaged quantity $\bar{\omega}_{r}(a)$ may be interpreted as an effective fraction of sites taking part in the switch ($\bar{\omega}_{w}(a)$ is similarly defined by considering $\ell_{w}$ instead of $\ell_{r}$). As shown in Fig. 3A, $\bar{\omega}_{r}(a)$ varies widely from one evolved system to the next, even among systems with nearly identical fitness value $\phi(a)$. In most but not all cases, the set of sites involved in the switch extends beyond the binding site. The size and shape of this extension varies again from case to case (Fig. 4). As in the simpler models, a switch evolves only under a constraint for exquisite discrimination controlled by the similarity between the two ligands $\ell_{r}$ and $\ell_{w}$: if $\ell_{w}<\ell_{0}<\ell_{r}$ no switch evolves while if $\ell_{0}<\ell_{w}<\ell_{r}$ or $\ell_{0}<\ell_{r}<\ell_{w}$ it always does (Fig. S3B).

Consistent with the proposed scenario, sites involved in the switch, and only those sites, are sensitive to perturbations. This is verified by applying to the evolved systems an additional local field $h^{\prime}_{i}=\pm 2$ at each site $i$ and estimating the effect $\Delta_{h_{i}}\phi$ of this perturbation on the fitness $\phi(a)$: as seen in Fig. 4, the sites responding to this local perturbation are the same as those participating to the conformational change. An overlap with sites sensitive to mutations ($\Delta_{a_{i}}\phi$) is also evident although less straightforward because a mutation at $i$ affects the couplings $J_{ij}(a_{i},a_{j})$ to neighboring sites in addition to the local field $h_{i}(a_{i})$. In any case, a selective pressure for exquisite discrimination is sufficient in this model to generate with high probability distant allosteric sites and long-range mutational effects.

To assess the implications for adaptation, we define a measure of evolvability $\varepsilon(a)$ that reports the phenotypic diversity of sequences differing from $a$ by a single mutations. Here, the relevant phenotype is the set of ligands that a sequence can discriminate. In the limit where $\ell_{r}$ and $\ell_{w}$ are very similar, $|\ell_{r}-\ell_{w}|\to 0$, there is typically either no or a single value of $\ell_{r}$ for which $\phi(a)>0$. This value $\ell_{r}(a)$ provides a synthetic characterization of the phenotype of $a$. A measure $\varepsilon(a)$ of the phenotypic diversity of the neighborhood of $a$ can thus be defined from the number of different values that $\ell_{r}(a^{\prime})$ takes when considering all the single mutants $a^{\prime}$ of $a$. In practice, we fix an interval of phenotypes of interest $[\ell_{\rm min},\ell_{\rm max}]=[-1,4]$, partition it into small subintervals of length $\delta=0.1$, and for each single mutant $a^{\prime}$ of a given sequence $a$ find, if it exists, the interval to which $\ell_{r}(a^{\prime})$ belongs: $\varepsilon(a)$ is then defined as the fraction of subintervals covered by the $(q-1)W(L+1)$ single-point mutants of $a$. As seen in Fig. 3B, $\varepsilon(a)$ correlates with $\bar{\omega}_{r}(a)$, consistent with the proposition that an extended conformational switch favors adaptation to new selective pressures; this correlation does not depend on the choice of $\delta$ to resolve different phenotypes (Fig. S4). In contrast, $\varepsilon(a)$ does not correlate with $\phi(a)$ (Fig. S5). Importantly, while $\varepsilon(a)$ is defined from the effect of single-point mutations only, it correlates with the capacity of a system to adapt to a change of selective pressure over multiple generations (Fig. S6).

In summary, the two-dimensional spin model illustrates how selection for exquisite discrimination gives rise to a conformation switch, how this switch may involve a varying number of sites and how the potential for allostery and evolvability increases with the number of coupled sites that control the switch.

III Conservation, coevolution and thermal stability

A relevant model of protein evolution should reproduce the salient statistical patterns present in multiple sequence alignments of natural protein families. Those families comprise sequences that evolved from a common ancestral sequence under presumably similar selective pressures. The simplest way to produce comparable alignments from the model is to take an evolved sequence as the ancestral sequence, generate from it $M$ independent trajectories, and collect the sequences at the end points into an alignment. Although this procedure assumes a trivial star-like phylogeny and strictly constant selection, we find it sufficient to produce statistical features comparable to those in natural alignments.

Given an alignment, a degree of evolutionary conservation can be defined at each site $i$ by the relative entropy

[TABLE]

where $f_{i}(a_{i})$ is the frequency at which amino acid $a$ is present at position $i$ in the alignment. The pattern of evolutionary conservation $D_{i}$ essentially reproduces the patterns of conformational changes $\omega_{i,r}$ and allosteric effects $\Delta_{h_{i}}\phi$ in the ancestral sequence (Fig. 5A, to be compared with the second column of Fig. 4, which characterizes the ancestral sequence from which the alignment is generated).

Beyond conservation at individual positions, we can also analyze coevolution between sites through the correlation matrix $C_{ij}(a_{i},b_{j})=f_{ij}(a_{i},b_{j})-f_{i}(a_{i})f_{j}(b_{j})$ where $f_{ij}(a_{i},b_{j})$ is the joint frequency of $(a_{i},b_{j})$ at sites $i$ and $j$. Different statistical patterns can be extracted from $C_{ij}(a_{i},b_{j})$ Rivoire:2013cy ; Cocco:2013el . The statistical coupling analysis (SCA) thus extracts conserved global modes Rivoire:2016bl while the direct coupling analysis (DCA) extracts pairs of strongly coupled sites Morcos:2011jg . Applied to natural sequence alignments, SCA identifies groups of evolutionary conserved and coevolving sites called sectors, which are found to be structurally connected and associated with core functions of proteins Rivoire:2016bl . DCA, on the other hand, infers a ranked list of coevolving pairs, the top ones of which are found to be in contact in the three-dimensional structure Morcos:2011jg . Applying SCA to alignments generated from the model, we recover as a sector essentially the same set of sites that participate to the conformational switch and are identified from evolutionary conservation (Fig. 5), comparable to what is obtained with natural alignments of proteins with a single sector. SCA combines conservation and correlations and here the sector is mostly controlled by conservation. In alignments of natural proteins, sequences of different specificities are often present which enhance the correlations. Applying DCA yields top coevolving pairs that are in structural contact (Fig. 5B); their number, however, is small compared to what is obtained from natural alignments of similar size and diversity.

One addition to the model is sufficient to correct for this discrepancy: a selective constraint on thermal stability. Our model indeed assumes that all sequences are folded into the same cylindrical shape without considering the possibility for some mutations to destabilize this structure. A simple way to integrate this possibility without explicitly modeling folding is to impose a maximal value $F^{*}$ for the free energy $F(a,\ell_{0})$: above this value, the system is considered unfolded. To account for the fact that natural proteins are only marginally stable Taverna:2002gj , we choose $F^{*}$ well below the typical values of $F(a,\ell_{0})$ in absence of stability constraint (Fig. S7) but well above what may be obtained by minimizing $F(a,\ell_{0})$ (Fig. S8). As shown in Fig. 5B, the constraint $F(a,\ell_{0})\leq F^{*}$ alone is sufficient to generate a large number of coevolving contacts, irrespectively of the exact value of $F^{*}$ (Fig. S9). Imposing jointly the two selective pressures, exquisite discrimination and thermal stability, reproduces both features of natural alignments, a localized and evolutionarily conserved sector controlling binding affinity and specificity, and a large number of coevolving and contacting pairs of sites distributed across the structure (Fig. 5C and Fig. S10).

IV Discussion

Proteins derive their many functions from a few key properties, amongst which the capacity to stably fold into three-dimensional structures (stability), to selectively bind distinct ligands and substrates (specificity), to be regulated through long-range effects (allostery) and to adapt to changing selective pressures (evolvability). These different properties have been characterized in a number of instances but understanding how they are encoded into amino acid sequences remains a challenge. One puzzling feature is the ubiquity of long-range effects: mutational and evolutionary studies concur to indicate that binding and catalysis are affected by substitutions of amino acids more than 10 Å from the active site licata1995long . Another puzzling feature is the co-existence of two types of coevolution patterns in multiple sequence alignments: structurally connected groups of coevolving amino acids called sectors Lockless:1999uf ; Halabi:2009jc as well as a variety of more isolated pairs of coevolving amino acids in contact in the three-dimensional structure Morcos:2011jg . Beyond case-by-case descriptions and statistical modeling, what explains the mutational and evolutionary patterns that we observe in protein sequences? To what extent are the different key properties of proteins independent of each other? Without necessarily seeking to retrace natural history, what parsimonious evolutionary scenarios may generate comparable sequences and phenotypes?

Here, we proposed and analyzed such a scenario, which is based on one primary selective pressure, the requirement to achieve fine molecular recognition. First, we showed in a variety of different physical contexts, how discriminative binding requires a finely tuned conformational switch. Second, we argued that fine tuning may involve an extended subpart of the system to benefit from an enlarged evolutionary space. We then showed that extended sectors indeed emerge in a physical model with short range interactions where selection operates only locally. Third, we showed that the sensitivity of this sector to perturbations implies a potential to evolve allosteric regulation and to adapt to new selective pressures. The model explains both the conformation changes observed in protein structures and the coevolving sectors observed in protein sequences. In this scenario, adaptability is not opposed to high-performance but comes as its natural consequence. The arguments are partly independent of the physical nature of the system, which only needs to be a structured network with short-range interactions subject to evolution. We therefore chose to illustrate the scenario with one of the simplest physical models, a spin model. Finally, we showed how the second type of coevolution patterns, distributed pairs of coevolving contacts, can be explained within the same framework by taking into account an additional selective pressure, thermal stability. Our evolutionary scenario is parsimonious in the sense that only two selective pressures, discriminative binding and thermal stability, are invoked, with long-range effects arising as a by-product.

This is in sharp contrast with several physical models recently been proposed for the evolution of allostery, which are based on a direct selection for long-range regulation Hemery:2015ei ; flechsig2017design ; rocks2017designing ; tlusty2017physical ; yan2017architecture . A general evolutionary argument and specific experiments, however, argue strongly against a direct selection for long-range effects in the initial steps of the evolution of allosteric regulation. From an evolutionary perspective, the emergence of allostery by direct selection for regulation is indeed considered implausible in the biological literature as it assumes the concerted evolution of two proteins, one becoming regulated and the other one regulating it Kuriyan:2007 . As for the evolution of other complex systems, the conundrum disappears if considering the pre-existence of components and properties that evolved under unrelated selective pressures. The evolution of allosteric regulation is, instead, easily explained if long-range effects pre-exist any direct selective pressure to evolve them. The proposal that proteins have an inherent propensity for allostery has been made previously Gunasekaran:2004 and is strongly supported by a number of experimental results, starting with the repeated discovery of serendipitous allosteric sites with non-physiological effectors Hardy:2004gy ; zorn:2010 . Latent allostery is also demonstrated by the successful engineering of allosteric regulation at several surface sites of a protein Lee:2008gd ; Reynolds:2011gs ; pincus2017evolution ; these sites were not known as regulatory targets but, consistent with our model, they are evolutionarily linked through a sector to the active site. These works, and more generally the ubiquity of long-range effects in proteins licata1995long ; Lockless:1999uf ; Daugherty:2000gz ; Morley:2005kc , unequivocaly reveal the prevalence of latent allostery. Finally, compelling evidence for pre-existing long-range effects as precursors of allosteric regulation is provided by an experimental study of MAP kinases Coyle:2013dn : an effector $A$ allosterically regulating a kinase $K$ in one species of yeast was shown to regulate evolutionary related kinases $K^{\prime}$ in other species where no such regulation or analog to $A$ had evolved, in line with a scenario where that $A$ had an allosteric effect on the common ancestor of $K$ and $K^{\prime}$ prior to any selection for allosteric regulation. While convincingly indicating latent allostery as an origin of allosteric regulation, these past studies do not explain, however, how latent allostery arises in the first place. Our model is the first, to our knowledge, to provide a possible explanation.

Our evolutionary scenario also contrasts with other explanations for the origin of evolvability. In some previous models, flexibility is linked to evolvability, but evolution under constant selective pressure leads to a phenomenon of canalization where flexibility and evolvability both become increasingly limited Ancel:2000gj . This is resolved in other models by considering an evolutionary history of temporally varying selecting pressures where, for instance, the nature of the ligand changes in the course of evolution. Parter:2008bb ; Hemery:2015ei ; Raman:2016gv . While such fluctuations are likely to be relevant to protein evolution, our model shows that they are not necessary to generate evolvable proteins. The two factors, selection for exquisite discrimination and fluctuating selective pressures, are however non-exclusive and may reinforce each other. In our model, a fluctuating selection may thus contribute to select for systems with more extended conformational switches. The two factors may in fact reflect the same principle: the presence of related but partly conflicting constraints. The relationship between these constraints may be more critical than their simultaneous or successive occurrence.

The contribution of conformational changes to molecular recognition has also been explained by a different mechanism called conformational proofreading savir2007conformational . In this model, a structural mismatch between the protein and its substrates favors correct recognition by penalizing binding to the wrong ligand more than it does to the right ligand. Here, we considered a stronger notion of molecular recognition called exquisite discrimination where binding to the wrong ligand is less favorable than binding to the solvent in addition to being less favorable than binding to the right ligand. We showed exquisite discrimination to be achievable through a particular type of conformational change, which can in particular take the form of a switch between two local minima of the potential (Fig. 1A). While conformational proofreading is particularly relevant to proteins involved in search processes vlaminck , exquisite discrimination may be more relevant to enzymes, which have to bind to a reactant but release its product. Despite differences, the two mechanisms imply a similar notion of fine-tuning, which we propose here to be a sufficient factor for evolving an extended set of coupled sites.

By examining the implications of discriminative recognition, our model complements the many studies focused on the folding and thermal stability of proteins. This includes models that analyzed the interplay between binding affinity and thermal stability Miller:1997es ; Williams:2001ua ; Manhart:2015eg . Our model illustrates a simple dichotomy. On one hand, selection for exquisite discrimination leads to a conserved set of coevolving positions structurally connected to the binding site, analogous to the sectors inferred from multiple sequence alignment of natural proteins Rivoire:2016bl . On the other hand, constraints on thermal stability lead to a large number of contacting coevolving sites distributed across the structure, analogous again to what is found in natural proteins Morcos:2011jg . Consistent with this scenario, the method that best infers structural contacts from coevolution (DCA) is also very effective at scoring sequences for their thermal stability Morcos14 ; Figliuzzi15 . Our model treats unfolded conformations implicitly but the generation of coevolving contacts from selection for thermal stability was previously demonstrated in lattice protein models where unfolded conformations are explicitly considered Jacquin:2016cl . At the origin of the two distinct statistical signatures are two essential differences between binding and stability constraints. Physically, binding involves a localized interaction while stability is a structurally distributed property. Evolutionarily, selective pressures on stability are typically less stringent than selective pressures on binding, as reflected by the marginal stability of most globular proteins Taverna:2002gj and the scarcity of mutations improving binding in wild-type proteins McLaughlinJr:2012hw . As a result, the stability problem has high entropy in sequence space, with different solutions involving different isolated interactions, while the binding problem has low entropy, with solutions all involving essentially the same clustered interactions.

Partial or global unfolding induced by perturbations at sites contributing to thermal stability, which may be located arbitrarily far from the active site, is in fact another mechanism by which protein activity may respond to distant perturbations without any explicit selection for allosteric regulation. Allosteric switches based on partial or global unfolding are indeed documented in several proteins Motlagh:2014kc and can provide a mechanism for adaptation Saavedra:2018ea . For such mechanism, no evolutionarily conserved pathway linking the allosteric site to the active site is necessary.

Thermal fluctuations around a folded state also provide a way by which perturbations may propagate through long distances, possibly without major conformational change Cooper:1984cn . This mechanism is evidenced in some proteins Popovych:2006 and is explained by elastic inhomogeneities mcleish:2013 . As these inhomogeneities are generic features of proteins, their presence may also precede any selection for allosteric regulation. Following our approach, it would be interesting to study this other scenario, taking into account the demonstrated relation to coevolving contacts Townsend:2015 . This scenario may depend on the nature of the physical interactions since, at least in some proteins, fluctuation-based allostery operates through elastic motions Townsend:2015 . Spin models, where only a few discrete states are accessible at each site, may instead be more suited to a description of side-chain fluctuations between a few rotameric states attached to a fixed backbone, which dominate fluctuations in some other proteins dubay2011long . Going beyond the generic effects analyzed in this work to study evolutionary scenarios that rely on more specific physical mechanisms is a natural avenue for future work.

While other evolutionary scenarios are possible and should be considered, the present scenario can already be confronted to experimental tests. First, it is possible to analyze the physical underpinning of the different coevolution patterns. Our model predicts that mutating sector positions impacts function (binding, catalysis) and possibly also thermal stability, while mutating coevolving contacts outside of sectors impacts mostly thermal stability. The idea that multiple sequence alignments contain correlations of different nature is already a key principle behind the statistical coupling analysis (SCA), where the most conserved correlations are up-weighted to highlight functional coevolution Lockless:1999uf ; Rivoire:2016bl . This approach is supported by several experimental studies, including the design of functional proteins by reproducing the patterns of functional coevolution Socolich:2005js ; Russ:2005kc . The same principles may be extended to design sequences enhancing or ignoring other types of correlations found in multiple sequence alignments. For instance, in the context of the direct coupling analysis (DCA), which also provides a generative model figliuzzi2018pairwise , we may predict that designing sequences at low statistical temperature will lead to more stable proteins.

A direct experimental test of the central idea that a selection for exquisite discrimination is sufficient to generate an extended sector is also possible. This can be done in two steps. First, proteins can be evolved to have desired specificities, using methods of directed evolution such as phage or yeast display Smith:1997vn ; Levin:2006 . These methods allow for the selection of proteins with high affinity to an arbitrary target ligand from populations of billions of different variants of a protein. This selection can be repeated with intervening steps of mutations and amplification, thus emulating in vitro evolution by natural selection. Specificity can be controlled by combining positive selection for binding to immobilized targets with negative selection for not binding to soluble targets that are washed away. The approach is powerful enough to generate proteins specific to a particular conformational state of a target protein Nizak:2003 or, on the contrary, cross-reactive to a range of variants of a target protein Garcia:2010 . Second, the presence of long-range effects can be assayed by screening and sequencing mutants of the evolved proteins Fowler:2014gq . Starting from an artificially evolved protein, a population of mutants can be produced and selected for specific binding: using high-throughput sequencing, the frequency of each mutant before and after selection can be known, which reveals which mutations have an effect on binding. Our prediction is that proteins evolved to be highly specific will exhibit more deleterious mutations far from their binding site. As specificity can be finely controlled during selection and as sensitivity to mutations can be assayed quantitatively and mapped to three-dimensional structures within a few angstroms, these experiments can not only test the prediction of our toy model but provide data for elaborating more accurate models.

Acknowledgements.

This work benefited from multiple inspiring discussions with Rama Ranganathan, Clément Nizak and Anton Zadorin, and from comments by Yaakov Kleeorin. It was supported by FRM AJE20160635870.

A. Fitness function for exquisite discrimination

Two ligands $\ell_{r}$ and $\ell_{w}$ with binding free energies $\Delta F_{r}(a)$ and $\Delta F_{w}(a)$ and chemical potentials $\mu_{r}$ and $\mu_{w}$ have, at thermal equilibrium, probabilities $p_{r}(a)$ and $p_{w}(a)$ to bind, with

[TABLE]

where $k=r$ or $w$. If binding to $\ell_{r}$ is desirable and binding to $\ell_{w}$ is not, the sequence $a$ has to maximize $p_{r}(a)$ and minimize $p_{w}(a)$. The problem may be formalized as the optimization of a fitness function $\phi(a)$ that is a decreasing function of $\Delta F_{r}(a)$ and an increasing function of $\Delta F_{w}(a)$.

The choice of $\phi(a)=\min(-\Delta F_{r},(a),\Delta F_{w}(a))$ leads to one of the simplest fitness functions that captures this trade-off. Imposing $\phi(a)>0$ is equivalent to imposing $F(a,\ell_{r})<F(a,\ell_{0})<F(a,\ell_{w})$. As another important factor controlling the difficulty of discrimination is the relative concentration of the two ligands. A natural generalization is therefore $\phi(a)=\min(-\Delta F_{r}(a)+\mu,\Delta F_{w}(a)-\mu)$, so that $\phi(a)>0$ is equivalent to imposing $F(a,\ell_{r})-\mu<F(a,\ell_{0})<F(a,\ell_{w})-\mu$. We show in Fig. S2 that taking $\mu=-1$ or $\mu=+1$ gives similar results than $\mu=0$.

To further demonstrate the robustness of the results to the choice of the fitness function, we also show that the less physical choice of $\phi(a)=-\Delta F_{r}(a)\times\Delta F_{w}(a)$, corresponding to the requirement that $\Delta F_{r},(a)$ and $\Delta F_{w}(a)$ should both have large absolute values but opposite signs, also gives similar results (Fig. S2).

B. Unidimensional models

The potentials $U(x,a,\ell)=U(x,h=\ell-a)$ of the models presented in Fig. 1 have a common symmetry:

[TABLE]

The free energy $F(h)=-\beta^{-1}\ln\int dx\ e^{-\beta U(x,h)}$ is therefore symmetric around $h=0$: $F(h)=F(-h)$, where $F(h)$ is a decreasing function of $h$ for $h>0$. Given $\ell_{0},\ell_{r},\ell_{w}$, three different cases thus arise when considering the existence of $a$ satisfying the condition $F(a,\ell_{r})<F(a,\ell_{0})<F(a,\ell_{w})$:

(i) if $\ell_{w}<\ell_{0}<\ell_{r}$ (respectively, $\ell_{r}<\ell_{0}<\ell_{w}$), the condition is satisfied by taking any $a<\ell_{w}$ (resp. any $a>\ell_{w}$) so that $\ell_{w}-a,\ell_{0}-a,\ell_{r}-a$ are of same sign.

(ii) if $\ell_{0}<\ell_{w}<\ell_{r}$ (respectively, $\ell_{r}<\ell_{w}<\ell_{0}$), the condition is satisfied by taking $\ell_{w}-a<-\ell_{0}+a<\ell_{r}-a$, i.e. $(\ell_{0}+\ell_{w})/2<a<(\ell_{0}+\ell_{r})/2$ (resp. $\ell_{w}-a<-\ell_{0}+a<\ell_{r}-a$, i.e. $(\ell_{0}+\ell_{r})/2<a<(\ell_{0}+\ell_{w})/2$) so that $\ell_{w}-a,\ell_{r}-a$ on one hand and $\ell_{0}-a$ on the other are of opposite signs.

(iii) if $\ell_{0}<\ell_{r}<\ell_{w}$ (respectively, $\ell_{w}<\ell_{r}<\ell_{0}$), the condition cannot be satisfied.

As $\langle x\rangle_{-h}=-\langle x\rangle_{h}$, (ii) involves a change of sign for the mean conformation while (i) does not. For the models of Fig. 1, solutions to (i) in fact exist that have arbitrarily small relative conformational changes, $(\langle x\rangle_{h_{r}}-\langle x\rangle_{h_{0}})/\langle x\rangle_{h_{0}}\to 0$, obtained for $a\to\pm\infty$.

The problem of exquisite discrimination is to be compared to the problem of maximizing the affinity to a single target $\ell_{r}$. As $\partial^{2}F(h)/\partial h^{2}=-(\langle x^{2}\rangle_{h}-\langle x\rangle_{h}^{2})\leq 0$, $F(h)$ is concave. In general, it is even strictly concave, the single-spin model at zero temperature where $F(h)=-|h|$ being a degenerate case. This implies that $F(a,\ell_{r})-F(a,\ell_{0})$ is minimized by taking $a\to+\infty$ when $\ell_{r}>\ell_{0}$ and $a\to-\infty$ when $\ell_{r}<\ell_{0}$. More generally, the condition $F(a,\ell_{r})<F(a,\ell_{0})$ may be seen as a limit of the condition $F(a,\ell_{r})<F(a,\ell_{0})<F(a,\ell_{w})$ when ${\rm sign}(\ell_{0})(\ell_{0}-\ell_{w})\to 0^{+}$, corresponding to case (i).

Single-spring elastic network

The model of Fig. 1A is defined by

[TABLE]

which verifies $U(x,-h)=U(-x,h)$. Here, $k>0$ represents the stiffness of the spring and $r>0$ its equilibrium length. While it is possible to write a general analytical formula for the free energy $F(h)$ at any $\beta$, it is more illuminating to consider the zero-temperature limit $\beta\to\infty$. In this limit, $F(h)=\min_{x}U(x,h)$ is reached for $x=\hat{x}(h)$ with

[TABLE]

$F(h)$ and $\hat{x}(h)$ are represented in Fig. 1A for $r=1$, $k=2$, $h_{0}=-1/2$, $h_{w}=1/4$, $h_{r}=3/4$.

Two-spring elastic network

The model of Fig. 1B is defined by

[TABLE]

which verifies $U(x,-h)=U(-x,h)$. Here, $k>0$ represents the common stiffness of the two springs, $r>0$ their common equilibrium length and $2d>0$ the distance between the two fixed points at which they are attached. In the zero temperature limit $\beta\to\infty$, the free energy is the minimum of $U(x,h)$, reached for $\hat{x}(h)$ solution to

[TABLE]

In Fig. 1B, the graph of $F(h)$ versus $h$ is obtained by parametrizing $F(h)$ and $h$ by $\hat{x}$: $F(h)=U(\hat{x},h(\hat{x}))$ and $h(\hat{x})=2k\hat{x}(1-r(\hat{x}^{2}+d^{2})^{-1/2})$. Fig. 1B corresponds to $k=1$, $r=1$, $h_{0}=-1/2$, $h_{w}=1/4$, $h_{r}=3/4$. Depending on the value of $r$, $U(x,h)$ may have either one or two local minima (Fig. S1).

Single-node harmonic model

A single-node Gaussian elastic network model is defined by

[TABLE]

It verifies $U(x,-h)=U(-x,h)$ only for $r=0$ but the model is physically equivalent to a model with potential $V(x,h)=U(x+r,h)+hr$ which verifies $V(x,-h)=V(-x,h)$ for any $r$. The conclusions are therefore the same as for the previous models. The free energy is

[TABLE]

It corresponds graphically to an inverted parabola, qualitatively similar to $F(h)$ for the models of Fig. 1.

Shape space model

In shape space models Perelson:1979gb , no flexibility is considered and the conformations $x$ are formally identified with the sequences $a$. The sequences $a$ and ligands $\ell$ are points in a $d$-dimensional Euclidean space and the binding energy between $a$ and $\ell$ is a function of their Euclidean distance, in the simplest case $F(a,\ell)=\|a-\ell\|_{2}$. In dimension $d=1$, this corresponds to $F(h)=|h|$ and exquisite discrimination is possible when $\ell_{0}<\ell_{r}<\ell_{w}$ provided $(\ell_{0}+\ell_{r})/2<a<(\ell_{0}+\ell_{w})/2$, similar to the criterion derived for elastic models.

Single-spin model

The model of Fig. 1C is defined by

[TABLE]

which verifies $U(x,-h)=U(-x,h)$. The free energy is

[TABLE]

In the zero-temperature limit $\beta\to\infty$, $F(h)=-|h|$, which is identical up to the sign to the free energy of a one-dimensional shape space model.

Let assume $\ell_{0}=0$ and $\ell_{r}>0$ to analyze the nature of the solutions satisfying $F_{r}<F_{0}<F_{w}$ where $F_{r}=F(a,\ell_{r})$, $F_{0}=F(a,\ell_{0})$ and $F_{w}=F(a,\ell_{w})$. We proceed by examining the different possible cases:

1. if $a>0$, then $F_{0}=-a$ and $F_{r}=-a-\ell_{r}$ satisfies $F_{r}<F_{0}$.

1.1. if $a<-\ell_{w}$, then $F_{w}=a+\ell_{w}$ and $F_{0}<F_{w}$ implies $-a<a+\ell_{w}$, i.e., $-\ell_{w}/2<a$. In total, $-\ell_{w}/2<a<-\ell_{w}$. This is possible provided $\ell_{w}<0$. In this case, $F_{r}-F_{0}=-\ell_{r}<0$ and $F_{w}-F_{0}=\ell_{w}+2a\leq-\ell_{w}$, whose maximum is reached for $\hat{a}=-\ell_{w}$, yielding $F_{w}-F_{0}=-\ell_{w}>0$. A solution exists in this case.

1.2. if $-\ell_{w}<a$, then $F_{w}=-a-\ell_{w}$ and $F_{0}<F_{w}$ implies $-a<-a-\ell_{w}$, i.e., $\ell_{w}<0$. In total, $0<-\ell_{w}<a$. This is possible provided $\ell_{w}<0$. In this case, $F_{r}-F_{0}=-\ell_{r}<0$, $F_{w}-F_{0}=-\ell_{w}>0$ independent of the value of $a$. A solution exists in this case.

2. if $a<0$, then $F_{0}=a$.

2.1. if $a<-\ell_{r}$, $F_{r}=a+\ell_{r}$ and $F_{r}<F_{0}$ implies $a+\ell_{r}<a$, which is inconsistent with the assumption that $\ell_{r}>0$. No solution exists in this case.

2.2. if $-\ell_{r}<a$, $F_{r}=-a-\ell_{r}$ and $F_{r}<F_{0}$ implies $-a-\ell_{r}<a$, i.e., $-\ell_{r}/2<a$.

2.2.1. if $a<-\ell_{w}$, $F_{w}=a+\ell_{w}$ and $F_{0}<F_{w}$ implies $a<a+\ell_{w}$, i.e., $\ell_{w}>0$. In total, $-\ell_{r}/2<a<-\ell_{w}<0$. This is possible provided $0<\ell_{w}<\ell_{r}/2$. In this case, $F_{0}-F_{r}=2a+\ell_{r}\leq\ell_{r}-2\ell_{w}$, $F_{w}-F_{0}=\ell_{w}$ with maximum reached for $\hat{a}=-\ell_{w}$.

2.2.2. if $-\ell_{w}<a$, $F_{w}=-a-\ell_{w}$ and $F_{0}<F_{w}$ implies $a<-a-\ell_{w}$, i.e., $a<-\ell_{w}/2$. In total, $\max(-\ell_{r}/2,-\ell_{w})<a<-\ell_{w}/2$. Possible provided $0<\ell_{w}<\ell_{r}$. $F_{0}-F_{r}=2a+\ell_{r}\leq\ell_{r}-\ell_{w}$, $F_{w}-F_{0}=-2a-\ell_{w}\leq\min(\ell_{r}-\ell_{w},\ell_{w})$, which cannot be both maximized. If maximizing $\min(F_{0}-F_{r},F_{w}-F_{0})$, the optimum is for $\hat{a}=-(\ell_{r}+\ell_{w})/2$ in which case $F_{0}-F_{r}=F_{w}-F_{0}=(\ell_{r}-\ell_{w})/2$.

All together, the nature of the solution depends on the ratio $\ell_{w}/\ell_{r}$:

– if $\ell_{w}/\ell_{r}<0$: two solutions, $\hat{x}_{0}=+$, $\hat{x}_{r}=+$, $\hat{x}_{w}=\pm$, with $\hat{x}_{w}=+$ for the solution optimizing $\phi(a)$;

– if $0<\ell_{w}/\ell_{r}<1/2$: two solutions $\hat{x}_{0}=-$, $\hat{x}_{r}=+$, $\hat{x}_{w}=\pm$, with $\hat{x}_{w}=-$ for the solution optimizing $\phi(a)$;

– if $1/2<\ell_{w}/\ell_{r}<1$: one solution $\hat{x}_{0}=-$, $\hat{x}_{r}=+$, $\hat{x}_{w}=+$;

– $1<\ell_{w}/\ell_{r}$: no solution.

A switch is thus involved whenever $\ell_{0}<\ell_{w}<\ell_{r}$.

C. Two-dimensional spin model

Given fields $h_{i}(a_{i})$ and couplings $J_{ij}(a_{i},a_{j})$ defined respectively on each node $i$ and each edge $ij$ of the lattice represented in Fig. 2A, and given an external field $\ell$ applied to the red node $b$ in Fig. 2A, the free energy $F(a,\ell)$ is computed exactly from Eq. (1) by transfer matrices using free boundary conditions for the spins at the left and right extremities. Performing the calculation from right to left allows for an efficient evaluation of $F(a,\ell)$ for different values of $\ell$ since in this case changing $\ell$ affects only the last layer.

The capacity of a system to discriminate between two ligands $\ell_{r}$ and $\ell_{w}$ is scored by $\phi(a)$ given by Eq. (2). To evolve systems under a selection for discrimination, a Metropolis Monte Carlo procedure is followed. It starts from a random sequence and generates a trajectory by iterating $T$ times the following steps: a site $i$ is chosen at random, an $a_{i}^{\prime}\neq a_{i}$ is chosen at random between the $q-1$ possibilities, and the substitution is accepted if $r<\exp(\gamma(\phi(a^{\prime})-\phi(a)))$ where $r$ is a random number uniformly drawn in $[0,1]$ and $a^{\prime}$ is sequence $a$ with $a_{i}^{\prime}$ substituted for $a_{i}$. The parameters $\gamma=100$ and $T=1000$ are taken. As shown in Fig. S12A, this number of iterations is sufficient for $\phi(a)$ to reach an equilibrium value.

This value is in most case close to the theoretical maximum that $\phi(a)$ may reach (Fig. 3). This maximum is the same for the two-dimensional spin model and for the single-spin model. It is reached when $F(h_{0}=a+\ell_{0})=(F(h_{r}=a+\ell_{r})+F(h_{w}=a+\ell_{w}))/2$ with $F(h)$ given by Eq. (S10). With the parameters taken in the main text, $\max_{a}\phi(a)\simeq 0.41$.

The different systems in Fig. 4 correspond to different choices of the function $h_{i}(a_{i})$ and $J_{ij}(a_{i},a_{j})$. Different solutions may also be obtained with the same mapping between sequences $a$ and fields $h_{i}(a_{i})$ and couplings $J_{ij}(a_{i},a_{j})$ but different initial conditions or/and different series of proposed mutations. The system in the second column of Fig. 4 (with $\bar{\omega}_{r}=0.11$) is evolved under an additional constraint for thermal stability to serve as a common ancestral sequence for the alignments of Fig. 5.

Thermal stability

A system is considered thermally stable if its free energy is below a fixed threshold, $F(a,\ell_{0})<F^{*}$. In Fig. 5, the value $F^{*}=-350$ is chosen for being smaller than the free energy of systems evolved without stability constraints (Fig. S7, top left panel) but larger than the free energy of systems optimized for thermal stability (Fig. S8).

To evolve stable and functional sequences, we first generate a stable sequence by a Metropolis Monte Carlo procedure where the scoring function is $-F_{0}(a)$ rather than $\phi(a)$ and stop the iterations as soon as the condition $F(a,\ell_{0})<F^{*}$ is fulfilled. Starting from this stable sequence, we then generate a trajectory as before with the only modification that mutations causing $F(a,\ell_{0})>F^{*}$ are systematically rejected (Fig. S11). The systems thus produced are marginally stable with free energies $F(a,\ell_{0})$ close to the threshold $F^{*}$ (Fig. S7).

Characterization

The mean conformations $\langle x_{i}\rangle_{a,\ell_{r}}$ entering the definition of $\omega_{i,r}(a)$ in Eq. (3) are computed exactly by transfer matrices. In Fig. 3, only the fittest evolved systems with $\phi(a)>0.35$ are considered, which represent 94% of 1000 evolved sequences. How $\bar{\omega}_{r}(a)$ relates to $\phi(a)$ is shown in Fig. S5B.

To quantify evolvability in Fig. 3B, we fix an interval of phenotypes of interest $[\ell_{\rm min},\ell_{\rm max}]=[-1,4]$, partition it into small subintervals of length $\delta=0.1$, and for each single mutant $a^{\prime}$ of a given sequence $a$ find, if it exists, the interval to which $\ell_{r}(a^{\prime})$ belongs. The fraction of subintervals covered by the $(q-1)W(L+1)$ single-point mutants of $a$ defines $\varepsilon(a)$, a measure of the phenotypic diversity its neighborhood in sequence space. How $\varepsilon(a)$ relates to $\phi(a)$ is shown in Fig. S5D. How the results depend on the choice of $\delta$ is shown in Fig. S4.

In Fig. 4, the sizes of the red dots in A are proportional to $\omega_{i,r}$ defined in Eq. (3), the sizes of the green dots in B are proportional to $\Delta_{h_{i}}\phi$, defined for each site $i$ as the maximal decrease in fitness $\phi$ when adding either $h_{i}^{\prime}=-2$ and $h_{i}^{\prime}=+2$ to $h_{i}(a_{i})$ at $i$, and the sizes of the blue dots in C are proportional to $\Delta_{a_{i}}\phi$, defined for each site $i$ as the mean loss in fitness when considering all possible $q-1$ mutations of $a_{i}$. The distribution of the effects of all single and double mutations are shown in Fig. S12B-C.

Multiple sequence alignments

The three alignments of Fig. 5 are all generated from the same ancestral sequence, a stable and functional sequence that corresponds to the second column of Fig. 4 (Fig. S11). The $M=1000$ sequences in each alignment are obtained from independent trajectories starting from this sequence.

For the first alignment (labelled “discrimination”), no constraints of stability is enforced (formally $F^{*}=\infty$) and selection is based on $\phi(a)$ given by Eq. (2). For the second alignment (“stability”), no constraint on discrimination is enforced (formally $\gamma=0$) and selection is limited to $F(a,\ell_{0})<F^{*}$. For the third alignment (“discrimination + stability”), the two selective constraints are jointly imposed. The distributions of $F(a,\ell_{0})$ and $\phi(a)$ in these three alignments are shown in Fig. S7. Taking a larger or smaller value of $F^{*}$ leads to similar conclusions (Fig. S9).

Inference of coevolution

Given an alignment, coevolving contacts are inferred by the pseudo-likelihood maximization method for direct coupling analysis Ekeberg:2013fq_re using an implementation in Julia juliaDCA with default parameters except for a reduced alphabet of $q=5$ amino acids. The algorithm outputs a ranked list of coevolving pairs. Contacts are defined based on the distance $d$ between nodes in the lattice of Fig. 2A: $d=1$ refers to two nodes connected by an edge and $d\geq 2$ to two nodes connected by $d$ edges through $d-1$ intermediate nodes.

For Fig. 4, the SCA matrix $\tilde{C}_{ij}$ is calculated as described in Box 1 of Ref. Rivoire:2016blre without any preprocessing except for a regularization parameter $\lambda=0.1$. Background frequencies are taken to be $1/q$ for all $q=5$ amino acids. Fig. S4 shows the mean values $\langle\tilde{C}_{ij}\rangle_{j}$ of each row $i$ of the matrix $\tilde{C}_{ij}$.