Principles for optimal cooperativity in allosteric materials

Allosteric proteins transmit a mechanical signal induced by binding a ligand. However, understanding the nature of the information transmitted and the architectures optimizing such transmission remains a challenge. Here we show using an {\it in-silico} evolution scheme and theoretical arguments that architectures optimized to be cooperative, which propagate efficiently energy, {qualitatively} differ from previously investigated materials optimized to propagate strain. Although we observe a large diversity of functioning cooperative architectures (including shear, hinge and twist designs), they all obey the same principle {of displaying a {\it mechanism}, i.e. an extended {soft} mode}. We show that its optimal frequency decreases with the spatial extension $L$ of the system as $L^{-d/2}$, where $d$ is the spatial dimension. For these optimal designs, cooperativity decays logarithmically with $L$ for $d=2$ and does not decay for $d=3$. Overall our approach leads to a natural explanation for several observations in allosteric proteins, and { indicates an experimental path to test if allosteric proteins lie close to optimality}.


INTRODUCTION
Many proteins are allosteric: binding a ligand at an allosteric site can affect the properties of a distant active site, sometimes located on the other side of the protein [1,2]. Predicting the existence of such allosteric pathways from protein structure alone would be of great interest [3,4], since they can be used as targets for drug design [5]. Solving this challenge requires to make progress on both physical and biological questions. First, how can such disordered materials [6] be designed to carry mechanical information specifically over long distances? Are there fundamental limits to what can be achieved? Second, what are allosteric pathways really optimized for? What kind of elastic information do they carry? A physical theory of allostery should address these points. It should also explain the following empirical facts: (i) Some allosteric proteins [7], including hemoglobin [8,9], essentially function as hinges, while others display a "shear" design where two rigid parts are connected by a weak plane [10]. This classification is however not exhaustive, as in various cases the response to binding a ligand cannot be described in term of a simple shear or hinge motion [11][12][13]. (ii) The response to binding often corresponds mostly to motion along few soft normal modes of the protein [14,15]. These modes tend to be conserved during evolution [16]. (iii) In some cases the allosteric functional effect at the active site is significant while the physical mean displacement induced by binding the ligand is small. It has been proposed that for these proteins binding can affect how particles near the active site fluctuate around their mean position, while changing little the latter [17][18][19][20].
Recently, allostery was investigated using in silico evolution schemes where a system evolves to perform a given function [21][22][23][24][25]. Most relevant here are schemes developed to solve inverse elastic problems [23][24][25], in a spirit similar to topology optimization used in engineering to design functional tools from compliant materials [26][27][28]. The task studied in [23][24][25] was to design a material whose response to a specific local strain applied on one of its sides (the allosteric site) leads to a displacement whose geometry is prescribed on the opposite side (the active site). Under broad conditions these algorithms find solutions that achieve such "geometric" tasks essentially perfectly. The corresponding architectures turn out to have surprising properties: their response almost vanishes in the bulk of the material and reappears near the active site [23]. This amplification of the elastic signal is caused by the emergence of a powerful lever, made of a soft elastic region surrounding the active site, where the system is just constrained enough to act as a solid [23,29]. Although there is great interest in finding whether such architectures exist in nature, an intriguing aspect of this approach is that it does not generate the well-known allosteric architectures such as the simple shear and hinge designs, in which the response remains of similar magnitude between the allosteric and the active sites.
Here we show that a simple modification of the task, where materials are optimized such that the binding at the allosteric site lowers the binding energy of another ligand at the active site, leads to different design principles. In the context of proteins, this task corresponds to maximizing the cooperativity of binding two ligands, a central feature of various allosteric proteins [1]. We find that there is a zoology of architectures achieving such cooperativity, but they always display a stiff structure (embedded in a softer elastic matrix) with a single very soft extended elastic mode or "mechanism". We lay out the principles behind such designs, and show theoretically that the soft mode frequency should be neither too large nor too small to optimize function: its optimal value decreases with the material size, and scales as L −d/2 in spatial dimension d. We prove that cooperativity then decays as ln −1 (L/c) for d = 2 and is even independent of L in larger spatial dimensions d ≥ 3, where L is the linear extension of the system and c the length scale on which binding takes place. This result is very different from a normal continuous elastic medium where cooperativity rapidly decays with distance as L −d . Overall the classification we provide leads to a natural explanation for the key aspects of allostery described in (i,ii,iii). It also shows that a path of large strain values connecting the allosteric and active site induced by binding is not necessary for cooperativity to occur, and it makes further testable predictions, including the locations where a shear or hinge design would be mostly affected by a mutation and conserved during evolution.

In-silico Evolution Scheme
Elastic networks: To model allosteric materials we consider elastic networks, often used to describe proteins [14][15][16]. Specifically, N = L d nodes are located on a lattice (slightly distorted periodically to avoid straight lines as discussed in Supplemental Material Section A and [30,31]), and among all N c links of nearest nodes, a subset of N s pairs are connected by harmonic springs of stiffness k = 1 , as indicated by lines in Fig. 1. We declare that σ α = 1 if a spring is present in the link α and σ α = 0 otherwise. Thus the network is entirely described by a connection vector |σ made of zeros and ones, whose dimension is the number of links N c . We define the average coordination number z ≡ 2N s /N and average connectionσ = N s /N c and keep them fixed during evolution. We find that our results do not depend qualitatively on z as long as z > z c = 2d, the rigidity limit derived by Maxwell [32].
Binding: Binding a ligand exerts forces locally that leads to an imposed local strain. To model this effect at the allosteric site, we choose four adjacent nodes on one side of the system (shown in purple in Fig. 1), and consider that binding at that site imposes a displacement |δR Al on these nodes, as indicated by purple arrows. (Strictly speaking, this description of binding assumes that the ligands are rigid. However we expect our results to hold true qualitatively as long as the ligands are not significantly softer than the protein itself). Minimizing the elastic energy in the entire system with these constraints then leads to a response |δR(σ) Al r that can be extended (see a formal expression for this response in Supplemental Material Section B and [23]). The corresponding energy cost associated with binding writes: where M is the stiffness matrix of the network (whose definition is recalled in Supplemental Material Section B) of dimension N d × N d, which depends on the network considered. The same procedure is used to model the binding of another ligand at the active site (indicated in blue in Fig. 1), allowing us to define a binding energy E Ac (σ). If the two binding events take place simultaneously, the same procedure leads to the derivation of a joint binding energy E Ac,Al (σ). Cooperativity: We seek to engineer materials in which binding at the allosteric site lowers the binding energy at the active site as much as possible, as illustrated in Fig. 2. In the absence of the ligand at the allosteric site, the binding energy at the active site is simply E Ac (σ), whereas if present it is E Ac,Al (σ) − E Al (σ). We seek to maximize the cooperative energy, simply defined as the difference between these terms: which also defines our fitness function. Cooperativity turns out to differ greatly from the geometric task in which a displacement imposed at one end of the material must elicit a given displacement at the other end [23][24][25] (see below and Supplemental Material for a detailed comparison). The architectures associated with the latter task are very asymmetric, in particular they are much softer near the active site than near the allosteric site [23]. By contrast, it is clear from our definition of cooperativity that both active and allosteric sites play a symmetric role. At an intuitive level, the difference can be understood by considering the limit of weak elastic coupling between allosteric and active sites for which one finds E coop ≈ F Ac |δR Al→Ac where |F Ac is the external force field generated by the substrate when it binds to the active site, and |δR Al→Ac is the displacement field induced at the active site by binding a ligand at the allosteric site. Maximizing cooperativity thus requires to have a large and specific response |δR Al→Ac (which is essentially what the geometric task accomplishes) and to have a large force scale |F Ac , which requires the material to be stiff near the active site. This additional constraint makes the cooperative task harder than the geometric one.
Evolutionary Dynamics: To generate cooperative architectures, we implement an evolution scheme which selects preferably networks with high fitness. Specifically, we use a Monte-Carlo algorithm where the relocation of individual springs is considered, i.e. |σ → |σ where a randomly chosen vacant link γ becomes occupied σ γ = 0 → σ γ = 1 and a randomly occupied link α becomes empty, σ α = 1 → σ α = 0. The new structure is selected with the probability p = min[1, exp( F (σ )−F (σ) Te )], where 1/T e is inverse evolutionary "temperature" characterizing the selection pressure.
We find that as the selection pressures increases and T e decreases, there is a rather sudden transition from nonworking networks with zero fitness to cooperative ones, as illustrated in inset of Fig.3. The fitness then appears to plateau, and in what follows we choose T e = 10 −4 where this plateau is reached. Interestingly, in this plateau region we find that the fitness landscape is glassy: there are many families of solutions that are not dynamically connected on the time scale of our runs, implying the presence of large fitness barriers. The families obtained in a given run are defined by the respective initial conditions, and do not display exactly the same fitness as shown in Fig.3. We checked that sequences are much more similar within a family than between different families. Indeed, in a single family the mean overlap between distinct configurations i and j, q ≡ σ i α σ j α α − σ 2 is high with q ≈ 0.36, while it is small q ≈ 0.03 for different families (• α averages over links, and • averages over configurations). Glassiness also implies that the architectures slowly evolve in time, but less and less so as time goes on. In what follows, we study architectures only in the last third of the run, when transient effects are weaker and fitness is nearly stationary. In total, we generated 25 families in d = 2 and 10 families in d = 3.

Analysis Toolbox
In this section we review useful observables characterizing allosteric architectures. Most of them are known in the protein literature, others are novel to the best of our knowledge.
Geometry of allosteric response: By computing the structure of proteins crystallized with and without the ligand bound on their allosteric site, one gets access to the internal response of the protein induced by binding, |δR Al r in our notations. As recently emphasized in this context [10], a key aspect of this response is its strain, which must be zero in parts of the proteins moving as rigid blocks. The strain thus captures where deformation is actually taking place. The strain tensor ↔ (i) can be directly computed from any displacement |δR = {δR i } where i labels particles or nodes, as shown in Supplemental Material Section C or Ref. [33]. Removing the trace leads to a local shear tensor where 1 is a d × d identity matrix. It's useful to define scalar observables to visualize the strain, in particular the shear intensity E shear (i) (not sensitive to compression or dilation) and the bulk intensity E bulk (i) (sensitive to it) as [10]: Rigidity of the structure: For elastic networks, as understood by Maxwell an important aspect of rigidity is the coordination number z(i), counting the local connectivity (number of springs) attached to a node i. This notion, sufficient in our model, can be extended to interactions relevant in proteins as discussed in [34].
Another commonly used observable is the B-factor or Debye-Waller factor [35]. It characterizes the mean square thermal fluctuations of the particle positions. In a harmonic approximation it can be expressed in terms of the vibrational modes (neglecting a temperaturedependent pre-factor): where the ωs and δR ω are frequencies and the corresponding vibrational modes, that are obtained from the diagonalisation of the stiffness matrix. B-factors however may not pick up the interesting flexibility of the structure. For example if a hinge connects two rigid parts, B-factors may be large in the rigid parts too as it is sensitive to rigid motions as well. Here we introduce an observable that would reveal the presence of a hinge, as it characterizes the thermal fluctuations of the strain (which is therefore zero by construction for rigid body). We call it the strain B-factor, which for harmonic dynamics follows: where E shear,ω and E bulk,ω are the shear and bulk intensities for a given mode δR ω (i), as defined from Eqs. (3). Spectral analysis: The response to binding can be decomposed into the vibrational modes [16], which form a complete orthogonal basis. We define the overlap: that satisfies ω q ω = 1.
The extendedness of the vibrational modes is characterized by the participation ratio, defined as: for normalized modes i δR ω (i) 2 = 1. Translations have a unity participation ratio. By contrast, if a mode only involves the motion of ∼ N 0 particles, then P ω ∼ N 0 /N . Conservation: We quantify the local conservation of the structure by considering the mean occupancy, defined over a period of observation τ : If there is no selection pressure on that link, we expect σ. We thus define the conservation Σ to quantify the deviation from this average [23]:

RESULTS
We now document examples of architectures generated by our scheme, focusing on shear, hinge and twist designs. We consider individual families: when average quantities are presented, they always correspond to a time average over the last third of our Monte-Carlo algorithm, as previously described. We then emphasize the features common to all these designs, to be explained in the next section.
Shear design: We start by the two-dimensional case where visualization is easier. If periodic boundary conditions are considered on the horizontal axis (cylindrical geometry), we find that all 25 architectures correspond to a shear design. This is illustrated in Fig. 4A showing the response to binding: except for a linear path connecting the allosteric and active sites, the motion is essentially that of a rigid body (pure rotations and translations). This is most obvious when plotting the map of the shear intensity E shear in Fig. 4B, which is essentially zero excepted along that path. Overall, the design is similar to that of the mint box illustrated in Fig. 4C, where strain also localizes on a hyperplane (a line for d = 2 and a plane fr d = 3). At the structural level, we find that the strain path corresponds to a softer region with lower coordination as shown in Fig. 4D and a larger strain Bfactor, as illustrated in Fig. 4E.
Hinge design. When open boundaries (instead of periodic ones) are used, we find that about 40 to 50 percent of the families lead to hinge architectures, and the rest display a shear design. In the former case, the response exemplified in Fig. 5A can be decomposed into the motion of two rigid bodies connected by a hinge. Again this is most apparent in the map of the shear intensity in Fig. 5B, showing that there is little strain excepted for two disconnected regions near the allosteric and active sites. There is thus no connecting path of high strain between these sites. This design is common in our daily life, as illustrated by the clothespin in Fig. 5C. At the structural level, the map of coordination shown in Fig. 5D and that of strain B factor shown in Fig. 5E display a "H" shape with two rather disconnected region being weakly coordinated with a high strain B-factor.
Twist design. In three dimensions we find a rich variety of architectures, whose structure and response are sometimes hard to describe. Here we present the simple case of a twist architecture, as illustrated in Fig. 6A with the Rubik Cube. To visualize this design, we consider the shear intensity in three sections parallel to the x-z plane as illustrated in Fig. 6B. We find that there is little strain except on the central plane connecting the al-losteric (purple) and active (blue) sites shown in Fig. 6C. There is not however a homogeneous shear on that plane: instead, the strain is low at its center and larger near the boundaries. Further evidence for the twist design appears in the allosteric response itself shown in Fig. 6D with the same slicing geometry: the two side planes show reverse rotating motions, whereas the middle plane shows a more complex displacement pattern. Once again, the structural analysis confirms this view: we find that the coordination is large and the strain B-factor is small except near the boundaries of the central plane, as shown in Figs. 6(E-H). The middle of the central plane thus acts as a well-connected joint around which two quite rigid bodies can rotate.
Universal features of cooperative designs: Our in-silico evolution scheme generates different designs, as illustrated with the examples above. However, all these designed architectures follow the same principles, which we list in the following. These principles are systematically tested by averaging on the 25 families found in two-dimensions with periodic boundaries in Fig.7. The same analysis holds for other boundary conditions and in three dimensions as well, as documented in Supplemental Material Section E: • The system separates into a rigid and a soft manifold, as observed in a class of proteins [36] and in protein models [20].
• The strain associated with the allosteric response is small in the rigid manifold (indicating rigid body or long-wavelength motion), while it is large in the soft manifold. Both properties are apparent in Fig. 7A, showing the two-dimensional density of nodes found with a given strain B-factor (reflecting the local rigidity) and strain intensity (reflecting the strain induced by the allosteric response). This histogram displays a branch of soft nodes, where the strain B-factor is large and positively correlated to the strain intensity.
• In all cases, the mutation cost is high precisely in these locations where the system is soft and where the strain intensity is large, as illustrated in Fig. 7B.
• Most importantly, the daily-life examples we provided all have a common point: they display a single mechanism, i.e. a very soft elastic mode. We observe that this is also true in our cooperative architectures: there is always a single soft and extended mode along which most of the response projects to. This fact is already apparent in the decomposition of the response on vibrational modes shown in Figs. 4H, Fig. 5H and Fig. 6K. It is studied systematically in Fig. 7C and Fig. 7D showing respectively the density of vibrational modes D(ω, P ω ) and the overlap q(ω, P ω ) as a function of both frequency ω and participation ratio P ω . Fig. 7C shows a peak of extended (large P ω ) modes at low ω, Fig. 7D shows that most of the response projects precisely on these modes. We find that essentially one mode governs the response. This result can also be visualized by classifying modes for each system by decreasing overlap q, and by representing the cumulative overlap (the sum of q ω for the r modes with the largest overlap) as a function of the rank r, as illustrated in Fig. 7E. In average, the first mode captures more than 90% of the response.
It is interesting to note that many properties of materials optimized to be cooperative, whose specific property is to display a single soft elastic mode controlling function, differ from materials studied previously optimized to propagate a given strain -below we will refer to both cases as "cooperative" and "geometric" designs. An extensive comparison is performed in Supplemental Mate-rial Section F and G. Salient differences include that: (a) the magnitude of the response is essentially constant in space in cooperative designs (it decays by five fold or more in geometric designs) (b) the cooperative design is symmetric: binding at the allosteric or at the active site leads to a very similar response (whereas elastic information cannot propagate from the active site to the allosteric site in geometric designs) (c) the cooperative design responds much more specifically than the geometric ones (in the latter case, imposing a strain anywhere in the material typically lead to a strong displacement at the active site) and (d) for geometric designs, the response does not correspond to a single soft elastic mode, but to a few of them, as already apparent in Fig. 7E.
To explain the universal features of cooperative designs, and to predict the frequency of the soft extended mode controlling the response, we now investigate the optimality of designs.

THEORY
Absence of design: We now argue that in a continuous elastic medium -where no design is involvedcooperativity decreases very rapidly with the distance L between the allosteric and active sites. Any imposed local strain can be decomposed into multipole moments (dipole and higher), and the slower decaying response in the far field -sufficiently distant from the source -is dipolar, since higher multipoles decay faster. To model the perturbation induced by ligand binding, we may thus consider without loss of generality two dipoles each of magnitude f c, where f is the applied force and c the distance over which these are exerted. Here we give a simple scaling argument for a medium with elastic modulus G. As mentioned earlier, for L c we have E coop ∼ dR|F where |dR is now the dipolar response induced by the first dipole, of magnitude dR(r) ∼ f c r d−1 G , and |F the force field of the second dipole. Since |F is dipolar its scalar product on |dR acts as a derivative taken at r = L, and one obtains e. a very rapid decay with distance. This result is confirmed numerically for the case of a crystalline network in the Section D of Supplemental Material.
Illustration of optimal cooperativity: Shear architecture. We now show that cooperativity can be greatly improved if the material presents a very soft extended mode. For illustration we consider the geometry of Fig.8 where a cylinder of elastic modulus G is cut on its length L, by a band of width c. This generates a zero mode corresponding to the rotation of a square. If displacements at the active or allosteric sites of size δ are imposed as illustrated in Fig.8, they will only couple to that mode (since it costs no energy), and lead to the same response. This statement will be true even if the band of width c is filled up with soft material of elastic modulus mode rank G w , as long as it is small enough (see below). Thus we have E Ac ≈ E Al ≈ E Ac,Al implying E coop ≈ E Al , which can be readily estimated as the amount of elastic energy stored in the soft band, i.e. E coop ∼ L c G w δ 2 . This results implies that E coop = 0 when the material presents a mechanism (i.e. G w = 0), but increases with G w . This argument eventually breaks down, however, when it becomes more favorable to deform the rigid material and to couple to other modes in the system. This takes place when the energy of deforming a continuous active allosteric c δ L FIG. 8: In a cylindrical geometry, a mechanism -or zero mode -can be constructed by slicing the cylinder, as can be achieved by creating a cut of length L and width c. One then obtains an object with the topology of a square, which now displays an additional zero mode corresponding to a rigid rotation. If the cut is filled up with a soft elastic material, the mode gets a finite frequency. As long as it is small, imposing a local displacement as indicated in the figure at the allosteric or at the active site will be dominated by this mode and will lead to essentially the same response. Thus E Ac ≈ E Al ≈ E Ac,Al and Ecoop ≈ E Al . medium of modulus G, E cont ∼ Gδ 2 / ln(L/c) becomes smaller than the energy associated with the soft mode L c G w δ 2 . Comparing these two expressions we get a crossover for G w = G * w with: For G w G * w , the role of the soft mode become negligible and the system will respond as a homogeneous elastic material (whose cooperativity is small as described above). Thus cooperativity will be maximal for G w ≈ G * w , leading to an optimal cooperativity of order: This result is confirmed numerically in the Section D of Supplemental Material. The small energy of the response to binding for large L described by Eq.11 implies the presence of a soft elastic mode, which is relevant experimentally. It can be detected in the vibrational spectrum of the protein, and implies large thermal fluctuations. Such fluctuations, in a harmonic approximation, are inversely proportional to the corresponding eigenvalue of the stiffness matrix, of order λ * ≈ E * coop /||δR|| 2 where ||δR|| 2 is the square norm of the allosteric response. For the shear mode considered ||δR|| 2 ∼ L 2 δ 2 since all particles are moving by a distance of order δ, leading to λ * ∼ 1/(L 2 ln(L/c)). For the vibrational spectrum such a small eigenvalue will lead to a low frequency ω * . Assuming for simplicity that all the particles have identical mass leads to: which is thus much softer that the lowest-frequency plane wave modes, of frequency 1/L [? ]. It is straightforward to extend these results to three dimensions in the geometry of a shear plane, where we find E * coop ∼ Gcδ 2 which does not decay with distance, and ω * ∼ L −3/2 which is now even much smaller than plane waves modes, thus justifying why the spectrum of our materials show an isolated soft extend mode at low frequency. These results are tested in Fig. 9 for d = 2, which confirms that cooperativity is optimal for a finite frequency of the soft extended mode.
Principles of optimal cooperativity: Overall, the common principle emerging from this study is that optimal cooperativity results from the following antagonist effects. On the one hand, the architectures are such that they nearly present an extended mechanism. Because this mode is much softer than others, an imposed strain strongly couples to it, thus allowing to transfer the elastic information over long distances. On the other hand, if this extended soft mode is too soft, the elastic costs associated with binding become too small, leading to a small cooperative energy. As a result, there is an optimal frequency scale for cooperativity.
This idea leads to a natural explanation for the empirical facts listed in the introduction. Indeed shear and hinge designs (i) are clear realizations of this principle, which implies the presence of a soft extended modes at low frequency, consistent with observation (ii).
We expect that our main result, i.e. the existence of an optimal vibrational frequency for cooperativity, will hold true when non-linearities are taken into account. This prediction can be tested using a combination of molecular dynamics (MD) and experiments. MD can be used to measure the effect of point mutations on the thermal fluctuations along the relevant normal mode, and experiments can measure the effect of the same mutation on cooperativity. In the spirit of Fig.9.B, we predict that there is an optimal magnitude of fluctuations for cooperativity to function properly. It would be very interesting to test if proteins function close to this optimum.
Fluctuation-driven cooperativity: Finally, as pointed out in [18,20], the existence of a soft extended modes of frequency ω * leads to the possibility of a cooperative effect with no mean displacement at play, once thermal effects are accounted for (iii). Indeed binding at the active site will hinder motion and increase the soft mode frequency, leading to an entropic cost that can be diminished if binding already took place at the allosteric site. Let us define ω Al , ω Ac and ω Ac,Al the frequencies of the soft mode after binding at the allosteric site, active site and both respectively. We can estimate these quantities as ω 2 Al = ω * 2 + e Al , ω 2 Ac = ω * 2 + e Ac and ω 2 Ac,Al = ω * 2 + e Al + e Ac where e Al (e Ac ) characterizes the additional energy required for the mode to move when a ligand is bound at the allosteric (active) site. Assuming harmonic dynamics, the entropy of a normal mode of frequency ω reads S = k B ln(k B T / ω). Using this expression, one can now estimate the cooperative free energy ∆∆F = −T ∆∆S = k B T ln(ω Ac,Al ω * /ω Al ω Ac ) = −k B T ln(1 − e Ac e Al /(ω * 2 + e Al )(ω * 2 + e Ac )) which can indeed be large if ω * 2 is small compared to both e Al and e Ac .

CONCLUSION AND OUTLOOK
We have used in-silico evolution to design materials which are highly cooperative. Strikingly, the architectures found differ greatly from materials optimized to propagate a geometrical information over long distances. The latter architectures are based on the emergence of a lever that amplifies the mechanical signal where it is desired, which may be relevant in proteins whose task is to trigger large motions when a ligand binds -e.g. to close an ion channel. By contrast, we predict that proteins optimized to be cooperative should display different ar-chitectures, including shear and hinge designs which are well-known in the literature. Intriguingly, we find that there is a great variety of possible functioning architectures, especially in the three dimensional case. However, they all function along the same principle: they nearly display an extended mechanism, whose frequency should be neither too large nor too small for optimal function to occurr.
Our approach rationalizes several empirical observations on allosteric proteins and it also makes testable predictions. In particular, we predict that a single soft extended mode contributes to function, whose frequency should decrease with protein size. We find that this prediction is hard to test stringently from a spectral decomposition of the allosteric response alone, because localized soft modes (typically near the surface of the system) can hybridize with the relevant mode if they lie at similar frequencies. As a result, the response appears to project on a few modes (despite the localized modes being irrelevant for function) instead of one.
Recent methods have been developed in computer science to clean-up spectra of localized modes -see e.g. [37] in the field of community detection. An exciting path forward is to adapt these methods to proteins, allowing one to test if a single extended mode indeed contributes to allostery. Ultimately, this suggests a mechanical approach to discover de novo allosteric proteins, as those in which a single extended mode lies at low frequency in the cleaned-up spectrum. Such an analysis would further predicts where mutations would affect function: we have observed that most damaging mutations hinder the allosteric response, and take place where the extended mode generates high shear.

Reference
[38] appears in the Supporting Material. . This material is based upon work performed using computational resources supported by the "Center for Scientific Computing at UCSB" and NSF Grant CNS-0960316, and by the "High Performance Computing at NYU".

SUPPLEMENTARY MATERIAL
A. Embedding lattice 2D triangular lattice. In our model, we introduce a slight distortion of the lattice to remove long straight lines that occur in a triangular lattice. Such straight lines are singular and lead to unphysical localized floppy modes orthogonal to them. One can remove them by imposing a random displacement on the nodes. Instead, we distort the lines without introducing frozen disorder. We group nodes in lattice by four, labeled as A B C D in Fig. S1. One group forms a cell of our distorted lattice. In each cell, node A stays in place, while nodes B, C, and D move by some distance δ: B along the direction perpendicular to BC, C along the direction perpendicular to CD, and D along the direction perpendicular to DB, as illustrated. We set δ to 0.2, where the straight lines are maximally reduced with this distortion.
3D face-centered cubic (FCC) lattice. We introduce a similar distortion to the FCC lattice. Again, we label the lattice nodes into four different types A B C D, as shown in one layer in the z direction. The nodes are labeled in such a way that all 12 nearest neighbors of a node are different from it. For example, the center node in the bottom panel of Fig. S1, labeled as C, is connecting to two As and two Bs (in solid lines) in the layer and four Ds with two other As and Bs (in dashed lines) out of the layer. To each layer in z direction, there are two other layers, and in those two layers, D and A, B are located at the same x and y. So we only see half of them (two D, one A and one B) connecting to C in a two dimensional projection along the z direction in Fig. S1. We move all As along negative y direction, all Bs along positive x direction, all Cs along (− 1

B. Linear response of elastic networks
Stiffness matrix. Consider a displacement field δ R i ≡ R i − R i0 , where R i0 is the position of the node i in the initial mechanical equilibrium. To the first order in δ R i , the distance among neighboring nodes, defined as r ij ≡ || R i − R j || between node i and j, changes by , wheren ij is the unit vector along link ij from j to i, is the structure matrix.
On the other hand, the force on a node is a composition of tensions, For linear springs on the neighboring connections, f ij = k ij δr ij , the response force to the displacement is, where the stiffness matrix M i,j = lm k lm S lm ,i S lm ,j , depends only on connection |σ and the link directions. The elastic energy corresponds to the displacement field |δR of dimension N d is Linear response to an imposed displacement. When we impose a displacement on the subset E of N E nodes, δR E , forces must be applied on these nodes. All other nodes adapt to a new mechanical equilibrium with no net forces on them, and follow a displacement δR r . Thus Eq.(S3) becomes, for this choice of basis: which leads to: When there are floppy modes in the network, linear equation (S6) may not be solvable. In that case, Q −1 should be understood as the pseudo-inverse so that the network does not respond along the floppy directions (corresponding singular values are zero in Q). Another possibility is to reduce singularity by imposing that each node also interacts with all its next nearest neighbors via weak springs of stiffness k w 1. Both methods lead to qualitatively identical results. For numerical costs, our results were computed using the second approach with k w = 10 −4 . So our stiffness matrix M = S t σ S σ +k w S t w S w . Regarding translations and rotations.
When binding a ligand, the translational and rotational degrees of freedom (TR) of the nodes are not determined. If we write the TR degrees of freedom at the imposed nodes as Ψ E , a dN E × d T R matrix, which is a set of vectors with d T R = 6 in d = 3 and d T R = 3 in d = 2, any imposed displacement giving the same shape change is then, where δR E 0 is purely determined by the shape change, δR E 0 · Ψ E = 0, and c is a parameter vector of dimension d T R to count TR contribution additional to the shape change. We can thus consider a new basis with a dN E by dN E transform matrix U to the original space on imposed nodes so that translations and rotations are isolated from the shape change defined by δ 0 . In this new basis, the forces F imposed on E obey total force and torque balance, The linear response problem thus becomes, Note that given the separation of the two subspaces the matrix U is of dimension dN E × dN E and the matrix I

C. Computing the local strain tensor in a network
In a continuous medium, a motion maps a point X in the reference configuration to a new point x in the current configuration, the strain tensor of the motion can thus be computed as, where a, b labels the spatial dimension. In a discrete medium as networks, the problem is to compute the partial derivative ↔ Λ = ∂ x/∂ X at node i for especially non-lattice structures. Ideally, for any neighbor j close enough in space, where ∆ X ij = R i0 − R j0 and ∆ x ij = R i − R j in our model. We have n b number of such equations for ↔ Λ i when n b neighbors are considered. So ↔ Λ i are usually over-determined when we consider all nearest neighbors (n b = 6 for a 2 × 2 matrix in triangular lattice, and n b = 12 for a 3 × 3 matrix in FCC lattice). Instead of solving Eq.(S15), we define a mean squared error function [33], where we have kept a weight function w j (i) of node j contribution to i in general. Specifically, we set w j (i) = 1 n b for all nearest neighbors to i on the original embedding where x L • are the contributions of the monopole at L. When c/L 1, the contribution outside of both the monopoles and the dipoles dominates the energy, For given displacement δ applied at the dipoles, showing that E coop decays as fast as L −d for two dipoles at a distance L from each other (with weak logarithmic corrections in d = 2). Numerical verification. In Fig. S2(A) we test our prediction for the cooperativity of a homogenous medium without any design (a distorted crystal with δ = 0.2), and confirm Eq. S30 for d = 2. In Fig. S2(B) we test our prediction for an optimal shear design, and confirm the very weal logarithmic decay of the cooperative energy in two dimensions, as described in Eq.11 in the main text. In the main text we have listed the principles underlying the cooperative architectures, and tested them in Fig.7 in two dimensions, with a periodic boundary. The same results hold with open boundaries in d = 2 (Fig.S4) and d = 3 (Fig.S3).

F. Geometric task
Definition: networks perform the geometric task by minimizing a cost function that measures the deviation of the allosteric response |δR Al r from a prescribed shape change located at the active site |δR Ac [23], where |U is a global translation and rotation, which does not change the shape at the active site. Here Ac corresponds to four sites defining the active site. We illustrate the method by studying the case d = 2 with periodic boundaries, as well as d = 3 with free boundaries. Our results are averaged over 25 runs with different initial conditions in d = 2 and 10 runs in d = 3. d = 2: This case is documented in Fig. S5, using the observables introduced in the main text. The architecture presents two main features. Most importantly, the response is non-monotonic and strongly amplified close to the active site, as shown in Fig. S5B. As demonstrated in [23,29], this effect is induced by the presence of a marginally connected region with z = 4 (Fig. S5D) which is thus very soft (Fig. S5E). It can be shown to act as a powerful lever [29]. It is also highly conserved and leads to high mutation costs (Fig. S5G,H). Another aspect of the observed design is the emergence of a shear mode close to the allosteric site, as can be seen from the mean response (Fig. S5A) and from the map of the shear intensity (Fig. S5F). This response is caused by the mergence of a weakly-coordinated band above the allosteric site (Fig. S5D).
Overall, the design is thus similar to that of a fruitpicker (Fig. S5C) where a stimulus (violet arrows) leads to a shear (black arrows) that couples to a head (black arrows), which acts as a lever and leads to a specific desired response. This design leads to a more complex spectral signature where several modes typically contribute to the response (Fig. S5I), instead of one as for the cooperative designs discussed in the main text. d = 3: The arguably most relevant case corresponds to d = 3 with open boundaries, and is illustrated in Fig. S6. As shown in Fig. S6A, we study the response by focusing on two sections: a vertical plane passing through both the allosteric and active sites, and a horizontal plane containing the active site. Once again, we find that the central aspect of the design is the emergence of a weakly-connected region with z ≈ 6 (the isostatic value) surrounding the active site (Fig. S6C), and leading to a very pronounced amplification of the response (Fig. S6B). This lever region is soft (Fig. S6D) and con-served (Fig. S6E,G). In that case, there is no evidence in the shear map of a hinge or shear motion in the material bulk (Fig. S6E). Fig. S6H shows that the lever design alone comes with a rather complex spectral decomposition of the response in which several modes contribute.

G. Comparison between cooperative and geometric designs
The key difference between the geometric and the cooperative designs is that the former develops a lever, while the latter doesn't. This fact leads to vastly different properties of the response to binding.
Amplification of the response: The map of the average magnitude of the response shows a non-monotonic behavior between the allosteric and the active sites (Fig. S7A) for geometric designs, not apparent for cooperative designs (as shown in Fig. S7B,C for both periodic and open boundaries respectively).
Symmetry of the response: For cooperative designs, the response to binding at the active site is very similar to binding at the allosteric site (Fig. S7E,F), because both type of stimuli mostly couple to the single soft elastic mode in the system. For the geometric de- toward imposing a stimulus at the surface of the material, at a site shown by the purple crosses which differs from the allosteric site where these materials were trained to respond. Blue crosses indicate the active site. A large response at the active site is always found in the geometric design (due to the presence of a lever) but not for the cooperative design.