Physical limits to biomechanical sensing
Abstract
Cells actively probe and respond to the stiffness of their surroundings. Since mechanosensory cells in connective tissue are surrounded by a disordered network of biopolymers, their in vivo mechanical environment can be extremely heterogeneous. Here, we investigate how this heterogeneity impacts mechanosensing by modeling the cell as an idealized local stiffness sensor inside a disordered fiber network. For all types of networks we study, including experimentallyimaged collagen and fibrin architectures, we find that measurements applied at different points throughout a given network yield a strikingly broad range of local stiffnesses, spanning roughly two decades. We verify via simulations and scaling arguments that this broad range of local stiffnesses is a generic property of disordered fiber networks, and show that the range can be further increased by tuning specific network features, including the presence of long fibers and the proximity to elastic transitions. These features additionally allow for a highly tunable dependence of stiffness on probe length. Finally, we show that to obtain optimal, reliable estimates of global tissue stiffness, a cell must adjust its size, shape, and position to integrate multiple stiffness measurements over extended regions of space.
Mechanical cues can govern cellular behavior in decisive ways (12). The elastic properties of a cell’s substrate have been shown to guide cell migration (3, 4) and determine cell fate (5, 6). Eukaryotic cells, including fibroblasts, mesenchymal stem cells, and cancer cells, attach to substrates via transmembrane protein complexes called focal adhesions, allowing the cell to sense stiffness (2, 710). Knockdown studies have established that this mechanosensing contributes both to motility and to the regulation of cell shape in threedimensional in vitro systems that closely resemble in vivo cellular environments (1013), where cells are surrounded by a loosely connected, disordered network of protein fibers such as collagen or fibrin (14). These biopolymers form a major structural component of the extracellular matrix (ECM), which serves as the physical scaffolding within which cells live and move. While it is clear that cells actively probe these extracellular networks, it remains unclear how ECM micromechanical properties impact mechanosensing.
Both in vitro experiments (1519) and theory (2024) have demonstrated that biopolymer networks exhibit rich macroscopic mechanical behavior, depending sensitively on network connectivity. However, because the size of a typical cell is comparable to the pore size of the ECM (25), any mechanical information must be inferred by locally probing an extremely heterogeneous material. Although a few studies have begun to characterize this microscopic response (26, 27), a theoretical understanding of how local mechanics are determined by the surrounding heterogeneous structure is still lacking.
How does the intrinsic heterogeneity of the ECM limit a cell’s ability to learn about its global environment from purely local mechanical measurements? For the case of chemical sensing, consideration of the fundamental physical limits dates back to Berg and Purcell’s consideration of noise due to the random arrival of diffusing particles (2830). Here, we take a similar approach to explore the fundamental limits of mechanosensing, where in contrast to chemical signals, the cues are static in time but distributed nonuniformly in space.
To quantify the physical limits of mechanosensing imposed by a cell’s disordered environment, we investigate a simple model consisting of two components: the ECM as an elastic network that deforms in response to external forces, and the cell as an idealized measurement device that probes the stiffness of its surroundings. We found that experimentallyimaged collagen and fibrin networks and randomlygenerated networks all yield a very broad range of modeled local stiffness responses, spanning roughly two decades. We trace the origin of this broad range of stiffnesses to two intrinsic features of disordered networks: (i) the local stiffness depends primarily on a small number of local fibers with consequently large variations, and (ii) these proximal fibers contribute to stiffness in a highly cooperative manner. The combined effect is a large uncertainty in global stiffness inference, which we demonstrate is only further increased by the presence of long fibers, proximity to elastic transitions, and correlations among nearby measurements. Finally, we argue that to obtain accurate estimates of global ECM stiffness, cells integrate multiple stiffness measurements over extended regions of space.
I Results
Cells in connective tissue can glean mechanical information about their surroundings by pulling on the individual biopolymers of the ECM. However, on the short length scale of a typical cell, the measured mechanical response is sensitive to the intrinsic structural disorder of the ECM. To investigate the role of local mechanical disorder in a physiologically relevant system, we considered collagen networks, which form the primary structural component of the ECM (14). We prepared a sample network by reconstituting fluorescentlylabeled collagen typeI monomers and imaged its threedimensional structure (31, Supporting Information). The network is loosely connected (with an average coordination number ) and highly heterogeneous at the cellular scale (with an average mesh size m). This reconstituted collagen architecture was used as an input to construct a mechanical network model where the fibers are treated as elastic beams that can bend and stretch (Fig. 1A). For simplicity, we modeled the stretching and bending of the beams, respectively, as springs and torsional springs connecting pointlike vertices with stretching modulus and bending modulus . Since biopolymers are expected to be much more pliable to bending than to stretching, we chose the stretching modulus such that . The bending modulus of the torsional springs was fitted to the experimental network using data from macroscopic rheology (Supporting Information).
Extracellular networks exhibit broad distributions of local stiffnesses spanning two decades. To study the mechanical response of the network to local forces that might be applied by cells, we defined a “local stiffness” as the linear displacement of two vertices in response to a dipole contractile force along the direction between the vertices (32). We then calculated this local stiffness by numerically solving the equations of force balance for the network. Strikingly, we found that local stiffness measurements yield a very broad range of values, spanning up to roughly two decades (Fig. 1C). Because of the large range of stiffnesses, in what follows we characterize this variability in terms of the geometric standard deviation (Eq. S15), which we find for this collagen network to be .
To explore the generality of this broad local stiffness distribution, we next considered a reconstituted fibrin network (Fig. S1E), as fibrin constitutes the main structural component of blood clots (14). We imaged a fibrin network prepared from a solution of fibrinogen and observed an average coordination number and an average mesh size m (Supporting Information). We found that the fibrin network also has a broad distribution of local stiffnesses, with very similar to that of the collagen network (Fig. S1F). The broad distribution of local stiffnesses for both the collagen and fibrin architectures suggests that this feature may be a general property of biological fiber networks.
To investigate the physical origins of the broad local stiffness distributions described above, we turned to idealized model networks. One way to generate model disordered networks consists of arranging vertices to lie on a regular lattice, such as a simplecubic lattice (SC) or facecenteredcubic lattice (FCC) in threedimensions (Fig. S2A, 22). Connecting these vertices randomly with a fixed probability results in a disordered lattice network. Such lattice networks have provided a useful starting point for characterizing the mechanical response of fiber networks in twodimensions (27). Interestingly, although threedimensional SC and FCC networks also yield a broad range of local stiffnesses (Fig. S4F, Fig. S2B), the geometric standard deviations of these distributions ( and , respectively) are considerably smaller than those of the experimental networks (Fig. 1C and Fig. S1F). A possible reason for this discrepancy is that the disordered lattice networks fail to capture important structural features of real networks at the scale of the cell, such as the random positions of the vertices.
To generate a model network that better matches the structural features
of the experimental networks, we begin by distributing vertices randomly throughout a volume. The density of these vertices is set roughly equal to those of the experimental networks. Pairs of vertices are then connected according to a probability function that depends on the intervertex distance. This probability function was chosen to produce an average coordination number and distribution of fiber lengths that match those of the experimental architectures (Supporting Information). The resulting modeled networks are referred
to as “Random Geometric Graphs” (RGG, Fig. 1B). We computed the local stiffness
distribution for the RGG network and found that the geometric standard deviation is , which is approximately equal to those of the experimental networks (Fig. 1D). Thus, in contrast to the lattice networks, the RGG network appears to quantitatively capture this important local mechanical property of the experimental networks.
The experimental collagen and fibrin networks are in a bendingdominated regime that lies between two macroscopic elastic transitions. What features of the network impact the very large range of local stiffnesses we found for all types of networks? Intuitively, more fiberdense regions should be stiffer. For a fixed density of vertices, the fiber density is roughly proportional to the network connectivity, defined as the average coordination number . We briefly review how the overall connectivity of a network affects its macroscopic mechanical properties (20, 22).
For high values of , the bulk elastic moduli of fiber networks are dominated by the stretching of fibers. As is lowered, the macroscopic response eventually undergoes a crossover from a stretchingdominated to a bendingdominated response, as the network can be deformed without stretching fibers but not without bending them. As is further lowered, another elastic transition occurs when the network ultimately loses macroscopic rigidity. Over the whole range, the macroscopic response depends strongly and nonlinearly on the network connectivity, with the strongest dependence near the two elastic transitions. We expect this sensitivity of the macroscopic stiffness to overall network connectivity to be reflected in a corresponding sensitive dependence of local stiffness on local network connectivity.
To determine how the rapid scaling of the macroscopic response with connectivity manifests locally, we varied the average coordination number of modeled networks for a fixed ratio of the bending modulus to the stretching modulus (with , Fig. 2A, Fig. S2C). At high connectivities, the entire local stiffness response is dominated by stretching interactions and scales with the stretching modulus . As the connectivity is lowered, the entire local stiffness distribution shifts to lower values. Specifically, the median stiffness decreases rapidly and nonlinearly with the average coordination . Near the stretchingbending crossover, the local stiffness distribution becomes bimodal with the emergence of a subset of probes for which the measured stiffness scales with the bending modulus . Below the stretchingbending crossover, the number of stretchdominated probes becomes negligible, and the local stiffness response enters a bendingdominated regime. Within this regime, the median stiffness resumes rapid decay. Finally, as the connectivity is brought below the rigidity transition, an increasing fraction of measurements yields zero stiffness as portions of the network become floppy. We thus find that the elastic transitions of the macroscopic response manifest locally as crossovers between stretchdominated, bendingdominated, and zerostiffness measurements. At these crossovers, the median stiffness decreases most rapidly, and the local stiffness distribution becomes bimodal, yielding a maximum in the geometric standard deviation.
The low connectivity of the experimental collagen and fibrin networks
suggests that they are situated in the bendingdominated regime. At these particular connectivities, the RGG network is poised approximately halfway between the elastic transitions. Here, the median stiffness scales least rapidly, and we therefore expect the proximity to the transitions to play a minimal role in governing the width of the RGG stiffness distribution. To verify this, we systematically varied the ratio of the bending modulus to the stretching modulus for these RGG networks (Fig. S1B,D). We found that the local stiffness measured by each individual probe is nonzero and scales with the bending modulus. This implies that the experimental networks are far from both elastic crossovers, implicating a different source for the broad local stiffness distributions.
The local stiffness depends only on proximal network structure. A very broad distribution of local stiffnesses appears to be a generic feature of disordered fiber networks. What is the origin of this large variance in local stiffness measurements? To investigate how a local stiffness measurement depends on the surrounding network structure, we calculated the stiffness loss upon removing a single fiber (Supporting Information). Intuitively, removal of fibers that are more proximal to the stiffness probe should have a greater effect on , which suggests that on average, the stiffness loss should decay as a function of the distance from the probe center to the center of the removed fiber. Indeed, for all types of networks we studied, we found that the average stiffness loss is consistent with decay (Fig. 1C,D, Fig. S1F, and Fig. S2B, Insets).
This apparent universality of the scaling of the stiffness loss suggests that the powerlaw decay should be calculable within continuum elasticity theory. We therefore calculated the effect of removing a single fiber in the vicinity of a local stiffness probe for the case of a uniform lattice network (Supporting Information). The defect created by removing a fiber at a distance from the probe perturbs the dipole strain field induced by the probe. We can treat this perturbation as an additional dipole strain field originating from the defect, with a magnitude proportional to the initial strain in the removed fiber. Since both this initial strain, and the consequent additional strain “reflected” back to the probe, decay as the strain field of a forcedipole, i.e. as (where is the dimension), the combined effect is an increase in the strain at the location of the probe .
The rapid decay of the stiffness loss due to fiber removal
in three dimensions suggests that the local stiffness is largely determined by the network structure in the immediate vicinity of the probe. Within this small local
region, all quantities are subject to large fluctuations; for example,
since a small region typically contains a small number of fibers,
the variance of fiber density will be large. The universal, rapid decay implies that
the mechanically relevant local structure will always be very
local, with consequently large fluctuations.
Broad stiffness distribution arises from local density fluctuations combined with strongly cooperative contributions of proximal fibers. To quantify the dependence of local stiffness on the surrounding network structure, we considered the number of local fibers , defined as the sum of fibers each weighted by a decaying function of its distance from the probe center (see above and Supporting Information). We find that and are well correlated for all types of networks studied, as shown in Fig. 2B,C. Importantly, has a strong, nonlinear dependence on . Specifically, the center of the marginal distribution of local stiffnesses at fixed increases more rapidly than linearly with . This indicates that local fibers influence local stiffness in a highly cooperative manner, i.e. combining multiple fibers typically results in much larger than additive changes to local stiffness.
To estimate the geometric standard deviation of local stiffness, we must account for these cooperative effects among fibers. We first note that the scaling of the median local stiffness with is consistent with that of the macroscopic shear modulus for all types of networks we studied (Fig. S4DF, S5DF, Insets). This suggests that much of the broad width of the stiffness distribution can be accounted for by the large variations in the local fiber density taken together with the strong, nonlinear dependence of the macroscopic shear modulus on overall fiber density. To test this notion, we estimated the distribution of local stiffnesses by taking transformed by the functional dependence of the macroscopic shear modulus, , where the modulus is that of a macroscopic network with an average number of local fibers given by . Upon accounting for the strong collective effects of fiber density in this manner, we found that the geometric standard deviation of provides a very good estimate for the actual, observed geometric standard deviation of local stiffness (Fig. 1C,D, Fig. S1F, Fig. S2B, Fig. S4DF, and Fig. S5DF). Consequently, the estimator correctly predicts the relative differences in the geometric standard deviation for the different types of networks, including the observation that the stiffness distributions for the experimental and RGG networks are much broader than for the disordered lattice networks.
Our prediction of the width of stiffness distributions is not sensitive to the particular form of the weighting function, provided the decay is rapid enough. For instance, a hard cutoff at twice the mesh size captures a majority of the broad width for the FCC network (Fig. S4E). However, a less rapidly decaying estimator includes a larger number of fibers in the local structure, and consequently the variance in fiber density is smaller. Thus, by comparison, an estimator with weights that decay as yields a smaller, much
poorer estimate of the width of the stiffness distribution (Fig. 1C,D, Fig. S1F, Fig. S2B, Fig. S4DF, and
Fig. S5DF).
Long fibers increase the dynamic range of local stiffness as a function of probe length. Above, we demonstrated how local structure impacts local stiffness for a fixed distance between probed vertices. Biologically, this distance corresponds to the scale over which cells measure stiffness. Since cells can adopt various sizes and shapes, this means that a cell might measure on different scales depending on the length of its “arms,” i.e. the distance between focal adhesions. We therefore consider the local stiffness as a function of the length of the forcedipole probe in our model.
For long probes, the geometric mean and geometric standard deviation of stiffness approach constant values (Fig. 3A,B). This asymptotic saturation occurs because for very long probes, the measured deformation is simply the sum of the deformations of two independent monopole probes, namely the deformation at a vertex in response to a force applied at that vertex. As the probe length is decreased, however, the measured stiffness increases by more than a decade before the probe length reaches the average mesh size, which is a much larger increase than predicted by continuum theory (Supporting Information). Surprisingly, the measured stiffening deviates noticeably from the continuum prediction even for probe lengths that are several times larger than the mesh size.
For the experimental and RGG networks, this anomalously large foldstiffening is present for all values of and larger near the elastic transitions (Fig. S8). In contrast, for the FCC network, the stiffness distribution has almost no dependence on probe length away from the elastic transitions (Fig. S7D). This implies that the dramatic foldstiffening of the experimental and RGG networks largely results from intrinsic features of these networks that are not present for the FCC network. One relevant feature that distinguishes the different types of networks is the distribution of fiber lengths. That is, the experimental and RGG networks both contain long
fibers connecting vertices several times more distant than the mesh
size, whereas for the FCC networks, the typical fiber length (given
by the number of consecutive, coaxial bonds uninterrupted by noncoaxial
crosslinking bonds) is equal to the mesh size, i.e. the lattice spacing. Indeed, a simple extension of continuum elasticity that includes structural heterogeneity on the scale of the fiber length can quantitatively capture the majority of the large foldstiffening observed for the experimental and RGG networks (Supporting Information). Thus, our results indicate that the polydispersity of fiber length in biopolymer networks can increase the foldstiffening as a function of probe length.
Size and shape of sampling region determines accuracy of stiffness inference. The intrinsic heterogeneity of fiber networks presents a challenge for cells attempting to glean information from stiffness cues. That is, tissues with different global stiffness properties may have significant overlap of their local stiffness distributions. In this case, a single local stiffness measurement would provide only a poor estimate of the global stiffness and identity of the tissue. In view of our conclusion that local stiffness distributions are generically broad, how accurately can cells infer the global stiffness of their surroundings from samples of local stiffness?
One possible strategy for cells to increase their inference accuracy is by averaging multiple stiffness measurements. With enough samples of local stiffness, this method would allow cells to reliably distinguish between global environments with different mechanical properties. We can visualize the local stiffness that a cell might infer at each point of the network by plotting the geometric mean stiffness of measurements obtained within a cellsized sphere centered on that point (Fig. 4B, Inset). The patterns in the stiffness landscape reflect the correlations of nearby stiffness measurements and extend over regions larger than typical cell sizes. More precisely, we find that the correlation function between the loglocal stiffnesses measured by two probes decays exponentially as a function of the distance between their centers, with decay lengths around 3 m (Fig. 4B and Fig. S10A, Insets). These spatial correlations of local stiffness arise because nearby stiffness measurements depend on shared local structure. Since correlations reduce the effective number of independent samples, we expect less accurate global inference if stiffness measurements are made closer together in space.
To quantitatively study the reduction in accuracy of cellular mechanosensing due to spatial correlations, we modeled cell inference as an idealized sampling and averaging. To be concrete, we considered the sampled stiffness to be the geometric mean of a random sample of three probes whose centers are contained within prolate spheroids of varying volume and aspect ratio (Fig. 4A, Inset). For all types of networks we studied, the shape of the sampling region directly impacts the uncertainty of stiffness inference (defined as the geometric standard deviation of the inferred stiffness, Fig. 4 and Fig. S10B). Specifically, the uncertainty is always reduced by increasing the volume of the sampling region as well as by increasing its aspect ratio. As the volume and aspect ratio of the sampling region are increased, typical samples of stiffness become increasingly uncorrelated, and, in the asymptotic limit, the uncertainty of the sampled stiffnesses approaches the geometric standard deviation for independent samples.
Ii Discussion
Mechanical cues can guide cell motility and differentiation (36). Indeed, cells are observed to reliably distinguish among twodimensional homogeneous substrates with bulk stiffnesses similar to brain ( kPa), muscle ( kPa), and bone ( kPa) (5). However, the threedimensional in vivo mechanical environment at the cellular scale can be extremely heterogeneous (14, 25). Here, we asked how this intrinsic heterogeneity impacts mechanosensing. Interestingly, we found that within macroscopically homogeneous but locally disordered collagen and fibrin networks, local probes yield a very broad range of stiffnesses, spanning roughly two decades (similar to the relative difference in the bulk stiffness of brain and bone). Moreover, the average measured stiffness is anomalously sensitive to probe length. We quantitatively captured these striking features of the experimental networks with modeled networks, which enabled us to elucidate their physical origins using a combination of simulations and scaling arguments.
We first established that the very broad range of local stiffnesses is universal across network types and spans a wide range of connectivities, including multiple elastic regimes. We then traced the origin of this pervasive broad distribution to variations in local structure. Specifically, stiffness probes are primarily sensitive to a very small region of local fibers, and these privileged fibers contribute to the measured response in a highly cooperative fashion. Finally, we found that the distribution of stiffnesses can be further broadened by tuning specific network features, including proximity to elastic transitions and geometrical disorder. While our results indicate that the collagen and fibrin networks are poised squarely in a regime where the response is dominated by fiber bending, these networks include strong geometrical disorder, including randomlyoriented fiber junctions and a polydisperse fiber length distribution. These structural features explain why the range of stiffnesses is larger for the experimental and random geometric graph (RGG) networks than for latticebased networks. The geometrical disorder of the experimental and RGG networks also introduces an anomalous foldstiffening, i.e. a dependence of stiffness on probe length that exceeds that predicted by continuum elasticity.
Our “idealmechanosensor” model does not address the internal mechanics of cells. Any internal noise in sensing or downstream signaling can only increase overall measurement uncertainty. There is evidence that cells can fix and thus regulate the relative uncertainty in measured stiffness by linearly modulating their applied stress to maintain a constant deformation (8). Notwithstanding, we have shown that even if cells sense nearly optimally, any single stiffness measurement is poorly informative of global tissue properties. This suggests that cells can benefit more from integrating the results of multiple stiffness measurements than from optimizing individual measurements, which may explain why some cells display more than a hundred focal adhesions (33). Yet even this strategy has diminishing returns, because nearby stiffness measurements probe the same underlying local structure and are therefore correlated. To extract useful information, cells must spread their measurements over extended regions of space, either by moving or by extending their shape. The benefits of an elongated shape are twofold, since measurements on larger scales are both more accurate and less spatially correlated. The biological relevance of this strategy is supported by the observation of highly polarized cells over five times longer than wide, including fibroblasts, mesenchymal stem cells, and cancer cells (5, 13, 33, 34).
Since cells and the ECM evolved together, the latter may also be tuned to facilitate mechanical inference. We uncovered two design principles for tuning the mechanical response: proximity to an elastic transition and the presence of long fibers. Both features result in an enhanced foldstiffening for shorter probes, which may be harnessed to provide specialized signals to different cell types that measure over different scales, albeit the foldstiffening comes at the cost of increased variance of the local stiffness. With regard to this variance, while we have only considered isotropic networks, actual tissues can have aligned fibers (12), which may yield a narrower stiffness distribution when probed along a fixed direction. This may explain why the alignment of the ECM has been observed to coordinate with cell polarization to promote migration (35).
Finally, probing stiffness beyond linear response could provide cells with additional information. How hard would a cell need to pull to access the nonlinear regime? Nonlinear effects will certainly become significant when fibers begin to buckle (36). For a single fiber equal in length to the mesh size, the Euler buckling thresholds are and nN using the bending moduli we inferred for the collagen and fibrin networks, respectively. Forces of several newtons are achievable by stronger cells, as well as by colonies that migrate collectively (37). This suggests that the nonlinear regime is an interesting direction for future study.
In summary, the disorder inherent in biological fiber networks places severe physical limits on the accuracy of cellular mechanosensing, suggesting that organisms must have evolved cellularscale strategies to cope with this uncertainty in vivo. Going forward, highthroughput gene deletion and mutation studies in conjunction with realistic patterned substrates (38) can help reveal the full array of internal components and pathways required for accurate mechanosensation.
Iii Networkbased model for mechanosensing
Fibers are modeled as simple elastic elements with a stretching modulus (units of energy/length) and bending modulus (units of energy length). Adding the energy contributions from the stretching and bending of all fibers gives the Hamiltonian for the total mechanical energy of the network to leading order for small deformations:
(1) 
where is the length of the fiber connecting
vertices and in the unperturbed reference state,
is the average unstretched length of fibers and ,
is the difference in vertex displacements , and
is the unit vector in the direction of
the fiber in the undeformed reference state. Here, the first sum corresponds
to stretching interactions and is taken over all fibers . The
second sum corresponds to bending interactions and is taken over all
connected pairs of fibers and .
The cell is modeled as an idealized stiffnessmeasuring device that
exerts a contractile, equal and opposite force of magnitude
on vertices and , which are separated by a distance in
the undeformed network. The vertices are chosen such that they are not directly connected to each other, in order to avoid trivial contributions to the measured stiffness due to the stretching characteristics of a single fiber. Provided the vertices and belong
to a rigid structure, in linear response the applied force deforms the network according
to the forcebalance equations , where is the forceconstant matrix of
the unperturbed Hamiltonian. The measured
local stiffness is defined as the magnitude of
the applied force divided by the relative displacement along the direction
of the force ). We obtain the local stiffnesses numerically by computing the deformations
of the vertices using the generalized inverse
of the forceconstant matrix (Supporting Information).
Acknowledgements.
We thank Yigal Meir for insightful discussions. This work was supported in part by the National Science Foundation Grants DMR1310266 (to L.M.J., S.M., and D.W), the Harvard Materials Research and Engineering Center DMR1420570 (to L.M.J., S.M., and D.W.), a LewisSigler fellowship (to C.P.B.), the German Excellence Initiative via the program “NanoSystems Initiative Münich” (NIM) and the Deutsche Forschungsgemeinschaft (DFG) via project B12 within the SFB1032 (to C.P.B. and F.B.), National Science Foundation Grants PHY1305525 and PHY1066293 (to C.P.B., F.B., and N.S.W.), and the hospitality of the Aspen Center for Physics (to C.P.B. and N.S.W.).A1em \titleformat
[runin] 1em
Supporting Information
Appendix S1 Mechanosensing model
s1.1 Fiber network model.
We study the mechanical properties of crosslinked biopolymer networks, such as those of the extracellular matrix (ECM), using a fiber network model, which consists of a collection of elastic elements that are connected at pointlike vertices. The elastic elements model the stretching and bending interactions of the constituent fibers. The mechanical energy in the fiber network model is given by (29):
(S1) 
where

is the stretching modulus,

is the bending modulus,

is the difference in positions of vertices and ,

is the length of the fiber connecting vertices and in the unperturbed reference state (i.e. the unstretched fiber length),

is the average unstretched length of fibers and ,

is the angle between fibers and that meet at vertex ,

is the angle between unit vectors and that point along the direction of the fibers connecting vertices in the unperturbed reference state.
The first sum corresponds to stretching interactions and is taken over all fibers . The second sum corresponds to bending interactions and is taken over all connected pairs of fibers and .
Small forces will generate deformations that are small compared to the lengths of fibers. In this case, it is convenient to express the positions of the vertices in terms of their deformations about their positions in the undeformed reference state as follows:
(S2) 
Upon expanding both sums in the mechanical energy and discarding higher order terms, we obtain the harmonic energy in the fiber network model, which is given by:
(S3) 
where is the relative deformation of nodes and .
s1.2 Idealized measurement device.
The cell is modeled as an idealized stiffnessmeasuring device that exerts force on vertices of the network. The effect of the applied force perturbs the mechanical energy as follows:
(S4) 
where the index is summed over all vertices. We assume that the ECM behaves as a viscoelastic solid, and that the forces applied by cells change slowly. Assuming forces applied by cells change slower than the timescale set by viscous damping, the network will be deformed in a quasistatic manner, i.e. the network reaches mechanical equilibrium for a given force. In this case, the deformation is completely determined by minimizing the full mechanical energy:
(S5) 
Expressing this equation in terms of the deformations and the forces leads to the equations of forcebalance in static equilibrium:
(S6) 
where is the forceconstant matrix of the unperturbed Hamiltonian:
(S7) 
We model a twopoint measurement performed by a cell as a force dipole, which is defined as a contractile force exerted on vertices and separated by a vector in the undeformed reference state. The vector of equal and opposite forces may be represented using Kronecker delta notation as follows:
(S8) 
where is the magnitude of the force applied to each vertex.
The local stiffness is defined as the linear response of the displacement of vertices and to the contractile force:
(S9) 
s1.3 Numerical procedure.
To calculate the local stiffness defined above, we compute the deformations of the vertices numerically by solving the equations of forcebalance in static equilibrium. These equations do not have a welldefined solution for all possible configurations of the applied force because the forceconstant matrix is singular, i.e. it contains eigenvectors with vanishing eigenvalues, or “zero modes.” If the force of a local stiffness probe couples to a zero mode, the resulting deformations will diverge and the local stiffness is undefined in linear response. Intuitively, this corresponds to probes that act on unconstrained, dangling portions of the network.
To solve the equations of forcebalance while dealing with the technical challenge provided by zero modes, we compute the generalized inverse of the forceconstant matrix. The generalized inverse allows us to first check whether a solution exists for a given force perturbation, and then to solve for the deformations by multiplying the force perturbation by the generalized inverse. Such an approach is efficient because it allows us to measure local stiffness over an entire network using only a single matrix multiplication operation per probe, which is computationally inexpensive compared with solving the system of linear equations anew for each local stiffness probe.
We obtain the generalized inverse by first computing the singular value decomposition (SVD) of the forceconstant matrix:
(S10) 
where and are orthogonal matrices and is a diagonal matrix consisting of singular values equal in number to the rank of the dynamical matrix. From this, the generalized inverse may be computed:
(S11) 
where is a diagonal matrix consisting of the reciprocals of the singular values . The generalized inverse provides a simple way to check whether a solution exists. That is, if a given force vector satisfies the following equation:
(S12) 
then the forcebalance equation has a solution. This guarantees that a given force configuration has no projection onto any zero mode. For welldefined force probes, which satisfy (S12), the resulting deformation of each vertex is finite and given by:
(S13) 
from which the local stiffness may be computed.
Appendix S2 Analysis of network architectures
s2.1 Experimental network architectures.
The network architectures for the fiber network models, i.e. the positions of the vertices and their pairwise connectivities, were obtained from experiment by analyzing images of reconstituted fiber networks. Collagen networks were prepared and imaged as described previously in (36). Briefly, a sample of 0.2 mg/ml protein in 1x DMEM with 25mM HEPES was prepared from a solution of collagen typeI monomers with a small fraction of fluorescentlylabelled monomers at 4 C. Network formation was induced by neutralizing the sample’s pH with 1M NaOH and incubating it at room temperature for four hours. The resulting network was imaged using a confocal microscope (Leica SP5, Wetzlar, Germany). We acquired a set of fluorescent images covering a threedimensional volume that was representative of the network structure. To determine the fiber positions and connectivity, the image stacks were thresholded and subsequently skeletonized, which resulted in a onevoxel thick line representation of all fibers. We define a branchpoint as the junction between three or more fibers. The number of fibers that join at a branchpoint defines the branchpoint’s connectivity . In our mechanical model, the vertices are positioned at the branchpoints and endpoints of fibers. By counting the number of vertices within the network volume, we found a vertex density of vertices/.
Fibrin networks were prepared and imaged as described previously (44). Briefly, a solution of human fibrinogen (Enzyme Research Labs, South Bend, IN) containing a small fraction of fluorescently labelled fibrinogen with a protein concentration of 0.2 mg/ml was prepared in buffer (150mM NaCl, 20mM CaCl, 20mM HEPES, pH 7.4). Network formation was induced by the addition of activated human alphathrombin (Enzyme Research Labs, South Bend, IN) to the fibrinogen solution (final thrombin concentration: 0.1 U/ml). After the sample was allowed to polymerize for 12 hours, the resulting fibrin network was subsequently imaged and the data was processed analogously to the collagen sample. The vertex density was around vertices/.
For the collagen and fibrin networks, we found average coordination numbers of and , respectively. Although these average coordination numbers are low, the networks are both macroscopically rigid, which is defined by a finite response to shear forces. Rheological measurements were performed on networks prepared under experimental conditions analogous to those of the imaged samples. For the collagen network, the measurements were done using an ARG2 rheometer (TA instruments, New Castle DE) equipped with a custommade plasticplate of 25mm and sealed with mineral oil (44). The macroscopic mechanical response was measured by applying a small oscillatory strain and measuring the resulting stress, which resulted in a shear modulus of Pa. The fibrin network was measured with a 40 mm / coneplate geometry. Samples were sealed with mineral oil and allowed to polymerize for at least 12 hours. The network was then perturbed in the same manner as for the collagen network and the shear modulus was found to be Pa.
Macroscopic rigidity requires macroscopic connectedness, i.e. the presence of a spanning cluster of vertices. In our analysis, we identified all vertices that belong to the largest cluster using a unionfind algorithm and removed all other vertices from the network. For both the collagen and the fibrin networks, we considered a spherical sample of radius m. Within this sample volume, the vertices are distributed approximately homogeneously on scales large compared to the mesh size , defined as the radius of a sphere whose volume is equal to the average volume per vertex:
(S14) 
For the experimental collagen and fibrin networks, we found m and m, respectively. On scales comparable to , the networks are intrinsically disordered due to heterogeneity in both the spatial positions of the vertices and the lengths of fibers. The lengths of the fibers are highly polydisperse, spanning over an order of magnitude in length with a substantial fraction of fibers that are very long compared to the average fiber length (Fig. S9A,B). The fiber length distribution peaks at small lengths (i.e. below the average fiber length) and decays roughly monotonically beyond the peak.
For each experimental network, the corresponding mechanical model has two undetermined parameters: the stretching modulus and the bending modulus . Since both parameters may be rescaled together resulting only in a trivial change of units, the only remaining nontrivial parameter is their ratio . Throughout the main text, we chose to capture the fact that individual fibers are much stiffer to stretching than to bending.
We compute the local stiffness distribution numerically using the procedure described above in Section 1.3; however, rather than considering an ensemble of individual probes over many network realizations, we sample many probes throughout the volume of a single network. These procedures are equivalent when the sample volumes are large compared to the scale over which local stiffness measurements are correlated, i.e. networks which are large enough to represent the ensemble via selfaveraging. We confirmed this equivalence by varying the sample volume and observing that the distribution converged below the network radius used in our analysis.
To determine the units of the modeled response, we first noted that both the experimental collagen and fibrin networks display a bendingdominated response. That is, for these networks, the macroscopic shear modulus and all local stiffness measurements are directly proportional to (Fig. S1AD). We therefore fit the bending modulus of our mechanical model using the experimentallymeasured macroscopic response. Specifically, we first calculated the shear modulus numerically using the mechanical model for an arbitrary value of , which we then scaled by the ratio of the experimental shear modulus to the calculated shear modulus. This resulted in fitted bending moduli of nNm and nNm for collagen and fibrin, respectively.
s2.2 Modeled network architectures.
To understand how the mechanical response emerges from the features of the network architecture, we considered idealized model architectures. For simplicity, the resulting structures were generated randomly such that the initial positions of the vertices and fibers were not spatially correlated. Throughout our analysis, we obtained our results for these networks by averaging over a large number of network realizations to reduce sampling noise.
Random geometric graph model Random geometric graphs (RGG) provide a way to generate model architectures with independent control over average coordination number, vertex density, and average fiber length. In particular, these parameters can be chosen to match those observed in experiment. In the RGG model, architectures are generated by first distributing vertices throughout a volume and then connecting them randomly according to a specified probability distribution. The position of each vertex is drawn from a uniform distribution over the network volume, and the probability distribution for connecting two vertices, or the “connection function,” depends only on the intervertex distance. For simplicity, we model the connection function as a power law with an exponential cutoff:
(S15) 
where is a length scale that describes the typical length of fibers. This results in networks that are homogeneous and isotropic over large scales. For such networks in three dimensions, the expected density of vertices that are a distance from a particular vertex scales as , so the expected distribution of fiber lengths in the threedimensional RGG model is .
Throughout our analysis, we considered networks with approximately vertices/ and we chose the length scale equal to m. We chose the average coordination number to match that of the experimental networks (except in Fig. 2A where we varied the average coordination number by randomly removing fibers from a network with a high initial average coordination number). This resulted in an average mesh size m. Due to the larger resulting density of these modeled networks in comparison to the experimental networks, we considered samples of radius m. We determined the bending modulus by fitting the geometric mean stiffness measured for the RGG network at m to that of the collagen network. This resulted in a bending modulus nNm and shear modulus . For the above choices of parameters, the RGG networks quantitatively capture the local mechanical response of the experimental networks, including the broad width of the local stiffness distribution, the large foldstiffening for short probes, and the twopoint stiffness correlation length (see Results). Because of the broad range of stiffnesses, throughout our work, we characterize the variability of the stiffness distribution using the geometric standard deviation, defined as the exponential of the standard deviation of loglocal stiffnesses:
(S16) 
Disordered lattice model Disordered lattice networks are an established class of modeled architectures that have successfully described many of the mechanical properties of biopolymer networks (20, 22, 23). In particular, lattice models are able to capture the macroscopic response and its striking dependence on network density. Because of their numerical convenience, twodimensional latticebased networks have provided a starting point for modeling the microscopic response of fiber networks (27). Here, we provide the first numerical analysis of the microscopic response of threedimensional disordered lattice networks and demonstrate that such models fail to quantitatively capture the broad width we observed for the experimental collagen and fibrin networks.
Disordered lattice networks are generated by placing vertices on an ordered lattice, for example facecenteredcubic (FCC) or simplecubic (SC) lattices in three dimensions, and connecting them randomly with a fixed probability . For a given choice of lattice, the connectivity is the sole parameter that impacts the overall network structure. In particular, the density of vertices cannot be varied independently from the average fiber length. Furthermore, since the vertices are initially positioned to lie on a regular lattice, fibers only meet at discrete, fixed angles in the unperturbed reference state. Disordered lattice networks are therefore significantly simpler than the experimental networks.
The local stiffness distribution for the FCC lattice network is shown in Fig. S2B. Here, we chose a bondprobability of to study the bending regime relevant to biological networks. Throughout our analysis, we set the stretching modulus and bond length . We considered local stiffness measurements for pairs of vertices separated by primitive lattice vectors of the form to avoid pathological effects associated with applying forces in the symmetry directions of the lattice.
For both the SC and FCC disordered lattice networks, we find the presence of a broad range of local stiffnesses spanning up to roughly two decades (Fig. S4F, Fig. S2B). This broad range of stiffnesses is generic over a wide range of network connectivities in the stretching and bending regimes and peaked near the crossovers (see Fig. S2C for FCC lattice results). Interestingly, although the range is qualitatively very large, the geometric standard deviation is substantially smaller for both the SC (, Fig. S4F) and FCC (, Fig. S2B) lattice networks than for the experimental (collagen: , Fig. 1C; fibrin: , Fig. S1F) and random geometric graph (RGG, , Fig. 1D) networks. In the following section, we study how this discrepancy arises from differences in network structure.
What is the physical origin of the qualitatively broad width we observed for the disordered lattice networks? As for the RGG network (Fig. S5DF), the disordered lattice networks also yield macroscopic shear moduli and median stiffnesses that depend sensitively on average coordination number (Fig. S4DF). This suggests that the broad width we observed for all types of networks we studied arises from intrinsic features common to all types of disordered networks. Indeed, as we observed for the RGG network (see Results), we found that for the disordered lattice networks: (1) the stiffness loss also obeys universal scaling (see Fig. S2B, Inset for FCC lattice results), and (2) the scaling of the marginal distribution of local stiffnesses with the number of local fibers is captured by the rapid, nonlinear scaling of the macroscopic shear modulus (Fig. S4DF, Insets).
We then estimated the geometric standard deviation for the SC and FCC networks in the same way as for the experimental and RGG networks (see next section, Dependence of local stiffness on local structure). As before, we estimate the local stiffness distribution as the distribution of the number of local fibers transformed by the macroscopic shear modulus. For both SC and FCC lattice networks, we found that counting the number of local fibers using rapidlydecaying weighting functions results in a good match between the estimated and actual geometric standard deviations (Table S1, Fig. S4DF). Furthermore, we found that the estimate does not depend significantly on the particular form of the weighting function used to calculate the number of local fibers. Specifically, the estimate of the geometric standard deviation for the FCC network does not depend strongly on whether we chose a rapidly decaying powerlaw or a hard cutoff at two mesh sizes (Fig. S4D,E). This suggests that obtaining a good estimate of the geometric standard deviation based on local structure does not require finelytuning the details of the weighting function, provided the chosen weights are given by a function that decays sufficiently rapidly with distance.
Appendix S3 Dependence of local stiffness on local structure
To gain insight into the variation of local stiffness, we determined how the presence of a single fiber in the vicinity of the stiffness probe impacts the local stiffness. We quantified the influence of a fiber by the stiffness loss , defined as the relative change in local stiffness upon removing the fiber from the network (Fig. S3). Intuitively, fibers which are more proximal should have a greater effect. We therefore expect the stiffness loss induced by removing a fiber to decay as a function of the distance from the probe center to the center of the fiber. In this section, we calculate how the stiffness loss scales with distance for an ordered lattice.
s3.1 Stiffness loss scaling.
We calculate the stiffness loss scaling by first deriving an exact formula for the change in local stiffness when a single fiber is removed from an ordered, stretchingonly network in dimensions, and then studying how the magnitude of the effect depends on the distance to the missing fiber from the center of the probe. The equation of forcebalance in static equilibrium (Eq. S6) becomes:
(S17) 
where the unperturbed forceconstant matrix , vertex deformations, and applied forces are represented using Dirac notation as follows:
(S18) 
(S19) 
(S20) 
where if vertices and are nearest neighbors and otherwise, and the vectors form a complete orthonormal set over vertices such that and .
These equations can be inverted by Fourier transformation to solve for the displacements:
(S21) 
where is the lattice Green’s tensor defined by and given by
(S22) 
in terms of the dynamical matrix:
(S23) 
Here, are unit vectors summed over the directions of the nearest neighbor vertices. The lattice Green’s tensor relates the force applied to vertex to the deformation at vertex , which provides a simple way to express the relative deformation of the vertices induced by the dipole:
(S24) 
This relative deformation determines the local stiffness defined in (Eq. S9) as the ratio of the applied force to the relative deformation .
We will now calculate the perturbed stiffness between vertices and when the bond connecting neighboring vertices and is removed, generalizing the analogous calculation for resistor networks performed by Cserti et al. (40). The equation of forcebalance for the perturbed network is
(S25) 
where is the forceconstant matrix for the perturbed network, i.e. the unperturbed forceconstant matrix minus the contributions associated with the missing bond :
(S26) 
(S27) 
where and . We solve for the deformations as before:
(S28) 
in terms of , the lattice Green’s tensor for the perturbed network defined by . By premultiplying this equation with , we obtain a recursive formula for :
(S29) 
We substitute the matrix representation of given above and iterate to obtain an explicit formula for in terms of an infinite series:
(S30) 
where we have defined the constant . We then sum the geometric series to obtain:
(S31) 
This is an explicit formula for the lattice Green’s tensor of the perturbed network in terms of the lattice Green’s tensor of the unperturbed network. The local stiffness measured between vertices and is now
(S32) 
with the perturbation to the measured deformation given by:
(S33) 
Using this exact formula for the local stiffness in the vicinity of a single removed fiber, we now approximate the stiffness loss, defined as the relative change in local stiffness upon removing the bond at a distance (see Fig. S3):
(S34) 
When the location of the missing fiber is far away from the probe, i.e. , the perturbation is small and scales as:
(S35) 
At these long distances, the lattice Green’s tensor approaches the continuum elastic Green’s tensor, which is given by the following (41):
(S36) 
where is the Young’s modulus, is the Poisson’s ratio, is the Kronecker delta function, is a unit vector parallel to the position vector , and is the magnitude of the position vector. Upon substituting the continuum elastic Green’s tensor in place of the lattice Green’s tensor and expanding the perturbation to third order in , we find the stiffness loss scales as
(S37) 
in dimensions. Comparing this universal result to the average stiffness loss measured for disordered networks (Fig. 1C,D, Fig. S1F, and Fig. S2B, Insets), we find that the observed stiffness loss is consistent with the scaling prediction for all types of networks considered.
s3.2 Numerical estimate of the local stiffness distribution.
The rapid decay of stiffness loss suggests that the local stiffness largely depends on a very small number of proximal fibers. To quantify this dependence, we estimate the local stiffness distribution by considering the number of local fibers , given by:
(S38) 
where is a rapidly decaying function that measures the influence of a fiber based on the distance from the center of the probe to the center of the removed fiber:
(S39) 
Here, for simplicity we characterize locality by the strength of the powerlaw decay beyond a shortrange cutoff , which we have introduced to eliminate unphysical divergences. Our introduction of a variable exponent will allow us to later check how the variance of (and consequently our estimate for the local stiffness distribution) depends on the strength of locality. For strong locality, i.e. equal to the rapid decay of stiffness loss in three dimensions, is strongly correlated with for all types of networks studied (Fig. S4AC and Fig. S5AC, Insets). This correlation suggests that explains much of the variance of local stiffness. However, the geometric mean of the marginal distribution of local stiffnesses scales rapidly and nonlinearly with . Consequently, the geometric standard deviation of is less than half of that of . This underestimation occurs because fibers have inherently cooperative contributions to the local stiffness that are not captured by the weighted sum in Eq. (S38), which implicitly assumes that fibers contribute to independently.
To understand how the strength of the collective effects scales with fiber density, we considered the macroscopic shear modulus as a function of average coordination number , or equivalently for a fixed density of vertices. We found that the scaling of the shear modulus with was consistent with the scaling of with (Fig. S4DF, Fig. S5DF). We therefore estimated the local stiffness distribution for all types of networks by transforming the numericallyobtained distribution of by the functional dependence of the corresponding shear modulus for a network with average number of fibers equal to , which resulted in a geometric standard deviation given by:
(S40) 
where
(S41) 
For strong locality, which for concreteness we consider throughout our analysis to be (equal to the rapid scaling of the stiffness loss), the estimator captures the majority of the geometric standard deviation for all types of networks we studied (Table 1, Fig. 1C,D, Fig. S1F, Fig. S2B, and Fig. S4F). Furthermore, the estimate does not depend sensitively on the details of the form of the weighting function, as long as the decay is rapid enough (Fig. S4D,E). Importantly, this estimator captures the large quantitative discrepancy between the width of the experimental and RGG networks versus that of the lattice networks (Table 1).
Collagen  Fibrin  RGG  FCC  SC  

0.54  0.63  0.56  0.37  0.29  
0.49  0.40  0.49  0.30  0.19 
[5 pt] Table 1: Numericallysimulated and theoreticallyestimated geometric standard deviations for all types of networks studied.
However, for weak locality, which we consider throughout our analysis as , the estimated geometric standard deviation is significantly smaller than the observed value (Fig. 1C,D, Fig. S1F, Fig. S2B, and Fig. S4F).
Appendix S4 Continuum theory for the mean stiffness
In this section, we present a simple model for the mean stiffness measured by a local stiffness probe interacting with an isotropic continuum elastic material in three dimensions.
s4.1 Effective continuum elasticmaterial model.
To gain insight into the average mechanical response of a disordered fiber network to a local forcedipole, we approximate the network as an isotropic continuum elastic material (41). The response of such an isotropic continuum elastic material can be completely characterized by two macroscopic moduli that appear in the elastic free energy:
(S42) 
where is the strain tensor for the deformation field , is the first Lame coefficient, and is the shear modulus. To determine these elastic parameters of our model, we calculate the shear modulus and the bulk modulus of the network.
To compute the shear modulus of a network, we simulate the deformation of a cubical region of network in response to a shear stress. We apply a shear stress via force monopoles applied to all nodes within buffer zones at the top and bottom of the network (see Fig. S6A). The deformations of the vertices are obtained by solving the forcebalance equation for static equilibrium using SVD (see Numerical procedure). The shear modulus is defined as
(S43) 
where is the sum of the magnitudes of all the forces applied to the network, is the area of the buffer region, and is the average shear strain obtained from the leastsquares linear fit of the deformation profile in the direction of the applied force.
To compute the bulk modulus, we simulate the deformation of a spherical region of network in response to a uniform compression. We apply a uniform compression via forcemonopoles applied to nodes at the edge of the network (see Fig. S6B). The deformations of the vertices are obtained as before by solving the forcebalance equation for static equilibrium using SVD. The bulk modulus is defined as
(S44) 
where is the sum of the magnitudes of all the forces applied to the network, is the distance from the center of the network to the center of the buffer zone, and is the average radial strain of the nodes within the bulk of the network. Here, the geometric prefactor arises from the spherical geometry.
s4.2 Continuum stiffnessmeasuring forcedipole model.
The local stiffness probe consists of two forcemonopoles separated by a distance which exert a contractile force of magnitude (see Figure S6C):
(S45) 
The deformation produced by each forcemonopole is given by the continuum elastic Green’s tensor (Eq. S36). The continuum elastic Green’s tensor relates applied forces to the resulting deformations of the medium as follows:
(S46) 
Thus, the combination of both forcemonopoles acting on the medium results in a deformation at each point within the medium:
(S47) 
We would like to define a local stiffness probe that measures the sum of the deformation at each monopole along the direction of the probe:
(S48) 
however, in contrast to the network model, here in the continuum model, the forcemonopoles induce diverging deformations at the points where the forces are applied.
To account for the unphysical divergences, we define the measured deformation at each point of force application as the average deformation within a spherical cutoff region of radius centered around the point, where is set equal to the mesh size. Intuitively, this choice of cutoff results in a modeled deformation that corresponds to a coarsegrained description of the deformation of a single vertex in the network. The total measured deformation is then given by:
(S49) 
We evaluate the integral assuming and find the local stiffness measured by the “continuum dipole” is:
(S50) 
where we have defined two constants that depend on the elastic moduli of the network:
(S51) 
(S52) 
The continuum theory results in monotonically decreasing stiffness as a function of probe length, with stiffness asymptotically approaching a constant value for . This occurs because the measured deformations are produced by two forcemonopoles, which add in a linear fashion. Since the selfdeformation of each monopole remains constant, and the contribution of the far monopole at decays to zero for , the local stiffness saturates for large probe lengths, i.e. it approaches a constant value determined by the combined deformations of two noninteracting monopoles.
To study the dynamic range of local stiffness as a function of probe length, we consider the foldstiffening, defined as the ratio of the geometric mean stiffness to the geometric mean stiffness for asymptotically long probes (i.e. ):
(S53) 
This foldstiffening predicted from continuum theory is roughly consistent with the foldstiffening for all types of discrete networks studied when the probe length is at least five times the mesh size (Fig. S7). Despite this match in the foldstiffening at long probe lengths, the observed local stiffness for short probe lengths is always much stiffer than the continuum prediction for the experimental and RGG networks (Fig. S7AC). This discrepancy occurs because the continuum theory does not adequately describe the deformation produced by a monopole at short distances for these networks. This results in a fixed difference between the forcemonopole selfdeformation predicted by continuum theory and the average deformation produced by a forcemonopole applied to the discrete networks.
s4.3 The effect of long fibers.
Below around five times the mesh size, the foldstiffening measured for the experimental and RGG networks begins to deviate rapidly from the continuum prediction (Fig. S7AC). In contrast, the continuum theory seems to capture the foldstiffening of the FCC network, which is fairly constant, even down to probe lengths close to the lattice spacing (Fig. S7D). What determines the length scale of the enhanced foldstiffening we observed for the experimental and RGG networks? In the previous sections, we saw that for fixed probe length, the range of stiffnesses is determined in part by the strong, nonlinear dependence of the marginal distribution of local stiffnesses on network connectivity, with the strongest dependence occurring close to the stretching and bending transitions. Does the scaling of the marginal distribution of local stiffnesses also depend on probe length? If so, the large foldstiffening we observed for the experimental and RGG networks could occur due to the proximity to an elastic transition.
To determine how network connectivity impacts the foldstiffening, we varied the average coordination number for a fixed ratio of the bending modulus to the stretching modulus . We found that for the RGG network, the dependence of foldstiffening on probe length becomes much stronger as the average coordination number is tuned to bring the networks closer to an elastic transition (Fig. S8). However, even far away from the transitions, there is still a large foldstiffening at short probe lengths. Since only the experimental and RGG networks have this large foldstiffening away from the transitions, this property must arise from intrinsic features of the networks that are not present for the FCC network. One feature of the network architecture that explicitly introduces a length scale is the distribution of fiber lengths (Fig. S9).
To predict the probe length for which we expect the crossover to shortrange stiffening to occur due to the presence of polydisperse fibers, we consider the distribution of fiber lengths for each network. For the disordered lattice networks, we define fiber length as the number of bonds in a filament that are consecutive coaxial segments uninterrupted by a noncoaxial connected segment. We approximate the distribution of fiber lengths as proportional to the probability of starting from a fiber endpoint and observing consecutive bonds followed by either a crosslinking (i.e. noncoaxial) bond or the absence of a bond, yielding the following distribution:
(S54) 
where is the coordination number (6 for the SC lattice and 12 for an FCC lattice). Thus, we find an exponentially decaying distribution with a decay length given by:
(S55) 
We see that for the lattice networks, the distribution of fibers decays exponentially. For the FCC network at connectivity , the decay length is around . This is approximately an order of magnitude shorter than the minimum probe length we considered (), which explains why we did not observe the foldstiffening to deviate from the continuum prediction away from the elastic transitions.
In contrast to the relatively monodisperse fiberlength distribution for the FCC network, the experimental and RGG networks contain long fibers that connect vertices separated by distances that are well beyond the mesh size. Beyond the mesh size , the fiberlength distribution decays roughly exponentially with a length scale approximately equal to the mesh size . Thus, we expect the long fibers to modify the foldstiffening out to intermediate probe lengths. Could this effect account for the larger foldstiffening we observed for the experimental and RGG networks?
The longranged structure in the network, e.g. fibers that are several times longer than the probe lengths we considered, precludes us from considering the free energy in the hydrodynamic limit. Thus, to properly describe these networks at these intermediate length scales, we should also consider terms in the free energy that consist, for instance, of gradients of the strain tensor. Although many higherorder terms could contribute at such scales, for simplicity we consider an extension that contains a single additional term (42). Specifically, for each term in the original elastic energy which contains the squaredsum of the diagonal components of the strain tensor or the squaredsum of the individual components , we also add the vector product of the gradient of each term with itself:
(S56) 
where is the partial derivative of with respect to the coordinate represented by index , and is an additional material parameter that corresponds to the extent of the microscopic nonlocality. By setting equal to the decay length obtained from the fiber length distributions of the networks, we show that this gradient elasticity theory captures the essential features shown in our simulations. The Green’s tensor for this nonlocal continuum theory is given by:
(S57) 
where, as in (42), we have defined the two functions
(S58) 
(S59) 
Incorporating the additional effect of the strain gradient into the continuum prediction for the experimental and RGG networks results in a foldstiffening that is approximately a decade at the mesh size , which is significantly larger than the amount predicted by conventional continuum elasticity and captures a majority of the anomalously large foldstiffening observed for the experimental and RGG networks (Fig. S7). In contrast, for the FCC network, the gradient prediction is nearly identical to the conventional continuum prediction, due to the much smaller value of relative to . These results suggest that accounting for the inherent nonlocality from highly polydisperse fibers is crucial to describe the response of biopolymer networks at the cellular scale.
Appendix S5 Generalized stiffness inference
How accurately can cells infer the global stiffness properties of their surroundings based solely on local measurements? As we showed above, the global properties are closely related to the overall behavior of the local stiffness distribution. Specifically, the median stiffness tracks the scaling with network connectivity of the macroscopic response. This suggests that a suitable average of several samples of local stiffness could provide useful information about the bulk properties of the network.
To quantify the accuracy of inference, we considered an idealized sampling process that consists of averaging local stiffness measurements obtained within prolate spheroids. For a given spheroid, the inferred stiffness is given by the geometric mean of measurements randomly chosen among all the probes whose centers are contained within its volume. We define the relative uncertainty of this estimate as the geometric standard deviation of the resulting distribution of inferred stiffnesses. This uncertainty is minimized for independent samples, for which is related to the logarithm of the geometric standard deviation of the distribution of local stiffnesses as follows:
(S60) 
In general, however, stiffness measurements are correlated in space because nearby stiffness measurements depend on shared local structure. To quantify these correlations, we consider the twopoint geometric correlation function for local stiffness measurements and , defined as follows:
(S61) 
We calculated the twopoint geometric correlation function for the experimental and RGG networks and found that were consistent with a decaying exponential form (Fig. 4B, lower Inset, Fig. S10A, Inset):
(S62) 
where is the separation between the centers of the stiffness probes, is the twopoint stiffness correlation length, and is the zerodistance geometric covariance of local stiffnesses whose centers lie at the same point ( is less than because the centers of distinct probes may coincide). The correlation length and the zerodistance geometric covariance are obtained by fitting to the modeled local stiffnesses. Due to the presence of these spatial correlations, the uncertainty of inference will be larger than the uncertainty predicted by assuming samples of local stiffness are independent. Provided that we can approximate the uncertainty using only the twopoint geometric correlation function, the geometric standard deviation (where is the position of sample for each of the samples) is given by:
(S63) 
For our idealized sampling process, the distances between the samples vary because the probes are positioned randomly throughout the network. In this case, each configuration provides an additive contribution to the variance with a weight proportional to its probability. Since the probes are distributed homogeneously and isotropically throughout space, the overall variance of the inference is simply given by integrating the positions of the samples over the sampling volume and normalizing by . Upon substituting the form of from above, we find an uncertainty given by:
(S64) 
For the simple case in which the sampling volume is a sphere of radius , this becomes: