Mechanics of allostery: contrasting the induced fit and population shift scenarios

TL;DR

This study investigates the mechanical principles underlying allosteric regulation in proteins, contrasting induced fit and population shift models, and identifies optimal elastic properties that support cooperative binding and protein function.

Contribution

The paper provides empirical evidence and a mechanical model analysis showing how elastic modes influence allosteric mechanisms and their robustness to mutations.

Findings

Allosteric conformational change occurs along a soft elastic mode with high shear regions.

Optimal stiffness for cooperative binding scales as 1/N, where N is the number of residues.

Population shift scenario is more robust to mutations affecting stiffness.

Abstract

In allosteric proteins, binding a ligand can affect function at a distant location, for example by changing the binding affinity of a substrate at the active site. The induced fit and population shift models, which differ by the assumed number of stable configurations, explain such cooperative binding from a thermodynamic viewpoint. Yet, understanding what mechanical principles constrain these models remains a challenge. Here we provide an empirical study on 34 proteins supporting the idea that allosteric conformational change generally occurs along a soft elastic mode presenting extended regions of high shear. We argue, based on a detailed analysis of how the energy profile along such a mode depends on binding, that in the induced fit scenario there is an optimal stiffness $k_a^*\sim 1/N$ for cooperative binding, where $N$ is the number of residues involved in the allosteric response. We find that the population shift scenario is more robust to mutation affecting stiffness, as binding becomes more and more cooperative with stiffness up to the same characteristic value $k_a^*$, beyond which cooperativity saturates instead of decaying. We confirm numerically these findings in a non-linear mechanical model. Dynamical considerations suggest that a stiffness of order $k_a^*$ is favorable in that scenario as well, supporting that for proper function proteins must evolve a functional elastic mode that is softer as their size increases. In consistency with this view, we find a significant anticorrelation between the stiffness of the allosteric response and protein size in our data set.

Figures

Tables

Table 1. Table 1: Table with the X-ray structures (PDB identifier) used in the analysis. The number of the common residues and of the considered chains are also reported, along with the participation ratio of the allosteric response and the maximal overlap found between the allosteric response and the normal modes. The PDB identifiers are taken from Daily and Gray ( 2007 ) and Mitchell  et al. ( 2016 ) (asterisk).

Protein type	Inactive (PDB)	Active (PDB)	Common residues	Chains	Overlap	Participation ratio
1.   arf6	1E0S	2J5X	160	A	0.27	0.16
2.   cdc42	1AN0	1NF3	183	AB	0.33	0.09
3.   rab11	1OIV	1OIW	162	A	0.32	0.07
4.   rac1	1HH4	1MH1	175	A	0.27	0.07
5.   ras	4Q21	6Q21	164	A	0.44	0.08
6.   rheb	1XTQ	1XTS	165	A	0.35	0.1
7.   rhoA	1FTN	1A2B	173	A	0.24	0.08
8.   sec4	1G16	1G17	152	A	0.23	0.08
9.   IGF-1R	1P40	1K3A	283	A	0.23	0.03
10. met repressor	1CMB	1CMA	204	AB	0.24	0.08
11. tet repressor	2TRT	1QPI	190	A	0.47	0.32
12. glcN-6-P deaminase	1CD5	1HOT	262	A	0.45	0.26
13. EF-Tu	1TUI	1EFT	393	A	0.61	0.30
14. $G_{t\alpha }$	1TAG	1TND	310	A	0.36	0.07
15. ERK2	1ERK	2ERK	347	A	0.19	0.06
16. IRK	1IRK	1IR3	296	A	0.31	0.05
17. lac repressor	1TLF	1EFA	536	AB	0.71	0.48
18. PurR	1DBQ	1WET	271	A	0.71	0.34
19. anthranilate synthase	1I7S	1I7Q	1402	ABCD	0.61	0.45
20. chorismate mutase	2CSM	1CSM	241	A	0.54	0.25
21. FBPase-1	1EYJ	1EYI	323	A	0.46	0.04
22. phosphofructokinase	6PFK	4PFK	315	A	0.32	0.05
23. PTB1B	1T48	1PTY	287	A	0.48	0.02
24. ATCase∗	6AT1	8AT1	2732	ABCDEFGHIJKL	0.64	0.55
25. hemoglobin	4HHB	1HHO	570	ACBD	0.66	0.52
26. NAD-malic enzyme	1QR6	1PJ2	2144	ABCD	0.72	0.57
27. phosphoglycerate DH	1PSD	1YBA	1580	ABCD	0.74	0.53
28. human serum albinum∗	1E78	2BXB	574	A	0.36	0.09
29. fixJ	1DBW	1D5W	238	AB	0.69	0.65
30. DAHP synthase	1KFL	1N8F	1340	ABCD	0.64	0.55
31. SpoIIAA	1H4Y	1H4X	106	A	0.40	0.22
32. CheY	3CHY	1FQW	124	A	0.34	0.18
33. glycogen phosphorylase	1GPB	7GPB	1626	AB	0.53	0.24
34. ATCase	1RAC	1D09	2774	ABCDEFGHIJKL	0.65	0.55

Equations

$$\Delta \Delta E=E_{10}+E_{01}-E_{11}$$

$$E_{00}(x)=k_a||d{R}_a|{|}^{2}f(x)$$

$$E_{00}(x)=k_a||d{R}_a|{|}^2\left(\frac{1}{8}{x}^4-\frac{1}{4}{x}^2+xb\right),$$

$$P=\frac{||d{R}_a|{|}^2}{N{\sum }_{i=1}^N||d{R}_a(i)|{|}^4}.$$

$$k_a=\frac{⟨d{R}_a|H|d{R}_a⟩}{||d{R}_a|{|}^2},$$

$$q_i=\frac{|⟨d{R}_a|v_i⟩|}{\sqrt{⟨d{R}_a|d{R}_a⟩⟨v_i|v_i⟩}}.$$

$$\log \Gamma \equiv \left[\gamma {\log }_{10}(\text{P})+{\log }_{10}(||E_{sh}||)\right]$$

$$V_{ij}=\frac{k}{2}(l_{ij}-l_{ij}^0{)}^2,$$

$$H_{i\mu ,j\nu }=\sum _{kl}{\frac{{\partial }^2{V}_{kl}}{\partial R_{i\mu }\partial R_{j\nu }}|}_{l_{kl}=l_{kl}^0}=\frac{k}{2}\sum _{kl}\frac{\partial l_{kl}}{\partial R_{i\mu }}\frac{\partial l_{kl}}{\partial R_{j\nu }},$$

$$k_a^{\ast }=\frac{⟨d{R}_a|H|d{R}_a⟩}{⟨d{R}_a|d{R}_a⟩}.$$

$${ϵ}_{ab}(\vec{X})=\frac{1}{2}\left(\frac{\partial \vec{x}}{\partial X_a}\cdot \frac{\partial \vec{x}}{\partial X_b}-{\delta }_{ab}\right),$$

$$\Delta {\vec{x}}_{ij}={\overleftrightarrow{\Lambda }}_i\cdot \Delta {\vec{X}}_{ij},$$

$$MSE(i)=\sum _j(\Delta {\vec{x}}_{ij}-{\overleftrightarrow{\Lambda }}_i\cdot \Delta {\vec{X}}_{ij}{)}^2{w}_j(i),$$

$${\overleftrightarrow{\Lambda }}_i=\sum _j\Delta {\vec{x}}_{ij}\Delta {\vec{X}}_{ij}{w}_j(i)\cdot {\left(\sum _j\Delta {\vec{X}}_{ij}\Delta {\vec{X}}_{ij}{w}_j(i)\right)}^{-1},$$

$$\overleftrightarrow{ϵ}(i)=\frac{1}{2}\left({\overleftrightarrow{\Lambda }}_i^t\cdot {\overleftrightarrow{\Lambda }}_i-\overleftrightarrow{\delta }\right),$$

$$E_{sh}(i)=\frac{1}{2}\sum _{l,m=1}^3[{\gamma }_{lm}(i){]}^2,$$

$$\log \Gamma \equiv \left[\gamma {\log }_{10}(\text{P})+{\log }_{10}(||E_{sh}||)\right]$$

$$E_w^{\alpha }=p_0^{\alpha }+p_1^{\alpha }{x}_{\alpha }+p_2^{\alpha }{x}_{\alpha }^2+p_3^{\alpha }{x}_{\alpha }^3+p_4^{\alpha }{x}_{\alpha }^4,$$

$$p_1^{\alpha }+2p_2^{\alpha }{x}_{Ac}+3p_3^{\alpha }{x}_{Ac}^2+4p_4^{\alpha }{x}_{Ac}^3=0,$$

$$p_4^{\alpha }=\frac{1}{2}{k}_w{a}^2,$$

$$p_1^{\alpha }(x_{In}-x_{Ac})+p_2^{\alpha }(x_{In}^2-x_{Ac}^2)+p_3^{\alpha }(x_{In}^3-x_{Ac}^3)$$

$$+p_4^{\alpha }(x_{In}^4-x_{Ac}^4)=k_w{a}^2{b}_w,$$

$$E_s^{\alpha }=\frac{1}{2}k(l_{\alpha }-l_0{)}^2,$$

$$E=\sum _{\alpha }{E}_s^{\alpha }+\sum _{{\alpha }^{\prime }}{E}_w^{{\alpha }^{\prime }}.$$

Full text

**Significance statement: **Many proteins are allosteric: binding a ligand affects their activity at a distant site. Understanding the principles allowing for such an action at a distance is both of fundamental and practical importance. From the thermodynamic viewpoint, two models have been proposed, according to which binding creates a new stable configuration or instead shifts the thermal equilibrium between existing states. We perform a mechanical analysis of these models and show that they are not equally robust to mutations. We argue that in both cases function can most properly occur along a soft elastic mode, whose stiffness decreases rapidly with protein size. We show data on 34 proteins substantiating this result, supporting a new principle for allosteric design.

Many proteins are allosteric: binding a ligand at one or several allosteric sites can regulate function at a distant site, a long-range communication often accompanied by large conformational changes Perutz et al. (1960); Perutz (1970). There is a considerable interest in predicting the amino acids involved in this communication, or “allosteric pathway”, from structure or sequence data Halabi et al. (2009); Amor et al. (2016), since they can be used as targets for drug design Nussinov and Tsai (2013). Yet, understanding the physical principles underlying such action at a distance in proteins remains a challenge Wodak et al. (2019); Thirumalai et al. (2018). From a thermodynamic standpoint, two distinct views have been proposed. In the induced fit scenario, exemplified by the Koshland-Némethy-Filmer (KNF) model Koshland et al. (1966), the protein essentially lies in one single state. The latter changes as binding occurs, leading to a conformational change. In an energy landscape picture such as that of Fig.1B, it corresponds to a displacement of the energy minimum upon binding. By contrast, in the population shift model, exemplified by the Monod-Wyman-Changeux (MWC) model Monod et al. (1965), two states are always present. Their relative stability can change sign upon binding, leading to an average conformational change. Although each of these models presumably applies to various proteins, they do not specify which designs allow for efficient action at a distance, and how robust these designs are to mutations Hopfield (1973).

In several proteins — see below for a systematic study — it has been observed that the allosteric response induced by binding occurs predominantly along one or few vibrational modes Kitao and Go (1999); Bahar et al. (1999); Xu et al. (2003); De Los Rios et al. (2005); Zheng et al. (2006). This result supports that in at least some proteins elasticity  — possibly non-linear —  is an appropriate language to describe allostery (in contrast to intrinsically disordered proteins that may be considered more as liquids than solids, for which the analysis proposed here would not hold). Very recently, there has been a considerable effort to use in-silico evolution Hemery and Rivoire (2015); Tlusty et al. (2017) to study how linear elastic materials can evolve to accomplish an allosteric task Yan et al. (2017); Rocks et al. (2017); Flechsig (2017); Yan et al. (2018); Dutta et al. (2018); Wang et al. (2018); Flechsig and Togashi (2018); Bravi et al. (2018); Rocks et al. (2019). In general, binding a ligand locally distorts the protein, which is modelled by imposing local displacements at some site, generating an extended elastic response that in turn determines fitness (chosen specifically to accomplish a given task). These models fall into the induced fit scenario, since in the framework of linear elasticity there is always a single minimum of energy. A particularly key allosteric function within proteins is the amount of cooperative binding, defined as the change of binding energy of a substrate at the active site caused by binding a ligand at the allosteric site. Materials optimized to display such a cooperativity over long-distances develop a single extended “mechanism” — a soft elastic mode, such as the motion of closing scissors — connecting the two binding sites Yan et al. (2018). It is found that the stiffness $k_{a}$ (i.e. the curvature of the energy) of this mode cannot be too large nor too small for cooperativity to occur, and that optimal design corresponds to $k^{*}_{a}\sim 1/N$ where $N$ is the number of particles in the system. If proteins are nearly optimal, mutations stiffening that mode should thus diminish cooperativity. Yet, these predictions are restricted to strictly linear elasticity, an approximation that presumably does not hold in the regime where most protein operate — certainly not in the population shift scenario.

In this work, we show that the population shift and the induced fit models are very different from a mechanical perspective. In the former case as well, function can be achieved by developing a mechanism or soft elastic mode, but cooperativity, instead of steadily decreasing, saturates to a constant value once the mode stiffness passes some characteristic value, whose scaling with $N$ is $k^{*}_{a}\sim 1/N$. We confirm this prediction in a non-linear elastic model of allostery. This result implies that cooperativity is more robust towards mutations increasing stiffness in the population shift scenario. Yet, displaying a stiffness much larger than $k^{*}_{a}$ implies a very long transition time between the two states, and is presumably prohibited, suggesting the hypothesis that allosteric proteins function with modes presenting a stiffness near the characteristic value $k_{a}^{*}$ in both cases. We test this proposition systematically using X-ray crystallographic data of 34 allosteric proteins. We first confirm that one or a few vibrational modes contribute to the allosteric response, and introduce a new observable establishing that this response presents unusually extended regions of large shear, as found previously for three proteins Mitchell et al. (2016). Next we confirm that the characteristic mode stiffness tends to decrease with the propagation length as we expect from the predicted scaling of $k_{a}^{*}$. Finally, we suggest systematic mutational studies to further test how mechanics constrains allostery.

Geometrical interpretation of mechanical constraints in induced fit allostery

As sketched in Fig.1A, a protein with two binding sites can be unbound (labeled “00”), bound to a single ligand (labeled “01” or “10”) or doubly bound (labeled “11”). We define by $E_{00},\,E_{10},\,E_{01},\,E_{11}$ the energy of the protein in these four situations (corresponding to the minimum energy of their energy landscape), and choose $E_{00}=0$ as the reference energy. Cooperativity is then defined as:

[TABLE]

To simplify notation below, we assume a symmetry between the states 01 and 10, in particular $E_{10}=E_{01}$. Our qualitative conclusions below however remain valid even if this symmetry does not hold.

Consider a protein displaying cooperativity thanks to the presence of a soft elastic mode. Let us denote by $x$ a variable indicating the motion along that mode (see our numerical model below for a concrete example in the context of a shear design), which varies from zero to unity as the protein undergoes its allosteric response. The energy profile $E_{00}(x)$ of the unbound state follows:

[TABLE]

where $||dR_{a}||^{2}=\sum_{i=1...N}||d{\vec{R}_{i}}||^{2}$ is the norm of the allosteric response $|dR_{a}\rangle\equiv\{d{\vec{R}_{a}(i)}\}$ and $d{\vec{R}_{a}(i)}$ is the vector displacement of the amino acid $i$. We expect that if the protein has $N$ residues, $||dR_{a}||^{2}$ is of order $Na^{2}$, where $a$ is the inter-atomic distance. We will confirm empirically the dependence of $||dR_{a}||^{2}$ with $N$ below. $k_{a}$ is the mode stiffness and $f(x)$ some function of order unity such that $f^{\prime\prime}(x=0)=1$. For a purely linear elastic material, $f(x)=x^{2}/2$. More generally for the induced fit scenario, $f(x)$ presents a single minimum, as illustrated in Fig.1B. In Yan et al. (2018) it was argued that if linear elasticity applies, then there is an optimal stiffness $k_{a}^{*}$ for cooperativity. We present a geometrical interpretation of this result, which extends it to the induced fit scenario in general, and will be useful to explain why the population shift scenario behaves differently. Local chemical interactions which lead to identical binding energies for the inactive and active states do not affect cooperativity, and thus need not be considered here. For concreteness, we assume that after some motion $x_{0}=1$ along that mode, the protein shape can accommodate perfectly both ligands without deforming. It implies that the energy profile of the bound states $E_{10}(x)$ and $E_{11}(x)$ satisfy $E_{10}(x_{0})=E_{11}(x_{0})=E_{00}(x_{0})$, as pictured in Fig.1B. However as $x$ departs from $x_{0}$, the protein shape does not match the ligands, imposing an elastic strain at the binding sites leading to an increase of elastic energy in the protein, that will trigger motion along other elastic modes. Assuming that the ligands are rigid and that each binding site involves on the order $n_{0}$ atoms that move by a distance $a$ as the protein undergoes its allosteric response, we have $E_{01}(x)-E_{00}(x)=n_{0}ka^{2}g(x-x_{0})$ where $k$ characterizes the stiffness of inter-atomic interactions and is large in comparison with the one of soft modes, i.e. $k\gg k_{a}$. Here $g$ is some dimensionless function vanishing quadratically in zero — but possibly non-linear at large arguments.

$E_{11}$, $E_{01}$ can be computed as the minimum of the curves $E_{11}(x)$ and $E_{01}(x)$ respectively, from which the cooperativity $\Delta\Delta E$ is readily computed. Two extreme cases occur, illustrated in Fig.1B:

(i) if $k_{a}N\ll n_{0}k$ (Fig.1B.I), as $x$ moves away from $x_{0}$, the elastic energy induced by binding $n_{0}ka^{2}g(x-x_{0})$ is very significant in comparison to the mode energy $E_{00}(x)\sim k_{a}Na^{2}f(x)$. Thus both $E_{10}(x)$ and $E_{11}(x)$ have a sharp minimum near $x_{0}$, with $E_{11}\approx E_{10}\approx k_{a}||dR_{a}||^{2}f(1)$. Thus $\Delta\Delta E=E_{10}+E_{01}-E_{11}\approx E_{10}\approx k_{a}||dR_{a}||^{2}f(1)$, which vanishes as $k_{a}\rightarrow 0$.

(ii) If $k_{a}N\gg n_{0}k$ is very large (Fig.1B.III), $n_{0}ka^{2}g(x-x_{0})$ is small in comparison to $E_{00}(x)$: $E_{00}(x)$, $E_{11}(x)$ and $E_{01}(x)$ are very close to each other, and must thus all present a minimum near $x=0$. Thus binding does not trigger motion along the soft mode, whose presence is useless. No extended modes couples the two binding sites and $E_{11}\approx E_{01}+E_{10}$, leading to $\Delta\Delta E\rightarrow 0$ as $k_{a}\rightarrow\infty$.

(iii) Optimal cooperativity is thus found at some intermediary $k_{a}^{*}\sim n_{0}k/N$, corresponding to Fig.1B.II. Note that the present argument for an optimal $k_{a}^{*}$ does not require the energy profile $f(x)$ to be an exact parabola, as long as it is monotonically growing in both directions around its minimum.

Mechanical aspects of the population shift model

MWC model

We recall some aspects of the MWC model. To simplify notation, we consider that the protein displays two symmetric binding sites as illustrated in Fig.1A. The protein is assumed to lie in two possible distinct configurations, “Inactive” ($In$) and “Active” ($Ac$). In the absence of binding, we take the energy of the inactive state as our reference (i.e. $E_{00}^{In}=0$) and denote the energy of the active state $E_{00}^{Ac}=E_{0}$. We assume that the active configuration has a well-suited geometry to bind each ligand, thus no elastic energy is spent for binding and $E_{10}^{Ac}=E_{01}^{Ac}=E_{11}^{Ac}=E_{0}$. By contrast, we assume that in the inactive state binding costs some energy $\Delta E$, leading to $E_{10}^{In}=E_{01}^{In}=\Delta E$ and $E_{11}^{In}=2\Delta E$. The last assumption of additivity of binding energies within a given configuration, i.e. for a frozen mode amplitude, is expected to be accurate if the two binding sites are distant enough, since no elastic coupling between them is expected in that case.

For each binding situation, the configuration (inactive or active) chosen is the one with the smallest energy, e.g. $E_{01}=\min(E^{In}_{01},E^{Ac}_{01})$. Computing the cooperativity one finds three different cases:

(i) if $E_{0}<\Delta E$, then the binding of one ligand is sufficient to drive the system in the active state, implying $E_{10}=E_{11}=E_{0}$ and $\Delta\Delta E=E_{0}$;

(ii) if $\Delta E<E_{0}<2\Delta E$, then the binding of one ligand is not sufficient and two ligands have to bind in order to drive the system in the active state. Consequently, $\Delta\Delta E=2\Delta E-E_{0}$;

(iii) if $E_{0}>2\Delta E$, then $\Delta\Delta E=0$ since the system stays in the inactive state even if two ligands are bound. The sketch of this behavior is shown in Fig.2, illustrating that the maximum cooperativity is found for $E_{0}=\Delta E$.

Mechanical consideration on the MWC model

Our observations (see the empirical section below) indicate that a significant fraction of allosteric proteins operate mainly along one normal mode of the elastic energy, supporting the idea that in these cases a favored path connects the inactive and active configurations. We expand the energy (in the absence of ligand) in terms of the motion $x$ along that path (in this two-states case, it is more convenient to chose a coordinate $x$ varying between $-1$ and $+1$ as the protein undergoes its allosteric response). Note that if non-linearities are present, this path is not along a single linear mode, but it bends in phase space. Our analysis below holds independently of such bending. We keep the minimal number of non-linear terms that allow to display two states (a polynomial of four degrees has five parameters, three of them can be fixed by redefining the reference energy, and changing the definition of $x$ by both a multiplicative and additive constants):

[TABLE]

where $b$ is a parameter reflecting how the energy profile is tilted towards the inactive state, $k_{a}$ characterizes the stiffness of the mode and $||dR_{a}||^{2}$ is the square norm of the allosteric response.

A typical profile following Eq.3 is shown in Fig.3A. We denote by inactive the lowest of the two minima and the other one by active. Note that (i) in the case $b=0$, Eq.3 describes two identical minima at $x=\pm 1$ of stiffness $k_{a}$, separated by an energy barrier $E_{b}=k_{a}||dR_{a}||^{2}/8$. (ii) When the parameter $b$ is positive, we have $E_{Ac}>E_{In}$ up to $b=b_{c}=1/(3\sqrt{3})$ where the active state becomes unstable and only a single stable state is left. (iii) At fixed $k_{a}$, the energy difference between the two states $E_{0}$ is maximal at $b=b_{c}$, where one finds $E_{0}=\frac{3}{8}k_{a}||dR_{a}||^{2}$. (iv) For small $b$, $E_{0}\simeq 2k_{a}||dR_{a}||^{2}b$. In what follows we focus on the case $b<b_{c}$ where the population shift model lies.

Next, we assume that the active configuration matches the shape of both ligands, so that binding events in that state cost no energy. However, by moving away from the configuration the shapes of the protein and ligands do not match anymore and the protein needs to deform elastically near the binding site. Again assuming that each binding site involves of order $n_{0}$ atoms, which move by a distance of order $a$ from the active to the inactive state, we have a binding energy $E_{10}(x)-E_{00}(x)=n_{0}ka^{2}g(x-x_{Ac})$ where $x_{Ac}\approx 1$ is the location of the active state along the path. This energy is exemplified by the difference between the blue and black curves in Fig.3A. Thus the binding energy in the inactive state follows $\Delta E=n_{0}ka^{2}g(x_{In}-x_{Ac})\approx 4n_{0}ka^{2}$, which is independent of the protein size and mode stiffness $k_{a}$. As explained before, if the two binding sites are distant enough, for a given mode amplitude the elastic costs of binding will simply add up: $E_{11}(x)-E_{00}(x)=2n_{0}ka^{2}g(x-x_{Ac})$ as shown in red in Fig.3A.

How the mode stiffness constrains cooperativity

To quantify this constraint we define the maximal cooperativity over all possible tilts $b$ given $k_{a}$: $\Delta\Delta E^{*}\equiv\max_{b}\left\{\Delta\Delta E(k_{a},b)\right\}$. We find two regimes:

(i) If $Nk_{a}\ll n_{0}k$, then the elastic costs associated with binding are very large compared to $E_{0}$. Both $E_{10}(x)$ and $E_{11}(x)$ are peaked close to $x_{Ac}$, as illustrated in Fig. 3B.a. Thus $E_{0}<\Delta E$, implying $\Delta\Delta E=E_{0}$ according to Fig. 2, which is maximized at $b=b_{c}$ leading to $\Delta\Delta E=\frac{3}{8}k_{a}||dR_{a}||^{2}\sim k_{a}Na^{2}$. Thus $\Delta\Delta E$ vanishes linearly at small $k_{a}$, as illustrated in Fig.3B. This result is qualitatively similar to the induced fit case, for which $\Delta\Delta E$ also vanishes linearly as shown in Fig.1.B.

(ii) If $Nk_{a}\gg n_{0}k$, then the elastic cost is very small in comparison to $E_{0}$: $E_{00}(x)$, $E_{10}(x)$ and $E_{11}(x)$ are almost identical as illustrated in Fig.3B.c. In that regime, cooperativity is optimized by chosing a small tilt fixing $E_{0}$ to $\Delta E$ according to Fig. 2, implying $\Delta\Delta E=\Delta E\sim n_{0}ka^{2}$ which is independent of $k_{a}$. This plateau behaviour is represented in Fig. 3B, and appears at $k^{*}_{a}\sim n_{0}k/N$.

This result represents a fundamental difference with the induced fit case, for which a large stiffness destroys cooperativity. Indeed in the induced fit scenario a large stiffness implies that the minimal energy is always found for $x\approx 0$ as illustrated in Fig. 1B.III, implying that binding does not move the protein along that mode, which is thus useless. This state of affairs is ultimately a geometric necessity stemming from the fact that the three curves $E_{00}(x)$, $E_{10}(x)$ and $E_{11}(x)$ must be very close to each other in that regime, and each present a single minimum. Consequently the positions of these minima must be very similar in the three cases, leading to $x\approx 0$ independently of binding ligands or not. This geometric necessity vanishes as soon as two minima are present.

Note that although $\Delta\Delta E$ asymptotes to a constant for $k_{a}\gg k_{a}^{*}$, the barrier $E_{b}$ between the inactive and active states grows linearly with $k_{a}$ in that limit. Large barriers would lead to undesirably slow transition rate between states, thus we expect that in practice $k_{a}$ lies reasonably close to $k_{a}^{*}$.

A mechanical model for population shift allostery

We seek to model that (i) the allosteric response often takes place mainly along a single vibrational mode. (ii) Various architectures can lead to allostery, including the well-known shear Mitchell et al. (2016); Gerstein et al. (1994) or hinge Xu et al. (2003) designs and others not falling in these categories Goodey and Benkovic (2008); McLaughlin Jr et al. (2012). Such a diversity is also found in in silico evolution schemes Yan et al. (2018). Yet, such synthetic architectures always present soft extended regions where most of the strain (i.e. relative motion) is located. Such an observation was made in a few proteins Mitchell et al. (2016) and will be generalized below. For proteins presenting two stable configurations, we expect these regions to present two possible ways of locally stacking amino acids well.

As an illustration, we consider the shear design in which the protein presents two three-dimensional rigid regions connected by a soft planar layer that can easily deform. The rigid regions consist of harmonic springs of stiffness $k$, shown in blue in Fig.4A. The soft layer consists of anharmonic springs (shown in red in Fig.4A), whose energy vs. extension curve is non-linear, chosen to have the same form as Eq.3. These non-linear springs have a characteristic stiffness $k_{w}$, and present two stable extensions at which they exert no force, whose relative energy is controlled by some bias $b_{w}$. These stable extensions are chosen such that two states of the protein as a whole, the inactive and active states shown in Fig.4A, present no contact forces by construction, and are thus local minima of the energy. Finally the protein presents two binding sites, at its top and bottom. At each site, binding imposes that the distance between two nodes (indicated by black lines in Fig.4A) equals the distance it naturally presents in the active state, thus favoring it. The details about the construction of the microscopic model are discussed in Supplementary Information.

We can now numerically compute the energy profiles $E_{00}(x)$, $E_{01}(x)$ and $E_{11}(x)$ as a function of the motion along the shear mode $x$, by imposing a shear displacement (i.e. a value of $x$) and letting the entire elastic energy of the material relax, except for that mode. The insets of Fig. 4 show our results. The value of the mode stiffness $k_{a}$ can be extracted from fitting Eq.3 to $E_{00}(x)$ and measuring the displacement norm $||dR_{a}||$, and can be increased by increasing $k_{w}$. From $E_{00}(x)$ one readily computes the energy difference $E_{0}$ between the inactive and active states, which can be increased by monitoring the microscopic bias $b_{w}$. From the minima of $E_{00}(x)$, $E_{01}(x)$ and $E_{11}(x)$ one immediately extracts the binding costs $\Delta E$ and $\Delta\Delta E$.

Fig.4B shows $\Delta\Delta E$ for two values of $k_{a}$, as the energy difference $E_{0}$ is increased. For large $k_{a}$, $\Delta\Delta E$ passes through a maximum $\Delta\Delta E^{*}=\Delta E$, whereas for small $k_{a}$, $\Delta\Delta E$ is smaller, and its maximum is fixed by the maximal achievable energy difference $E_{0}$. $\Delta\Delta E^{*}$ as a function of $k_{a}$ is shown for different system sizes $N$ in Fig.4C, confirming the presence of two regimes with a cross-over at some $k_{a}^{*}\sim 1/N$ as shown in Fig.4D. Overall, these observations validate our theoretical predictions on the dependence of $k_{a}^{*}$ with the number of residues, and on how the stiffness $k_{a}$ qualitatively affects cooperativity.

Empirical study of 34 allosteric proteins

Allosteric Response

We identify a set of 34 allosteric proteins in the Protein Data Bank (PDB), for which both the active (ligand bound at the allosteric site) and inactive (no ligand bound at the allosteric site) crystalline X-ray structures are available. Their PDB identifiers are taken from Daily and Gray (2007); Mitchell et al. (2016) and reported in Supplementary Information. The set is diverse in functionality, and includes enzymes (13), G-proteins (10), Kinases (3), response regulators (3), DNA-binding proteins (4) and the Human Serum Albumin, among which 12 protein complexes are present. We can thus estimate the allosteric response $|dR_{a}\rangle$ as the displacement field between the inactive and active structures (after having aligned them via the software Pymol 2.1.1 Schrödinger, LLC (2015)). Here we focus on the motion of the $N$ amino-acids, located by the position of their $\alpha$-carbon. As an illustration, the allosteric response of a given protein (the elongation factor Tu) is shown in black arrows in Fig. 5A.v From the allosteric response $|dR_{a}\rangle$, one can readily estimate: (i) the magnitude of the displacement $||dR_{a}||^{2}$. (ii) The fraction of the protein involved in the response. For any displacement field, this fraction is usually estimated via the participation ratio Bell and Dean (1970):

[TABLE]

(iii) A measure of how much relative displacement takes place around atom $i$. Following Mitchell et al. (2016) we consider the shear pseudo-energy $E_{sh}(i)$ quantifying the amount of strain — essentially a measure of the relative displacement between adjacent atoms — at residue $i$, whose precise definition is given in Supplementary Information. $E_{sh}(i)=0$ indicates that the protein moves as a rigid body near atom $i$, and by contrast is large where atoms slide rapidly past each other. $E_{sh}(i)$ is shown in color in Fig. 5A for the protein Tu, illustrating that two parts of the protein are rigidly moving (and counter-rotate), while the central region displays significant pseudo-energy $E_{sh}(i)$, which is reminiscent of a hinge design.

Elastic Networks Analysis

To estimate protein elasticity we use elastic network models (ENM) Atilgan et al. (2001), in which harmonic springs of identical stiffness are placed between all $N$ $\alpha$-carbons laying below a chosen cutoff radius $R_{c}\in[8-12]$ Å. ENM is obviously a crude approximation of real atomic interactions. Yet, its simplicity allows for the systematic study of various proteins, and it has been successful in capturing normal modes relevant for the function of some proteins. The dependence of the results on the value of $R_{c}$ are indicated by error bars in Fig. 5B and discussed in Supplementary Information for Fig. 6B. This procedure defines an elastic energy from which the matrix of the second derivatives, i.e. the Hessian matrix $H$, can be computed. From $H$, one can readily estimate the stiffness $k_{a}$ of the allosteric response as the curvature of the elastic energy in that direction:

[TABLE]

Finally, the eigenvectors of $H$ define the $3N$ vibrational modes of the protein $\{|v_{i}\rangle\}_{i=1\dots 3N}$. Following Kitao and Go (1999); Bahar et al. (1999); Xu et al. (2003); De Los Rios et al. (2005); Zheng et al. (2006), the overalp $q_{i}$ between the allosteric response and mode $i$ characterizes their similarity ($q_{i}=1$ implies that there are identical):

[TABLE]

Geometry of the allosteric response

Fig. 5B reports the maximal overlap $q^{*}\equiv\max_{i}q_{i}$, as a function of the participation ratio $P$. Our observation indicate that $q^{*}$ is in general large (in half of the case larger than $0.45$), supporting further that allostery indeed occurs mainly along one mode Kitao and Go (1999); Bahar et al. (1999); Xu et al. (2003); De Los Rios et al. (2005); Zheng et al. (2006). Interestingly, this effect is stronger when most sites of the protein are involved in the allosteric response ($P$ large). Interestingly, some protein complexes (for example ATCase, see the values of $q^{*}$ in Table 1 of Supplementary Material) also display a large projection of their allosteric response on a single mode, supporting that mechanical considerations can be valuable in such cases as well.

We now provide systematic data supporting that the allosteric response presents extended regions of large shear energy Mitchell et al. (2016). More specifically, we argue that while some vibrational modes can present significant shear (e.g. localized modes capturing the motion of dandling loops) and other can be extended (such as plane-wave-like modes), the allosteric response is unique in presenting both aspects, thus revealing a specific design principle. To quantify this effect, we introduce the quantity, that can be defined on any displacement field:

[TABLE]

where $||E_{sh}||$ is the total magnitude of the shear energy, i.e. $||E_{sh}||=(\sum_{i}E_{sh}(i)^{2})^{1/2}$. $\log\Gamma$ is large if the displacement is extended and if the shear energy is large. The factor $\gamma$ characterizes the trade-off between these two features. Here we choose $\gamma=3.5$ reflecting the fact that for vibrational modes, we find that $P$ varies about 3.5 times less in relative terms than $||E_{sh}||$ as shown in Supplementary Information. Thus for $\gamma=3.5$, the spatial extension and the amount of shear equally affect $\log\Gamma$. Fig. 5C shows $\log\Gamma$ averaged over the $17$ proteins with $q^{*}>45\%$ for the allosteric response (yellow line), the mode with maximum overlap (blue line) and the first $75$ vibrational modes (having subtracted the one with largest overlap) as function of the mode rank. We find that $\Gamma$ is typically 160 times larger for the allosteric response than for vibrational modes, a very significant difference underlying the specific geometry of the allosteric response.

Scaling of response stiffness $k_{a}$ with protein size

We can now test our conjecture that the allosteric stiffness $k_{a}$ is close to $k_{a}^{*}\sim 1/N$ where cooperativity saturates, which implies in particular an anti-correlation between $k_{a}$ and protein size. In our theoretical estimate of $k_{a}^{*}$, we have assumed that the allosteric response magnitude is linear in the protein size, i.e. $||dR_{a}||^{2}\sim N$. It is a natural assumption since the larger the protein, the more likely its response involves many residues. The relationship between these two quantities is tested in 34 proteins in Fig.6A. We indeed find a strong correlation between the logarithms of $||dR_{a}||$ and $N$ (Pearson coefficient $r=0.76$ with p-value $p=1.64\times 10^{-7}$. Pearson coefficients are computed on the logarithmic values via the Matlab R2017b function corr). Overall, data are consistent with our assumption of proportionality.

Finally, we plot the allosteric response stiffness $k_{a}$ measured according to Eq.S3 in terms of $N$ for all proteins in Fig. 6B. A key finding is the fair anti-correlation between the logarithms of these two quantities (Pearson coefficient $r=-0.64$ with p-value $p=4.00\times 10^{-5}$), supporting the idea that larger allosteric proteins need to evolve a softer elastic mode to accomplish function, as expected from our analysis. The signal that we observe is encouraging, especially given the possible sources of noise in the analysis, among which, in particular, the use of ENM to build the Hessian matrix.

Conclusion

We have provided systematic evidence that the allosteric response occurs along one soft elastic mode, and we have introduced a novel observable $\Gamma$ to establish that this response generally displays unusually extended regions of high shear strain. These observations support that for many proteins, elasticity is a useful starting point to describe allostery. We have revisited the two classical thermodynamic models of allostery from this perspective, and provided a detailed study of how the energy profile along the soft mode evolves with binding. We find that induced fit and population shift models qualitatively differ. In the induced fit model, there is an optimal stiffness $k_{a}^{*}\sim 1/N$ associated to that mode beyond which the cooperative binding energy eventually decreases to zero. The population shift model is more robust to mutations affecting stiffness, and $k_{a}^{*}\sim 1/N$ simply marks a cross-over beyond which cooperativity saturates and the transition time between configurations rapidly explodes. We introduced a novel non-linear elastic model for allostery supporting these views. Our key result is that proper function is achieved if proteins evolve an elastic mode whose softness must rapidly decrease with size, a prediction supported by the anti-correlation observed between these quantities.

Systematic mutation scan on one single protein, in which binding assays to measure cooperativity are combined with single molecule experiments or ultrafast laser pulses to estimate the stiffness of the allosteric response, would be extremely useful to test the predicted relationship between these quantities. Molecular dynamics experiments could further test how the energy profile along the soft elastic mode evolves with binding. Elucidating such an interplay between thermodynamics and mechanics in proteins would be valuable in a variety of tasks, including de novo protein design and the discovery of novel allosteric pathways.

Acknowledgements.

We thank P. Barth, B. Bravi, P. De Los Rios, D. Malinverni, R. B. Phillips and M. Popović for discussions. R.R. is supported by the Swiss National Science Foundation under Grant No. 200021-165509/1. M.W. thanks the Swiss National Science Foundation (Grant No. 200021-165509) and the Simons Foundation (Grant #454953 Matthieu Wyart). This material is based upon work performed using computational resources supported by the “Center for Scientific Computing at UCSB” and NSF Grant CNS-0960316, by the “High Performance Computing at NY” and by the “High Performance Computing at EPFL”.

I Supplemental Material

A. Mechanics of allosteric proteins from X-ray structural data

Dataset

We consider 34 allosteric proteins with both inactive and active X-ray structures available Daily and Gray (2007); Mitchell et al. (2016). The Protein Data Bank (PDB) entries of these proteins are listed in Table 1; the number of common residues between the inactive and active structures and the number of chains that we consider are also reported in the table. We remove the first two and last two residues for every structure as they correspond to the fluctuating starting and ending points of the protein chain.

Elastic network model

From the X-ray structures we compute the Hessian matrix $H$ using the anisotropic network model introduced in Atilgan et al. (2001). We consider the positions of residues by looking only at alpha-carbon $C_{\alpha}$ atoms; ${\bf R}_{i}$ for alpha-carbon $i$. The model assumes a harmonic interaction between two residues at distance smaller than a cutoff distance $R_{c}$,

[TABLE]

where $k$ is the spring constant (fixed to unity), $l_{ij}=||{\bf R}_{i}-{\bf R}_{j}||$ is the distance between residues $i$ and $j$, and $l_{ij}^{0}<R_{c}$ is the distance at equilibrium. Building the Hessian is then straightforward by taking second derivatives of the potential $V_{ij}$ with respect to the coordinates of residues and evaluated at equilibrium $l_{ij}=l_{ij}^{0}$,

[TABLE]

where $\mu,\nu=1,2,3$ label the spatial dimension of the atoms.

The definition of a cutoff distance is an empirical fitting parameter. In practice, we computed Hessian matrices at nine cutoff values equi-distributed in the range $[8-12]$ Å. We define the stiffness of the allosteric response $|dR_{a}\rangle$ as

[TABLE]

This quantity changes systematically with $R_{c}$ as shown in Fig.S1A: by increasing $R_{c}$, $k_{a}^{*}$ also increases, which leads to all points in Fig. S1B moving together in the same direction when changing $R_{c}$. In the main text we reported the mean value of $k_{a}^{*}$ over the nine choices of $R_{c}$, and in Fig. S1B we show the range of variability of $k_{a}^{*}$ when $R_{c}\in[8-12]$ Å.

Computing the local strain tensor in a protein

In a continuous medium, a motion maps a point $\vec{X}$ in the reference configuration to a new point ${\vec{x}}$ in the current configuration, the strain tensor of the motion can thus be computed as

[TABLE]

where $a$, $b$ labels the spatial dimension.

In a discrete medium like proteins which are a collection of atoms (or residues as we consider), computing the partial derivative $\overset{\leftrightarrow}{\Lambda}=\partial\vec{x}/\partial\vec{X}$ at residue $i$ is not straightforward. Ideally, for any neighboring residue $j$ close enough in space

[TABLE]

where $\Delta\vec{X}_{ij}=\vec{R}_{i0}-\vec{R}_{j0}$ and $\Delta\vec{x}_{ij}=\vec{R}_{i}-\vec{R}_{j}$ in our setting, where $\vec{R}_{i}$ is the position of residue $i$ taken from the X-ray structure. We have $n_{b}$ number of such equations for $\overset{\leftrightarrow}{\Lambda}_{i}$ when $n_{b}$ neighbors are considered. So $\overset{\leftrightarrow}{\Lambda}_{i}$ are usually over-determined when we consider all neighbors below a certain cutoff distance $R_{c}$ (we choose $R_{c1}=8.5$ Å for first nearest neighbors and $R_{c2}=10.5$ Å for second nearest neighbors). Instead of solving Eq. (S5), we define a mean squared error function Gullett et al. (2007)

[TABLE]

where we have kept a weight function $w_{j}(i)$ of node $j$ contribution to $i$ in general. Specifically, we set as in Mitchell et al. (2016) $w_{j}(i)=1$ for all nearest neighbors to $i$ ($R_{ij}<R_{c1})$, $w_{j}(i)=1-\dfrac{R_{ij}-R_{c1}}{R_{c2}-R_{c1}}$ for $R_{c1}<R_{ij}<R_{c2}$ and $w_{j}=0$ otherwise. By minimizing the mean squared error with respect to $\overset{\leftrightarrow}{\Lambda}_{i}$, we have

[TABLE]

and

[TABLE]

where $\overset{\leftrightarrow}{\delta}$ is the identity tensor.

The shear pseudo-energy Mitchell et al. (2016), a vector field whose components contain a measure of the relative motion of each residue, can be defined from the strain tensor $\overset{\leftrightarrow}{\epsilon}(i)$ computed above

[TABLE]

where the local shear tensor $\overset{\leftrightarrow}{\gamma}(i)=\overset{\leftrightarrow}{\epsilon}(i)-(1/3)\text{Tr}[\overset{\leftrightarrow}{\epsilon}(i)]\overset{\leftrightarrow}{\delta}$ depends only on the displacement between the two conformations via the strain tensor $\overset{\leftrightarrow}{\epsilon}(i)$ and $\overset{\leftrightarrow}{\delta}$.

Empirical definition of $\Gamma$

In the main text we introduced the observable $\Gamma$ combining information on the participation ratio (P) and the amount of shear pseudo-energy ($||E_{sh}||$) in a given mode as

[TABLE]

with $\gamma=3.5$. In Fig. S2.A we show the range of values of participation ratio (dashed line) compared with the one of shear pseudo-energy (continuous line).

The measured relative variation from the two curves is $\gamma^{*}\simeq 3.46$ and we choose $\gamma=3.5$ in the definition of $\Gamma$. Fig. S2.B is an alternative way of representing the result discussed in the main text. We consider a scatter plot of $\log\Gamma$ and participation ratio for the three quantities considered (modes with maximum overlap, other modes and allosteric response) and we see that the allosteric response although displaying a large participation ratio it has also a large $\Gamma$, meaning that it involves also a considerable amount of shear.

B. Microscopic model

To construct a microscopic model with two global minima that correspond to the “Inactive” ($In$) and “Active” ($Ac$) states, we need to impose that all nonlinear springs to stay exactly at one of their minimum when the network is in either $In$ or $Ac$ state.

This can be achieved by having each nonlinear spring $\alpha$ following a quartic potential,

[TABLE]

where $x_{\alpha}=l_{\alpha}/a$, with $l_{\alpha}$ the distance of two particles connected by the interaction $\alpha$ and $a$ a parameter capturing the typical local deformation of adjacent particles between $In$ and $Ac$ states. The above compromise can thus be satisfied by choosing the right prefactors $p_{0}^{\alpha}$ to $p_{4}^{\alpha}$ Given a nonlinear spring in the $In$ and $Ac$ states with rest length, respectively, $x_{In}$ and $x_{Ac}$, the springs are in a local minimum, if (i)

[TABLE]

and (ii) same for $x_{In}$. Only the relative value of energy is important, (iii) we can set all $p_{0}^{\alpha}$ to zero. (iv) We assume the stiffnesses of the nonlinear springs to be the same,

[TABLE]

where $k_{w}$ is the parameter capturing the stiffness of nonlinear interactions that has a unit of spring stiffness. (v) Finally, we define

[TABLE]

where $b_{w}$ is a parameter capturing the energy difference between the $Ac$ and $In$ states. We now have five equations (i-v) for the five parameters ($p_{0}^{\alpha}$ to $p_{4}^{\alpha}$) for each of the nonlinear spring $\alpha$. Hence, we can define the quartic potentials for all nonlinear springs.

Specifically, we consider an elastic network embedded in a Face-Centered-Cubic (FCC) lattice in three dimensions (3D) as shown in Fig. S3-A,B. The $In$ and $Ac$ states are connected by a shear mode, which is enabled by rotations of two rigid blocks on the left and right against each other, as shown in Fig. S3-C. Such a model can be realized by having harmonic springs with spring stiffness $k$ in the two rigid blocks and the nonlinear springs $k_{w}\ll k$ at the boundary of the two.

The elastic energy of a linear (strong) spring $\alpha$ is then

[TABLE]

where $l_{0}$ is the rest length of particle separation on the lattice. While, the energy of a nonlinear (weak) spring $E_{w}^{\alpha^{\prime}}$ obeys Eq. (S10) with parameters solved by the particle distances in two states. The total energy of the system is then obtained by summing the energies over all two kinds of springs,

[TABLE]