Designing allosteryinspired response in mechanical networks
Abstract
Recent advances in designing metamaterials have demonstrated that global mechanical properties of disordered spring networks can be tuned by selectively modifying only a small subset of bonds. Here, using a computationallyefficient approach, we extend this idea in order to tune more general properties of networks. With nearly complete success, we are able to produce a strain between any pair of target nodes in a network in response to an applied source strain on any other pair of nodes by removing only 1% of the bonds. We are also able to control multiple pairs of target nodes, each with a different individual response, from a single source, and to tune multiple independent source/target responses simultaneously into a network. We have fabricated physical networks in macroscopic two and threedimensional systems that exhibit these responses. This targeted behavior is reminiscent of the longrange coupled conformational changes that often occur during allostery in proteins. The ease with which we create these responses may give insight into why allostery is a common means for the regulation of activity in biological molecules.
Introduction
The ability to tune the response of mechanical networks has significant applications for designing metamaterials with unique properties. For example, the ratio of the shear modulus to the bulk modulus can be tuned by over 16 orders of magnitude by removing only 2% of the bonds in an ideal spring network Goodrich2015a (). Such a pruning procedure allows one to create a network that has a Poisson ratio anywhere between the auxetic limit () and the incompressible limit ( in dimensions). In another example, the average coordination number of a network controls the width of a failure zone under compression or extension Driscoll2016 (). Both these results are specific to tuning the global responses of a material. However, many applications rely on targeting a local response to a local perturbation applied some distance away. For example, allostery in a protein is the process by which a molecule binding locally to one site affects the activity at a second distant site Ribeiro2016 (). Often this process involves the coupling of conformational changes between the two sites Daily2007 (). Here we ask whether disordered networks, which generically do not exhibit this behavior, can be tuned to develop a specific allosteric structural response by pruning bonds.
We introduce a formalism for calculating how each bond contributes to the mechanical response anywhere in the network to an arbitrary applied source strain. This allows us to develop algorithms to control how the strain between an arbitrarily chosen pair of target nodes responds to the strain applied between an arbitrary pair of source nodes. In the simplest case, bonds are removed sequentially until the desired target strain is reached. For almost all of the initial networks studied, only a small fraction of the bonds need to be removed in order to achieve success. As was the case in tuning the bulk and shear moduli, we can achieve the desired response in a number of ways by pruning different bonds. We have extended our approach to manipulate multiple targets simultaneously from a single source, as well as to create independent responses to different locally applied strains in the same network.
We demonstrate the success of this method by reproducing our theoretical networks in macroscopic physical systems constructed in two dimensions by laser cutting a planar sheet and in three dimensions by using 3D printing technology. Thus, we have created a new class of mechanical metamaterials with specific allosteric functions.
Our central result is the ease and precision with which allosteric conformational responses can be created with only minimal changes to the network structure. This finding can be viewed as a first step towards understanding why allosteric behavior is so common in biopolymers Gunasekaran2004 (). It has been emphasized that the ability to control allosteric responses in folded proteins could lead to significant advances in drug design Nussinov2013 (); Guarnera2016 (). While much work has focused on identifying, understanding and controlling preexisting allosteric properties, the question of how to introduce new allosteric functions is relatively unexplored Dokholyan2016 ().
Theoretical Approach
Our networks are generated from random configurations of soft spheres in three dimensions or discs in two dimensions with periodic boundary conditions that have been brought to a local energy minimum using standard jamming algorithms OHern2003 (); Liu2010 (); the spheres overlap and are in mechanical equilibrium. We convert a jammed packing into a spring network by joining the centers of each pair of overlapping particles with an unstretched centralforce spring. We chose this ensemble because it is disordered and provides initial networks with properties – such as elastic moduli – that depend on the coordination of the network in ways that are understood Ellenbroek2009 (); Ellenbroek2015 (); Goodrich2015a (). We can work either with the entire system that is periodically continued in space or with a finite region with free boundaries that is cut from the initial network.
Starting with a network with nodes and bonds in dimensions, our aim is to tune the strain between a pair of target nodes in response to the strain applied between two source nodes. (The two nodes comprising each of the target or source are chosen so that they are not initially connected by a bond; see SI). We create a specific response in our system by tuning the strain ratio to a desired value . At each step, we calculate how would change in response to the removal of each bond in the network and then remove the bond which minimizes the difference between and .
We define the response of the network as a vector of bond extensions of length . We use the orthonormal set of bond basis vectors to access the extensions of particular bonds. The extension of bond is attained by , while the strain on the bond is then where is the bond’s equilibrium length. The equilibrium matrix relates the bond extensions to the vector of node displacements of length through the relation Calladine1978 (). Additionally, we define the vector of bond tensions where the tension on bond is . The equilibrium matrix also relates to the vector of net forces on the nodes by .
The matrix encodes the structure of the network and can therefore be used to define a convenient basis for bond tensions. We find this basis by performing a singular value decomposition of and calculating its right singular vectors Pellegrino1993 (). Our result provides a complete basis of size which is composed of two mutually orthonormal subbases. The first subbasis corresponds to the singular values of that are zero; each of these vectors contains a set of tensions on the bonds that do not result in net forces on the nodes. These are commonly known as the states of selfstress (SSS). The second subbasis corresponds to the positive singular values of ; these vectors contain tensions that correspond to net forces on nodes. We call these the states of compatible stress (SCS). Using this basis, we can calculate analytically the change in bond extensions when a bond is removed (see Methods). This gives us a significant advantage over methods that typically require solving a system of equations at each iteration.
Computational Results
We apply our tuning approach to networks with free boundaries in both two and three dimensions (see Methods). We characterize the connectivity of our networks by the excess coordination number . Here is the average number of bonds per node and is the minimum number of bonds needed for rigidity in a network with free boundary conditions Goodrich2012 (). For each trial, a pair of source nodes was chosen randomly on the network’s surface, along with a pair of target nodes located on the surface at the opposing pole. (Note that we could have chosen anywhere in the network for the location of the source and target.) In two dimensions we chose networks that on average had nodes and bonds before tuning, with . In three dimensions networks had on average nodes, bonds and . Prior to pruning, the average strain ratio of the networks in two dimensions was and in three dimensions was for the systems sizes and values we studied. The response of each network was tuned by sequentially removing bonds until the difference between the actual and desired strain ratios, and respectively, was less than 1%.
To demonstrate the ability of our approach to tune the response, we show results for . Note that () corresponds to a larger (smaller) separation between the target nodes when the source nodes are pulled apart. Figure 1 shows a typical result for a twodimensional network: in Fig 1(A), the strain ratio has been tuned to with just 6 (out of 407) bonds removed; Figure 1(B) shows the same network tuned to with a different set of 6 removed bonds. The red lines in each figure indicate the bonds that were pruned. Animations of the full nonlinear responses of these networks are provided in Videos S1 and S2 of the Supporting Information. We note that some of the removed bonds are the same for both and .
The average strain ratio versus the number of removed bonds is shown in Fig. 2(A). Remarkably few bonds need to be removed in order to achieve strain ratios of . In two dimensions only about 5 bonds out of about 400 were removed on average (1%); similarly, in three dimensions only about 4 bonds out of about 740 were removed on average (0.5%). Fig. 2(B) shows the fraction of networks that cannot be tuned successfully to within 1% of a given strain ratio. The failure rate is less than 2% for strain ratios of up to in two dimensions and less than 1% in three dimensions. Therefore, not only does our algorithm allow for precise control of the response, it also works the vast majority of the time. The failure rate increases significantly for , but here we are considering only the linear response of the network. Extremely large values of necessitate an extremely small input strain at the source and may therefore not be physically relevant.
The failure rate is insensitive to except at very small values. In the small regime the failure rate is higher because very few bonds can be removed without compromising the rigidity of the system. If we increase the bond connectivity to for networks in two dimensions, the failure rate remains very low, but bonds are removed in a thin region connecting the target and source. This narrowing of the “damage” region is reminiscent of the results of Ref. Driscoll2016 (), in which bonds above a threshold stress were broken, or of Ref. Goodrich2015a (), in which bonds that contribute the most to either the bulk or shear modulus were successively pruned.
Figure 2(C) shows the distribution of the number of bonds that must be removed to tune a network to within 1% of a desired strain ratio for , , and . These distributions are broad and the mean shifts upwards as increases. The inset shows that the distributions collapse when normalized by the average number of removed bonds . Note that we do not achieve the desired strain ratio simply by tuning the entire free surface of the network to have large strain ratios; the response between the designated target is large while the response between other pairs of nodes is essentially unaffected by the source strain (see Figure S1 in the SI).
Figure 3 shows we are also able to create systems with more complicated responses. In Fig. 3(A), a network with one pair of source nodes controls three targets, each of which has been tuned to a different strain ratio. Figures 3(BI) and (BII) show a single network with two independent sources and targets whose responses were tuned simultaneously and independently of one another. When a strain is applied to the first pair of source nodes, its target responds strongly while the other target does not respond at all. Likewise, when the strain is applied to the second pair of source nodes, its target responds while the first target does not. These networks have ; the failure rate for creating these more complicated responses is generally higher for lower values of in two dimensions.
Experimental Results
Figure 4(A) shows an image of a twodimensional network created by laser cutting a flat sheet. The network is the same as the one shown in Fig. 1(A). The zoomedin areas show the strain response at the target nodes along with the applied strain at the source nodes. The measured value of is consistent with the theoretical prediction. Figure 4(B) shows an image of a threedimensional network created by 3D printing. In this case, the network was designed to have a strain ratio of . The insets again show the relative strains between the pairs of target and source nodes.
In order to obtain a quantitative analysis of how well the physical realizations agree with the simulated networks, we measure the strain on every bond in the twodimensional example when the distance between the source nodes is varied. A majority of the bonds do not change their length appreciably. We therefore focus only on the distance between nodes that were connected by bonds (labeled ) that were removed as the network was tuned. As one might expect, these are the most sensitive to the applied source strain. We calculate, for those changes in distances, the Pearson correlation coefficient between the experiments and the simulations:
Here () is defined as the fractional change due to the source strain in the distance between nodes initially connected by bond as measured in experiments (computer simulations). The standard deviations of and are and , respectively. We find that when averaged over 4 experimental realizations of different designed networks, . This indicates that the experiments are a very accurate realization of the theoretical models.
In contrast to our theoretical models, our experimental systems have bondbending forces that tend to restore angles between bonds to their preferred values. There is also a possibility of buckling out of the plane of twodimensional networks, along with nonlinear effects that are present in real systems undergoing finite strains. In spite of these differences, experiments show that the responses generated by our linear response theory survive into the nonlinear regime probed by the experiments in a significant fraction of the realizations studied.
Discussion
We have shown that it is strikingly easy to tune allosteric deformation responses into an arbitrary spring network by removing only a small fraction of the bonds. Not only can we tune the strain ratio to large negative or positive values for the same network, but we achieve strain ratios of order with almost 100% success. Our theoretical approach can also be extended to more general responses. We can control multiple pairs of target nodes simultaneously with the same pair of source nodes and we can tune multiple independent source/target responses simultaneously into a network. We have also achieved similarly excellent results for tuning responses in periodicallycontinued systems.
The approach we have described here performs a discrete optimization of the response. We have also tuned the response using a standard numerical optimization technique (e.g., gradient descent), by varying the stiffnesses of all the bonds continuously. This bruteforce method is less efficient but equally successful in producing a desired response, and has the advantage of being able to tune nonlinear behavior. Our approach can also be generalized to other types of bond manipulation such as introducing new bonds.
Our theoretical approach provides a framework for understanding and controlling the response of networks relevant to a wide range of fields. For example, networks with builtin localized, longdistance responses could be a novel way of designing architectural structures based on disordered frameworks that have added functionalities. In addition, our theoretical approach can be generalized to other problems such as origami, where one may wish to tune the fold structure so that the system folds in a specific way in response to locally applied external forces Fuchi2015 (). This problem is similar to ours, except that folds are added instead of bonds being removed. Ref. Fuchi2015 () introduces an optimization technique in which fold rigidities vary continuously. This technique is computationally expensive because the network response must be recalculated with each optimization step. A generalization of our theoretical approach to origami, using language similar to that of Ref. SchenkMarkandGuest2011 (), could lead to a more efficient algorithm.
The network responses we create are reminiscent of the localized, longrangecorrelated deformations which characterize allostery in proteins. In fact, folded proteins have long been modeled as elastic networks Bahar2010 () and the response to localized forces in the resulting networks has been studied Atilgan2010 (). Our results demonstrate the ease with which allosteric conformational changes in networks can be achieved by removing a very small set of bonds. As such, it suggests why allostery is so common in large biological molecules.
Similarly, our finding that networks can be tuned to have a variety of different responses may help elucidate multifunctional behavior Favia2009 () and multiple allosterically interacting sites Yuan2015 () in proteins. It has also been observed that small changes in a protein’s covalent structure can often change its biochemical function Jeffery2016 (). One might ask whether our method could be extended to develop a systematic way to determine which intraprotein interactions to modify to create new allosteric functions. Our success in constructing experimental systems in spite of nonlinear and bondbending effects suggests that results are often robust even outside the simple linear regime. However, proteins are thermal whereas our networks are athermal structures. Statistical fluctuations in the structure of proteins has been shown to play an important role in allosteric functionality Tsai2008 (); Motlagh2014 (). It is thus important to investigate how thermal effects can influence the ability to design a desired response. In addition, protein contact networks generally contain prestressed bonds, as well as bondbending and twisting constraints, while our theoretical networks are constructed in the absence of such effects Edwards2012 (); Thorpe2001 (); Srivastava2012 ().
Further work needs to be done to understand why removing specific bonds achieves the desired response. Our method of identifying the elements of the stress basis associated with individual bonds indicates that these stress states are fundamental to this understanding. The dependence on network size and node connectivity also needs to be understood in greater detail. The limits of our algorithm are not yet known, including the number of targets that can be controlled and the number of independent responses that can be tuned for networks of a given size and coordination. To understand experimental systems ranging from proteins to the macroscopic networks we have fabricated, we must extend the theory to include temperature, dynamics, prestress, bond bending, and nonlinear effects due to finite strains. Our approach provides a starting point for addressing these issues.
Materials and Methods
Computed Networks and Choice of Source and Target Nodes
To create a finite network, we choose a cutoff radius from the center of our box and remove all bonds that cross that surface. This process often creates zero energy modes at the boundary of our network. Since we require rigid networks, we remove nodes associated with these modes. We calculate zero modes by performing a spectral decomposition of the dynamical matrix. For each zero mode calculated this way, we identify the node with the largest displacement amplitude and remove it. We then recalculate the zero modes and repeat this process until no zero modes exist. This method of removing zero modes works in any dimension and does not require an arbitrary threshold for whether a node contributes to a zero mode or not. Our final networks are approximately discshaped in two dimensions or ballshaped in three dimensions with nodes and bonds.
We choose the pair of source nodes to lie on the exposed surface of the networks. The pair of target nodes is chosen to be on the opposing pole of the network surface. When choosing a pair of nodes, we also ensure that they are not connected by a bond. This done to avoid surface bonds whose tensions do not couple the the rest of the network. However, since our formalism relies on applying tensions and measuring the strains of bonds, we introduce a bond of zero stiffness, called a “ghost” bond, between each pair of nodes for convenience (see SI).
Further Details of Theoretical Approach
Our approach tunes the ratio of the target strain to the source strain by removing bonds sequentially, one at a time. First, we define the cost function which measures the difference between the network’s response and the desired response . This is given by
(1) 
where indexes the targets and their corresponding sources (For example, in Figs. 1, Fig. 3(A) and Fig. 3(B), respectively.) . Target/source pairs may be defined for the same network response, or for separate independent responses for the same network with different applied source strains. With each step, we choose to remove the bond which creates the largest decrease in .
To decide which bond to remove, we must calculate how the removal of each bond changes . First we define the vectors of bond extensions and bond tensions in response to the externally applied strain, each of length . In order to access the extensions and tensions on individual bonds, we define the complete orthonormal bond basis where indexes the bonds. The extension on bond can then be found, , along with the bond tension, . The strain of bond is where is the bond’s equilibrium length. The tension and extension are related by a form of Hooke’s law,
(2) 
where the flexibility matrix is defined as . Here we choose the stiffness of bond to be where is the bond’s material modulus with units of energy per unit length.
In addition to the bond tensions and extensions, we can define the vectors of node displacements and net forces on nodes . The equilibrium matrix relates quantities defined on the bonds to those defined on the nodes through the expressions and Calladine1978 (). In general , is a rectangular matrix with rows and columns. The total energy can then be written
(3) 
where the Hessian matrix is a matrix. In the presence of an externally applied set of tensions , the minimum energy configuration satisfies
(4) 
To calculate the change in the displacements if a bond were removed, the naive approach would be to set the stiffness to zero for that bond in the flexibility matrix and to solve this equation. However, performing this matrix inversion to test the removal of each bond can be prohibitively expensive with a computational cost of , so we have developed a more efficient approach. Note that here we calculate the response to applied tensions, not the strains we need to calculate . However, since we are only interested in the ratio of the target strain to the source strain and are working in the linear regime, we do not need to explicitly apply a strain.
We use the equilibrium matrix to define a convenient basis of the bond tensions and extensions. Performing a singular value decomposition of gives access to is right singular vectors Pellegrino1993 (). This yields two mutually orthonormal subbases of vectors that together form a complete basis of size . The first subbasis is comprised of vectors with singular values of that are zero; that is, tensions that do not result in net forces on the nodes. These are commonly known as the states of selfstress (SSS), and we denote them as where indicates the particular basis vector. These vectors can also be interpreted as incompatible extensions, or extensions that do not correspond to valid displacements. The second subbasis is comprised of vectors with positive singular values of ; tensions that correspond to net forces on nodes, or extensions that are compatible with node displacement. We call these vectors the states of compatible stress (SCS) and denote them as where indicates the basis vector. In total there are SCS basis vectors and SSS basis vectors which total to .
Using the SSS and SCS subbases, our goal is to find the change in when the stiffness of a given bond is modified. We start by decomposing the bond tensions and extensions,
(5)  
(6) 
Now suppose we apply some external tension to the bonds, . The part of the external tension that projects onto the SCS basis will be balanced by tensions in the bonds, so that . Additionally, the bond extensions that project onto the incompatible extensions, or SSS basis, should be zero because they are unphysical, . Inserting our decompositions of the tension and extension into (2), we get
(7) 
If we project this equation onto the SCS vector , we get a system of equations,
(8)  
(9) 
where is an square matrix. If we invert this system of equations to solve for the extensions, we get
(10) 
The full extension is then
(11) 
In general, calculating the matrix inverse is computationally intensive since it is a square matrix of size . In order to calculate the change in upon the removal of each individual bond, it would appear necessary to invert this large matrix times with a runtime of . This would not necessarily be more efficient than solving (4). To improve this, we can convert freely between a system where all the bonds have different stiffnesses and one where all stiffnesses are unity (see SI). Suppose that all bonds are the same stiffness except for bond which has stiffness . Now we define the unique SCS basis vector:
(12) 
This SCS is closely related to the unique SSS defined in Ref. Sussman2016 (). We can now rotate the SCS basis so that one of the SCS vectors is , making sure to reorthonormalize the rest of the basis with respect to this unique SCS. The benefit of this rotation is that now only the unique SCS contains a nonzero element for bond . The matrix can then be simplified to
(13) 
where we have defined . The resulting extensions are
(14) 
If bond is removed, goes from to giving us
(15) 
From this equation we can calculate the changes in both and and therefore the change . This result can also be derived by inverting (4) and using the ShermanMorrison formula to calculate the change in the inverse of the Hessian Sherman1950 (). Note that this calculation does not include the zero stiffnesses of the ghost bonds, which cannot be mapped to unity with the rest of the system. A generalization of (15) is needed in order to take this into account (see SI).
The next step is to calculate (15) (or its generalization found in the SI) for the removal of each bond. We choose the bond which minimizes in (1) upon removal. One restriction is that we do not choose bonds which introduce zero modes (see SI). Finally, once a bond is chosen, we recalculate the SCS and SSS subbases with the bond removed (see SI).
A summary of our tuning algorithm contains the following steps:

Transform to a system where all bonds initially have the same stiffnesses and add a ghost bond of zero stiffness for each pair of target and source nodes.

Use the equilibrium matrix to calculate the initial SCS and SSS bases.

Calculate the initial extensions of the source and target bonds in response to the applied tension using (11). Note that initially . Use this result to calculate the initial .

For each bond, use the general form of (15) found in Eq. SX of the SI to calculate the change in if that bond were to be removed.

Remove the bond that minimizes in (1). Recalculate the SCS and SSS subbases with the bond removed.
We repeat (3)  (5) until or the process fails.
There are three potential sources of failure represented in Fig. 2(B): cannot be lowered below by removing any bond, no bonds can be removed without creating zero modes, or the numerical error in exceeds . This third source of failure arises because numerical error is introduced as bonds are removed. In order to ensure that our results are accurate, we compare to the value obtained from the solution of (4) with the given set of pruned bonds removed. If the absolute value of the difference exceeds 1%, we call it a failure. Our results constitute an upper bound on the failure rate, which could potentially be reduced by using more accurate techniques to decrease numerical error or more sophisticated minimization algorithms.
Experimental Networks
We create experimental realizations of the theoreticallydesigned networks in both two and three dimensions. These networks consist of physical struts connecting nodes. Inevitably, there are forces that tend to restore angles between bonds to their preferred values and therefore resist any rotation of the bonds at the nodes. We do not include such bondbending forces in our calculations, but we design our networks to minimize such effects.
To make twodimensional networks, we obtain the positions of the nodes and struts from our design algorithm. Next, we laser cut the shape of the network from a silicone rubber sheet. To reduce outofplane buckling, we use 1.6 mm thick polysiloxane sheets with a Shore value of A90. The ratio of strut length to width within the plane of the network is approximately 10:1. To minimize the forces due to bondbending, the struts are made half as wide at their ends  close to each node  as they are at their centers (see Fig. 4(A)). This difference in width at different points ensures that the struts will deform primarily near the nodes rather than buckling at their centers.
To make threedimensional networks, we determine the positions of nodes and struts from the computer simulations and fabricate the networks using 3D printing technology. The proprietary material is a mixture of rubber (simulating styrene based thermoplastic elastomers) and rigid plastic (simulating acrylonitrile butadiene styrene, ABS) with a Shore value of A85. The dimensions of each strut have a ratio of approximately 1:1:11. As in our twodimensional networks, the struts are made thinner at their ends in order to alleviate bondbending.
Acknowledgements
We thank S. Leibler, T. Tlusty and M. Mitchell for discussions that inspired this work, and the Simons Center for Theoretical Biology at the Institute for Advanced Study for its hospitality to AJL. We also thank D. Hexner, A. Murugan, J. Onuchic and D. Thirumalai for instructive discussions. This research was supported by the US Department of Energy, Office of Basic Energy Sciences, Division of Materials Sciences and Engineering under Awards DEFG0205ER46199 (AJL) and DEFG0203ER46088 (SRN) and by a National Science Foundation Graduate fellowship (JWR). This work was partially supported by grants from the Simons Foundation (305547 and 327939 to AJL) and by the National Institute of Standards and Technology under Award 60NANB15D055 (SRN).
References
 [1] Carl P. Goodrich, Andrea J. Liu, and Sidney R. Nagel. The principle of independent bondlevel response: Tuning by pruning to exploit disorder for global behavior. Physical Review Letters, 114(22):225501, 2015.
 [2] Michelle M. Driscoll, Bryan Ginge Chen, Thomas H. Beuman, Stephan Ulrich, Sidney R. Nagel, and Vincenzo Vitelli. Tunable failure: control of rupture through rigidity. arXiv, 2016.
 [3] Andre A. S. T. Ribeiro and Vanessa Ortiz. A chemical perspective on allostery. Chemical Reviews, 116(11):6488–6502, 2016.
 [4] Michael D. Daily and Jeffrey J. Gray. Local motions in a benchmark of allosteric proteins. Proteins: Structure, Function and Bioinformatics, 67(2):385–399, 2007.
 [5] K. Gunasekaran, Buyong Ma, and Ruth Nussinov. Is allostery an intrinsic property of all dynamic proteins? Proteins: Structure, Function and Bioinformatics, 57(3):433–443, 2004.
 [6] Ruth Nussinov and ChungJung Tsai. Allostery in disease and in drug discovery. Cell, 153(2):293–305, 2013.
 [7] Enrico Guarnera and Igor N. Berezovsky. Allosteric sites: Remote control in regulation of protein activity. Current Opinion in Structural Biology, 37:1–8, 2016.
 [8] Nikolay V. Dokholyan. Controlling allosteric networks in proteins. Chemical Reviews, 116(11):6463–6487, 2016.
 [9] Corey S. O’hern, Leonardo E. Silbert, Andrea J. Liu, and Sidney R. Nagel. Jamming at zero temperature and zero applied stress: The epitome of disorder. Physical Review E, 68(1):011306, 2003.
 [10] Andrea J. Liu and Sidney R. Nagel. The jamming transition and the marginally jammed solid. Annual Review of Condensed Matter Physics, 1(1):347–369, 2010.
 [11] Wouter G. Ellenbroek, Zorana Zeravcic, Wim van Saarloos, and Martin van Hecke. Nonaffine response: Jammed packings vs. spring networks. EPL (Europhysics Letters), 87(3):34004, 2009.
 [12] Wouter G. Ellenbroek, Varda F. Hagh, Avishek Kumar, M. F. Thorpe, and Martin van Hecke. Rigidity loss in disordered systems: Three scenarios. Physical Review Letters, 114(13):135501, 2015.
 [13] C. R. Calladine. Buckminster Fuller’s ”Tensegrity” structures and Clerk Maxwell’s rules for the construction of stiff frames. International Journal of Solids and Structures, 14(2):161–172, 1978.
 [14] S. Pellegrino. Structural computations with the singular value decomposition of the equilibrium matrix. International Journal of Solids and Structures, 30(21):3025–3035, 1993.
 [15] Carl P. Goodrich, Andrea J. Liu, and Sidney R. Nagel. Finitesize scaling at the jamming transition. Physical Review Letters, 109(9):095704, 2012.
 [16] Kazuko Fuchi, Philip R. Buskohl, Giorgio Bazzan, Michael F. Durstock, Gregory W. Reich, Richard A. Vaia, and James J. Joo. Origami actuator design and networking through crease topology optimization. Journal of Mechanical Design, 137(9):091401, 2015.
 [17] Mark Schenk and Simon D. Guest. Origami folding: A structural engineering approach. In Origami 5: Fifth International Meeting of Origami Science Mathematics and Education, pages 291–304. 2016.
 [18] Ivet Bahar, Timothy R. Lezon, LeeWei Yang, and Eran Eyal. Global dynamics of proteins: Bridging between structure and function. Annual Review of Biophysics and Biomolecular Structure, 9(39):23–42, 2010.
 [19] C. Atilgan, Z. N. Gerek, S. B. Ozkan, and A. R. Atilgan. Manipulation of conformational change in proteins by singleresidue perturbations. Biophysical Journal, 99(3):933–943, 2010.
 [20] Angelo D. Favia, Janet M. Thornton, and Irene Nobeli. Protein promiscuity and its implications for biotechnology. Nature biotechnology, 27(2):157–167, 2009.
 [21] Yue Yuan, Ming F. Tam, Virgil Simplaceanu, and Chien Ho. New look at hemoglobin allostery. Chemical Reviews, 115(4):1702–1724, 2015.
 [22] Constance J. Jeffery. Protein species and moonlighting proteins: Very small changes in a protein’s covalent structure can change its biochemical function. Journal of Proteomics, 134:19–24, 2016.
 [23] ChungJung Tsai, Antonio del Sol, and Ruth Nussinov. Allostery: Absence of a Change in Shape Does Not Imply that Allostery Is Not at Play. Journal of Molecular Biology, 378(1):1–11, 2008.
 [24] Hesam N. Motlagh, James O. Wrabl, Jing Li, and Vincent J. Hilser. The ensemble nature of allostery. Nature, 508(7496):331–339, 2014.
 [25] Scott A. Edwards, Johannes Wagner, and Frauke Gräter. Dynamic prestress in a globular protein. PLoS Computational Biology, 8(5):e1002509, 2012.
 [26] M. F. Thorpe, Ming Lei, A. J. Rader, Donald J. Jacobs, and Leslie A. Kuhn. Protein flexibility and dynamics using constraint theory. Journal of molecular graphics & modelling, 19(1):60–69, 2001.
 [27] Amit Srivastava, Roee Ben Halevi, Alexander Veksler, and Rony Granek. Tensorial elastic network model for protein dynamics: Integration of the anisotropic network model with bondbending and twist elasticities. Proteins: Structure, Function and Bioinformatics, 80(12):2692–2700, 2012.
 [28] Daniel M Sussman, Carl P Goodrich, and Andrea J Liu. Spatial structure of states of self stress in jammed systems. Soft Matter, 12:3982–3990, 2016.
 [29] Jack Sherman and Winifred J. Morrison. Adjustment of an inverse matrix corresponding to a change in one element of a given matrix. The Annals of Mathematical Statistics, 21(1):124–127, 1950.
Supporting Information
Ghost bonds
On the surface of our networks there are many nodes with exactly bonds in dimensions. Any bond attached to one of these nodes is uncoupled from the rest of the network  applying a tension to one of these bonds does not communicate any tensions or extensions to the rest of the network to linear order. Likewise, no extensions can be measured on these bonds when a tension is applied elsewhere in the network. Therefore, we avoid choosing pairs of source or target nodes that are connected by uncoupled bonds. This is done by ensuring that neither the pair of source nor target nodes share a bond.
However, all calculations involve the bonds, so in order to apply a tension or measure an extension between two nodes, it is convenient if they share a bond. To apply our approach, we introduce a “ghost” bond of zero stiffness between each pair of source or target nodes. These bonds do not affect our results, but allow us to work without any direct reference to the nodes.
Creating identical bond stiffnesses
In order to calculate Eq. (15) in the main text, it was necessary to work in a system where all bonds had identical stiffness. However, we do not want to be restricted to systems that satisfy this special requirement. The bonds in our experimental systems all have the same material modulus , but their equilibrium lengths differ, resulting in bonds with nonidentical stiffnesses . To handle this, we start with a system in which the bond stiffnesses are all different and map it onto an equivalent system in which all the default bond stiffnesses are identical. This is done by introducing a flexibility matrix (as defined below Eq. 2 in the main text) and scaling the equilibrium matrix so that . (Note: We can only scale out stiffnesses that are nonzero.) The energy can then be written in terms of :
(16) 
where the scaled flexibility matrix is proportional to the identity matrix except for any entries that are zero. This energy is the same as that in Eq. (3) in the main text, so the minimum energy configurations should have a onetoone correspondence. Scaled extension or tension vectors are related to the unscaled versions by
(17)  
(18) 
Thus we have implicitly performed all of our calculations on the scaled system and have converted back when calculating .
Modifying multiple bonds
Here we extend Eq. (15) to allow for multiple bonds that do not have identical stiffness . Suppose that all the bond stiffnesses are identically , except for a small subset of bonds which we call . We say a bond if and only if . This set includes any ghost bonds with zero stiffness, along with bonds that are being tested for removal or modification. As discussed in the main text, we typically include just three bonds in – the source and target ghost bonds of zerostiffness, along with the bond tested for removal.
Our goal now is to rotate our SCS subbasis so that as few basis vectors as possible project onto the bonds in . We will denote this new rotated SCS subbasis . We define a special set of basis vectors such that if and only if for some . Typically, the size of will equal the size of .
In order to calculate our rotated SCS subbasis, , we first find the unique SCS for each bond in , which we denote shown in Eq. (12). These unique SCS vectors then are orthonormalized using a modified GrahamSchmidt algorithm. The result is the set of basis vectors described previously. The remainder of the rotated SCS basis is found by using the modified GrahamSchmidt algorithm to orthonormalize the original SCS basis with respect to the set of vectors , throwing out any vectors that are completely zeroed out. The result is our set of orthonormal rotated SCS vectors. However, it will be shown that only the vectors in will be necessary for our solution.
Each basis vector that is not in has zero projection onto bonds that are in , i.e. if , then for all . This means that if either or , then and are orthogonal over a reduced basis such that
(19) 
This new basis now gives us the means to rewrite for or ,
(20)  
(21)  
(22)  
(23) 
The total matrix is then
(24) 
where we have defined the submatrix for . We see that the matrix inversion problem is simplified to just inverting .
(25) 
Since the size of is very small (in our case typically a set of size ), calculating the inverse of this matrix is very fast. The extension can now be represented as
(26) 
The change in extension on bond when the stiffnesses are modified is then
(27) 
Note that the solution only depends on the basis vectors in . This means that only this small number vectors must calculated and the rest my be neglected.
Avoiding the introduction of zero modes
We impose the constraint that we do not introduce any zero energy modes into the system when we remove bonds. We ensure this by only removing bonds which contribute to the SSS subbasis. By MaxwellCalladine counting, , where is the number of zero modes [13]. This means that if we remove a bond, we can either add a zero mode (increase ) or remove a SSS (decrease ). If a bond is removed that contributes to the SSS subbasis, a unique SSS will also be removed and no zero mode will be created [28]. This unique SSS, which we define as for bond , is calculated analogously to the unique SCS, , shown in Eq. (12). We find that
(28) 
As long as , then bond contributes to the SSS subbasis.
Removing bonds from SSS and SCS subbases
When a bond is removed, we remove the unique SSS vector from our SSS subbasis by subtracting off its projection onto each SSS basis vector. We also remove the entries for bond from all vectors in both our SSS and SCS subbases. The two bases are then reorthonormalized using a modified GrahamSchmidt algorithm. This procedure dominates the computational complexity of our algorithm with a computational cost of . However, this is still significantly faster than solving Eq. (4) each time a bond is removed, which has a cost of .
Animations of nonlinear response
Videos S1 and S2 show animations of the responses of the networks in Fig. 1 in the main text. Although our algorithm only considers and controls the linear response, we show the full nonlinear deformations. The resulting strain ratio in the nonlinear regime is within a factor of two of the desired linear result up to source strains of .
To calculate the nonlinear response, we start with a tuned network and minimize the nonlinear configurational energy for increments of the source strain from 0% to 40%. In the network’s undeformed state, each node has an initial position vector . If the network is deformed in some way, then each node will have a new position
(29) 
where is the node’s displacement vector. If two nodes share a bond, then the bond vector going from node to node is with magnitude , while the deformed bond vector is
(30)  
(31)  
(32) 
where is the bond extension vector. The configuration energy summed over all bonds with centralforce harmonic potentials is then
where and is the stiffness of bond . We minimize this energy numerically with respect to the displacement vectors under the constraint that the source strain is a specified value.