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

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.

Figure 1

Elementary models illustrating how exquisite discrimination implies evolutionary fine-tuning and a conformational change. (a) Model consisting of a single particle attached to a fixed point by a spring and subject to a constant force $h$. The potential is $U(x,h)=k{\left(\right|x|-r)}^2/2-hx$ and the free energy $F\left(h\right)$ is shown as a function of $h$ for $r=1,k=2,\beta =\infty$. The state $\widehat{x}\left(h\right)$ of the system here is simply the minimum of $U(x,h)$. The constant force $h$ is jointly controlled by the evolutionary and environmental degrees of freedom $a$ and $\ell$ through $h=\ell -a$. To achieve exquisite discrimination in the form $F(a,{\ell }_r)<F(a,{\ell }_0)<F(a,{\ell }_w)$ the variable $a$ must be chosen so that $F(h_0={\ell }_0-a)$ lies in between $F(h_w={\ell }_w-a)$ and $F(h_r={\ell }_r-a)$. When ${\ell }_0<{\ell }_w<{\ell }_r$, this requires $0\le h_w<-h_0<h_r$. In this case, binding induces a conformational switch from $\widehat{x}\left(h_0\right)<0$ to $\widehat{x}\left(h_r\right)>0$. (b) Model with two springs under similar evolutionary constraints. Here, the conformational change is not a switch between two states of $U(x,h_0)$, although it is for other values of the parameters (Appendix pp2-s2). (c) Spin model where $U(x,h)=-hx$, exhibiting the same phenomenology in a minimal setting where $x$ takes only two values $\pm 1$.

Figure 2

Two-dimensional spin model. (a) Lattice on which the model is defined. Each node $i$ carries a physical variable $x_i\in \{-1,+1\}$ and an evolutionary variable $a_i\in \{1,\dots ,q\}$. One node at a boundary, here in red, defines the binding site, which may be subject to an external field $\ell$ representing a ligand. The thinner lines from top to bottom nodes implement the periodic boundary conditions, giving to the structure the geometry of a cylinder. Here and in the following figures, $L=10$ and $W=5$. (b) The evolutionary variables $a_i$ and $a_j$ define fields $h_i\left(a_i\right)\text{and}{h}_j\left(a_j\right)$ and couplings $J_{ij}(a_i,a_j)$ to which all physical variables $x_i$ and $x_j$ connected by a link in the lattice are subject. The potential $U(x,a,\ell )$ is given by Eq. (1).

Figure 3

Extent of the conformational change and evolvability for different realizations of the two-dimensional spin model. (a) Distribution of ${\overline{\omega }}_r$, the effective fraction of sites changing conformation upon binding, over 1000 different evolved systems $a$. Only the fittest systems with $\varphi >0.35$ are considered. As the distribution of $\varphi$ is picked around its theoretical maximum $\varphi \simeq 0.41$, these top systems represent $>94\%$ of the total. For these systems, ${\overline{\omega }}_r$ is nearly independent of $\varphi$ (Appendix pp3-s2). Inset: Cumulative distribution of $\varphi$ over all 1000 systems. The red dotted line indicates the theoretical maximum. (b) For the same systems, the measure of evolvability $\varepsilon$ is positively correlated with ${\overline{\omega }}_r$, although independent of $\varphi$ (Appendix pp3-s2): the more extended the conformational switch, the wider the range of phenotypes available in the neighborhood of a sequence, irrespective of its fitness.

Figure 4

Five evolved systems with equivalent fitness $\varphi \simeq 0.4$ but switches of different extensions ${\overline{\omega }}_r$, as indicated on the top. (a) Conformational change ${\omega }_{i,r}$ at each position $i$ when ${\ell }_r$ is substituted for ${\ell }_0$. (b) Maximal fitness cost ${\mathrm{\Delta }}_{h_i}\varphi$ when adding an external field $h_i^{\prime }=\pm 2$. (c) Average fitness cost ${\mathrm{\Delta }}_{a_i}\varphi$ of mutating $a_i$. In each case, the size of the dots is proportional to the quantity of interest. The first two rows are nearly identical, indicating that the same sites that change conformation are allosteric. The relation to sites displaying a strong mutational effect (last row) is also apparent although less straightforward due to the impact of mutations on both $h_i\left(a_i\right)$ and $J_{ij}(a_i,a_j)$. These examples illustrate how long range effects may arise from local selection for exquisite discrimination [the binding site is always the middle node on the left edge as in Fig. 2].

Figure 5

Conservation and coevolution in multiple sequence alignments. Three alignments of 1000 sequences were generated from the same ancestral sequence, corresponding to the second column of Fig. 4, but under different selective pressures: (a) exquisite discrimination alone, (b) stability alone, and (c) the two selections jointly. The first row indicates how the measure of evolutionary conservation $D_i$ defined in Eq. (4) varies from site to site. The second row indicates the mean value of the SCA correlation matrix ${\langle {\tilde{C}}_{ij}\rangle }_j$ for each site $i$, a quantity that combines conservation and coevolution [21]. The third row displays the performance of contact prediction by plmDCA [23]: pairs of sites are ranked from most to least coevolving and for the top $r$ pairs the fraction in contact in the cylindrical structure is reported as a function of $r$, with contact defined either as distance $d=1$ (in blue), $d=2$ (orange), or $d=3$ (green). Selection for discrimination produces a set of evolutionary conserved and correlated sites that correspond to sites showing a strong mutational effect in the ancestral sequence (second column, third row of Fig. 4), but only few coevolving contacts. Selection for stability produces a different pattern, without strongly conserved sites but with many coevolving contacts. Finally, joint selection for discrimination and stability generates both a set of conserved sites and many coevolving contacts.

Figure 6

Graphs of Fig. 3 for alternative choices of the fitness function $\varphi \left(a\right)$. (a) $\varphi \left(a\right)=min(-\mathrm{\Delta }{F}_r\left(a\right)+\mu ,\mathrm{\Delta }{F}_w\left(a\right)-\mu )$ with $\mu =1$. (b) $\varphi \left(a\right)=min(-\mathrm{\Delta }{F}_r\left(a\right)+\mu ,\mathrm{\Delta }{F}_w\left(a\right)-\mu )$ with $\mu =-1$. (c) $\varphi \left(a\right)=-\mathrm{\Delta }{F}_r,\left(a\right)\times \mathrm{\Delta }{F}_w\left(a\right)$. In any case, the results are qualitatively similar to those of Fig. 3, which corresponds to $\varphi \left(a\right)=min(-\mathrm{\Delta }{F}_r\left(a\right)+\mu ,\mathrm{\Delta }{F}_w\left(a\right)-\mu )$ with $\mu =0$.

Figure 7

Potential $U(x,h)$ for the model of Fig. 1 with springs of different equilibrium length $r$. $U(x,h)$ is given by (B4) with here $k=1,d=1,h_0=-1/2,h_w=1/4,h_r=3/4$ such that ${min}_{x}U(x,h_r)<{min}_{x}U(x,h_0)<{min}_{x}U(x,h_w)$, i.e., the global minimum of the blue curve is in between the global minima of the red and green curves. For $r\le 1,U(x,h)$ has a single minimum while for $r>1$ it has two minima. Figure 1 corresponds to $r=1$.

Figure 8

Further characterization of the five systems presented in Fig. 4. (a) Evolution of the fitness $\varphi \left(a\right)$ as a function of the number of iterations in the Metropolis Monte Carlo algorithm. The systems presented in Fig. 4 are the end points of these trajectories ($T=1000$). The initial conditions are random sequences, except for the second system (${\overline{\omega }}_r=0.11$), which, for the purpose of Fig. 5, is started from a stable sequence satisfying $F(a,{\ell }_0)<F^{*}=-350$ (Fig. 11). All systems converge to $\varphi \simeq 0.4$, near the theoretical maximum ${max}_a\varphi \left(a\right)=0.41$, although with different extensions ${\overline{\omega }}_r$ of their sector. (b) Distributions of mutational effects ${\mathrm{\Delta }}_{a_i}\varphi$ when considering all $(q-1)W(L+1)=275$ possible single mutations. In Fig. 4, the sizes of the blue dots in the last row indicate for each $i$ the mean value of ${\mathrm{\Delta }}_{a_i}\varphi$ over the $q-1$ possible mutations. The distributions share a peak around zero (neutral mutations) and a more or less extended negative tail, commensurate with ${\overline{\omega }}_r$. The peak around ${\mathrm{\Delta }}_{a_i}\varphi \simeq -2.4$ that is observed in some cases corresponds to mutants that are effectively equivalent to a single-spin model with $a^{\prime }=0$, in which case $\varphi \left(a^{\prime }\right)=min\left(\right|{\ell }_r|,-|{\ell }_w\left|\right)=-{\ell }_w=-2$.

Figure 9

Distribution of the free energy $F(a,{\ell }_0)$ and fitness $\varphi \left(a\right)$ over the sequences of the three alignments studied in Fig. 5. The first column corresponds to an alignment selected for exquisite discrimination: the sequences have relatively high free energy $F(a,{\ell }_0)>F^{*}=-350$ and nearly optimal fitness $\varphi \left(a\right)\simeq 0.41$. The second column corresponds to an alignment generated under the constraint for stability $F(a,{\ell }_0)<F^{*}=-350$: the sequences have a free energy $F(a,{\ell }_0)$ close to the stability threshold $F^{*}$ but low fitness $\varphi \left(a\right)<0$. The third column corresponds to an alignment selected for exquisite discrimination under a constraint for stability: the sequences have both a free energy below $F^{*}$ and a near optimal fitness.

Figure 10

Free energy of systems optimized for stability, i.e., evolved with $-F(a,{\ell }_0)$ as fitness function rather than $\varphi \left(a\right)$. The resulting free energies reach $F_{\min }\simeq -450$, far below the free energies of the system selected for discrimination without stability constraint [$F(a,{\ell }_0)>F_0=-325$ in Fig. 9]. To obtain marginally stable systems, the stability threshold $F^{*}$ is chosen to satisfy $F_{\min }<F^{*}<F_0$. In Fig. 5, $F^{*}=-350$ with other choices being considered in Fig. 16.

Figure 11

Evolution of a stable and discriminating system. First, the system is evolved by Metropolis Monte Carlo with selection for thermal stability, using $-F(a,{\ell }_0)$ as fitness function. The iterations are stopped as soon as the required minimal stability, here $F(a,{\ell }_0)<F^{*}=-350$, is reached (red dotted lines). $T=1000$ additional iterations are then performed with selection for discrimination and constraint on stability, using $\varphi \left(a\right)$ as fitness function and excluding any mutation leading to $F(a,{\ell }_0)>F^{*}$. (a) Evolution of $F(a,{\ell }_0)$. (b) Evolution of $\varphi \left(a\right)$. The end point of this evolution is the system characterized in the second column of Fig. 4, which serves as the common ancestral sequence for the three alignments of Fig. 5.

Figure 12

Dependence on the ligand ${\ell }_w$. (a) Binding free energies $\mathrm{\Delta }{F}_r$ and $\mathrm{\Delta }{F}_w$ as a function of ${\ell }_w$ for ${\ell }_0=0$ and ${\ell }_r=3$ (Figs. 3 and 4 correspond to ${\ell }_w=2$). For ${\ell }_w<{\ell }_0$, evolved systems verify $\varphi =min(-\mathrm{\Delta }{F}_r,\mathrm{\Delta }{F}_w)=\mathrm{\Delta }{F}_w<-\mathrm{\Delta }{F}_r$, while, for ${\ell }_w>{\ell }_0,\varphi =min(-\mathrm{\Delta }{F}_r,\mathrm{\Delta }{F}_w)=-\mathrm{\Delta }{F}_r=\mathrm{\Delta }{F}_w$ and the two constraints are equally at play. When ${\ell }_w>{\ell }_r$, achieving $\varphi >0$ is not possible. (b) Extensions of the conformational switches ${\overline{\omega }}_r$ and ${\overline{\omega }}_w$ as function of ${\ell }_w$ for the same systems. For ${\ell }_w<{\ell }_0$, no switch occurs (${\overline{\omega }}_r={\overline{\omega }}_w=0$), while for ${\ell }_w>{\ell }_0$ switches are obtained that display a wide range of extensions. For each value of ${\ell }_w$, the dots represent a mean and the error bars represent a standard deviation over 100 different evolved systems [Fig. 3 shows the full distribution of ${\overline{\omega }}_r$ for ${\ell }_w=2$].

Figure 13

Relationships between switch extension, evolvability, and fitness. (a, b) Switch extension ${\overline{\omega }}_r\left(a\right)$ as a function of the fitness $\varphi \left(a\right)$ over 1000 different evolved systems, showing either the systems with $\varphi \left(a\right)>0$ [98% of the total, panel (a)] or those with $\varphi \left(a\right)>0.35$ [94% of the total, panel (b)]. For these fittest systems, on which Fig. 3 is based, ${\overline{\omega }}_r\left(a\right)$ is nearly independent of $\varphi \left(a\right)$. (c, d) Evolvability $\varepsilon \left(a\right)$ as a function of the fitness $\varphi \left(a\right)$ for the same systems. For the fittest systems [$\varphi \left(a\right)>0.35$, panel (d)], $\varepsilon \left(a\right)$ is also nearly independent of $\varphi \left(a\right)$.

Figure 14

Graphs of Fig. 3 for alternative choices of the interval $\delta$ used to resolve the different phenotypes. The scale ( $y$ axis) differs in each case but the same trend is observed irrespective of $\delta$: systems with larger ${\overline{\omega }}_r$ have significantly larger $\epsilon$.

Figure 15

Evolutionary response to a change of selective pressures of systems with different sector sizes and evolvability $\varepsilon$. (a) Evolutionary trajectories consecutive to a change of selective pressure at $t=0$ from $({\ell }_r,{\ell }_w)=(3,2)$ to $({\ell }_r,{\ell }_w)=(2.5,1.5)$, starting from the five systems of Fig. 4, here labeled 1 to 5. The graph shows the average fitness $\varphi$ and its standard deviation over 100 independent trajectories as a function of the number $t$ of generations. It illustrates how the systems with a large sector (3,4,5) adapt more efficiently than the systems with a small sector (1,2). (b) Same as in panel (a) but for a change of selective pressure from $({\ell }_r,{\ell }_w)=(3,2)$ to (3.5,2.5). (c) More systematic analysis starting from the $>940$ systems with $\varphi >0.35$ considered in Figs. 13 and 13 and performing at $t=0$ the same change of selective pressure as in panel (a). The fitnesses of the systems after $t=50$ and 100 generations are compared to their evolvability $\varepsilon$. This shows that $\varepsilon$, despite being defined based on the effect of single mutations only, is informative of the capacity of a system to adapt over multiple generations. (d) Same as in panel (c) for the change of selective pressure used in panel (b).

Figure 16

Varying the stability threshold $F^{*}$ similarly to Figs. 5 and 9 but considering only alignments produced under both selection for exquisite discrimination and constraint for stability $F(a,{\ell }_0)<F^{*}$ (‘discrimination+stability”). The three alignments are generated from a common sequence with $F(a,{\ell }_0)<-400$ that is different from the ancestral sequence of alignments in Fig. 5, and then evolved using different values of $F^{*}$, as indicated on the top. In each case, a sector appears as evolutionary conserved. Coevolving contacts are inferred in every case, although in larger number when the constraint for stability is higher ($F^{*}$ smaller).

Figure 17

Evolutionary conservation of the top 50 pairs returned by DCA for the three alignments of Fig. 5. Blue, distribution of the conservation $D_i$ over all sites; orange, distribution of $D_i$ over the 100 sites involved in the top 50 DCA pairs (sites involved in multiple pairs are counted multiple times). Except for the extremely conserved sites of the alignment selected for discrimination and stability (last graph), the two distributions are very similar. This indicates that the coevolving pairs predicted by DCA are distributed all over the structure, independently of the sector inferred by SCA that corresponds to the most conserved sites (Fig. 4).