A

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 experimentally-imaged 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.

12

Mechanical cues can govern cellular behavior in decisive ways (1-2). 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, 7-10). Knockdown studies have established that this mechanosensing contributes both to motility and to the regulation of cell shape in three-dimensional in vitro systems that closely resemble in vivo cellular environments (10-13), 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 (15-19) and theory (20-24) 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 (28-30). 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 experimentally-imaged collagen and fibrin networks and randomly-generated 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 fluorescently-labeled collagen type-I monomers and imaged its three-dimensional 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 point-like 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).

Figure 1: Force-dipole stiffness distribution. (A and B) Examples of local stiffness sensing by force dipoles. Modeled deformation under stress from a local force dipole of length m (green arrows) of (A) experimental collagen network and (B) RGG network. Magnitude of fiber deformations indicated by color (small deformations blue, large deformations red). (C and D) Distribution of local stiffnesses , defined as the linear response of local deformation to a force dipole of length m for (C) collagen network and (D) RGG network. Geometric standard deviation of local stiffness indicated by bars (actual distribution black, estimated distribution assuming strong locality red, estimated distribution assuming weak locality orange; see Supporting Information for details). Insets: stiffness loss , defined as the relative change in local stiffness upon perturbing a network by removing a single fiber, versus distance of center of removed fiber from probe center. For collagen, probe length m and removed-fiber length m, and for RGG, probe length m and removed-fiber length m. Error bars are smaller than the size of data points. Dashed lines show asymptotic scaling from continuum theory, which predicts for (see Methods for details).

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.

Figure 2: Collective effect of network structure on local stiffness. (A) Distribution of local stiffnesses for RGG network versus average coordination number of vertices, (force dipole length m, m). (B and C) Joint distribution of local stiffness and number of local fibers , defined as a weighted fraction of local bonds, with bond weight below a short-range cutoff and decaying as beyond, where (B) for RGG network (m) and (C) for collagen network (m). Dashed lines show the macroscopic shear moduli of the networks as a function of average number of local fibers .

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 simple-cubic lattice (SC) or face-centered-cubic lattice (FCC) in three-dimensions (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 two-dimensions (27). Interestingly, although three-dimensional 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 inter-vertex 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 bending-dominated 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 fiber-dense 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 stretching-dominated to a bending-dominated 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 stretching-bending 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 stretching-bending crossover, the number of stretch-dominated probes becomes negligible, and the local stiffness response enters a bending-dominated 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 stretch-dominated, bending-dominated, and zero-stiffness 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 bending-dominated 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 power-law 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 force-dipole, 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. S4D-F, S5D-F, 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. S4D-F, and Fig. S5D-F). 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. S4D-F, and Fig. S5D-F).

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 force-dipole 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 fold-stiffening 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 fold-stiffening 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 non-coaxial 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 fold-stiffening 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 fold-stiffening as a function of probe length.

Figure 3: (A and B) Fold-stiffening , defined as the geometric mean stiffness divided by the geometric mean stiffness for asymptotically long probes (calculated using local stiffnesses of monopole probes), versus probe length for (A) collagen network and (B) RGG network. Error bars show geometric standard deviation of local stiffness; solid curves (blue) from continuum theory (Supporting Information).

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 cell-sized 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 log-local 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.

Figure 4: Minimum uncertainty of stiffness inference. (A and B) Geometric standard deviation of the geometric mean of a random sample of local stiffnesses measured by force dipoles whose centers lie within spheres (A, Inset: blue) and prolate spheroids of equivalent volume (A, Inset: green, aspect ratio 2:1, and red, aspect ratio 3:1) versus radius of spheres for (A) experimental collagen network and (B) experimental fibrin network. Dashed black lines shows the geometric standard deviation of the geometric mean of independent local stiffness measurements. Upper Inset of B shows a two-dimensional slice of the fibrin network stiffness landscape, defined at each point as the geometric mean of local stiffnesses measured by all force dipoles whose centers lie within a sphere of radius m centered on the point. Lower Inset of B shows the two-point geometric correlation for local stiffness, defined as the covariance between log-local stiffnesses and divided by the logarithm of the geometric standard deviation of local stiffness squared , versus distance for collagen (orange) and fibrin (cyan) networks.

Ii Discussion

Mechanical cues can guide cell motility and differentiation (3-6). Indeed, cells are observed to reliably distinguish among two-dimensional homogeneous substrates with bulk stiffnesses similar to brain ( kPa), muscle ( kPa), and bone ( kPa) (5). However, the three-dimensional 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 randomly-oriented 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 lattice-based networks. The geometrical disorder of the experimental and RGG networks also introduces an anomalous fold-stiffening, i.e. a dependence of stiffness on probe length that exceeds that predicted by continuum elasticity.

Our “ideal-mechanosensor” 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 fold-stiffening for shorter probes, which may be harnessed to provide specialized signals to different cell types that measure over different scales, albeit the fold-stiffening 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 cellular-scale strategies to cope with this uncertainty in vivo. Going forward, high-throughput 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 Network-based 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 stiffness-measuring 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 force-balance equations , where is the force-constant 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 force-constant matrix (Supporting Information).

Acknowledgements.
We thank Yigal Meir for insightful discussions. This work was supported in part by the National Science Foundation Grants DMR-1310266 (to L.M.J., S.M., and D.W), the Harvard Materials Research and Engineering Center DMR-1420570 (to L.M.J., S.M., and D.W.), a Lewis-Sigler 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 SFB-1032 (to C.P.B. and F.B.), National Science Foundation Grants PHY-1305525 and PHY-1066293 (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.).
\titleformat

Appendix A

A1em \titleformat

A.1

[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 point-like 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 stiffness-measuring 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 force-balance in static equilibrium:

(S6)

where is the force-constant matrix of the unperturbed Hamiltonian:

(S7)

We model a two-point 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 force-balance in static equilibrium. These equations do not have a well-defined solution for all possible configurations of the applied force because the force-constant 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 force-balance while dealing with the technical challenge provided by zero modes, we compute the generalized inverse of the force-constant 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 force-constant 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 force-balance equation has a solution. This guarantees that a given force configuration has no projection onto any zero mode. For well-defined 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.

Figure S1: Analysis of experimental network architectures. (A and C) Macroscopic shear modulus versus bending modulus for (A) experimental collagen network and (C) experimental fibrin network. (B and D) Distribution of local stiffnesses for a broad range of bending moduli for (B) experimental collagen network (: nNm purple, blue, green, orange, and brown) and (D) experimental fibrin network (: nNm purple, blue, green, orange, and brown). Insets: comparison of stiffness measurements for two different values of bending modulus, and , scaled by bending modulus for (B) collagen network ( nNm) and (D) fibrin network ( nNm). Dashed line shows linear fit. (E) Modeled deformation under stress from a local force dipole of length m (green arrows) of experimental fibrin network. Magnitude of fiber deformations indicated by color (small deformations blue, large deformations red). (F) Distribution of local stiffnesses (force dipole of length m) for experimental fibrin network. Geometric standard deviation of local stiffness indicated by bars (actual distribution black, estimated distribution assuming strong locality red, estimated distribution assuming weak locality orange. Inset: stiffness loss , defined as the relative change in local stiffness upon perturbing a network by removing a single fiber, versus distance of center of removed fiber from probe center (probe length m and removed-fiber length m). Error bars are smaller than the size of data points. Dashed lines show asymptotic scaling from continuum theory, which predicts for .

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 type-I monomers with a small fraction of fluorescently-labelled 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 three-dimensional 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 one-voxel thick line representation of all fibers. We define a branch-point as the junction between three or more fibers. The number of fibers that join at a branch-point defines the branch-point’s connectivity . In our mechanical model, the vertices are positioned at the branch-points and end-points 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 alpha-thrombin (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 AR-G2 rheometer (TA instruments, New Castle DE) equipped with a custom-made plastic-plate 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 / cone-plate 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 union-find 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 non-trivial 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 self-averaging. 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 bending-dominated response. That is, for these networks, the macroscopic shear modulus and all local stiffness measurements are directly proportional to (Fig. S1A-D). We therefore fit the bending modulus of our mechanical model using the experimentally-measured 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 inter-vertex 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 three-dimensional 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 fold-stiffening for short probes, and the two-point 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 log-local 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, two-dimensional lattice-based 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 three-dimensional 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 face-centered-cubic (FCC) or simple-cubic (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 bond-probability 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. S5D-F), the disordered lattice networks also yield macroscopic shear moduli and median stiffnesses that depend sensitively on average coordination number (Fig. S4D-F). 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. S4D-F, 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 rapidly-decaying weighting functions results in a good match between the estimated and actual geometric standard deviations (Table S1, Fig. S4D-F). 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 power-law 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 finely-tuning the details of the weighting function, provided the chosen weights are given by a function that decays sufficiently rapidly with distance.

Figure S2: Analysis of FCC lattice network. (A) Modeled deformation under stress from a local force dipole (green arrows) of disordered face-centered-cubic lattice network (FCC, , bond occupancy , and bending modulus ). Magnitude of fiber deformations indicated by color (small deformations blue, large deformations red). (B) Distribution of local stiffnesses for FCC network. Geometric standard deviation of local stiffness indicated by bars (actual distribution black, estimated distribution assuming strong locality red, estimated distribution assuming weak locality orange). Inset: stiffness loss , defined as the relative change in local stiffness upon perturbing a network by removing a single fiber, versus distance of center of removed fiber from probe center (). Error bars are smaller than the size of data points. Dashed lines show asymptotic scaling from continuum theory, which predicts for . (C) Distribution of local stiffnesses for FCC network versus average coordination number of vertices, , or equivalently, bond occupancy (force dipole length , bending modulus .

Appendix S3 Dependence of local stiffness on local structure

Figure S3: Stiffness loss due to fiber removal. Schematic illustration of the geometry used to calculate the fiber-removal effect. The fiber-removal effect is the relative change in local stiffness when a single fiber is removed at a distance .

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, stretching-only 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 force-balance in static equilibrium (Eq. S6) becomes:

(S17)

where the unperturbed force-constant 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 force-balance for the perturbed network is

(S25)

where is the force-constant matrix for the perturbed network, i.e. the unperturbed force-constant 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.

Figure S4: Dependence of local stiffness on local structure for disordered lattice networks. (A-C) Joint distribution of local stiffness and number of local fibers , defined as a sum of local bonds, for (A) FCC network with bond weight below and decaying as beyond, (B) FCC network with bond weight below and zero beyond, and (C) simple-cubic network with bond weight below and decaying as beyond. (D-F) Distribution of local stiffnesses for (D and E) FCC network (, bond occupancy , and bending modulus where is the stretching modulus), and (F) simple-cubic network (, bond occupancy , and bending modulus ). Dashed lines show estimated distributions (assuming strong locality as in A-C red, weak locality orange, where weak locality is defined by bond weights that decay as beyond the cutoff (or zero beyond the hard cutoff for B). Insets: geometric mean of marginal stiffness distribution versus . Error bars are given by the geometric standard deviations divided by the square root of the number of data points; red curves shows the macroscopic shear modulus versus the macroscopic average of the number of local fibers.

Figure S5: Dependence of local stiffness on local structure for experimental and RGG networks. (A-C) Joint distribution of local stiffness and number of local fibers , defined as a sum of local bonds, for (A) collagen network with bond weight below (m for collagen) and decaying as beyond, (B) fibrin network with bond weight below (m for fibrin) and decaying as beyond, and (C) RGG network with bond weight below (m for RGG) and decaying as beyond. (D-F) Distribution of local stiffnesses for (D) collagen network network (m), (E) fibrin network (m), and (F) RGG network (m). Dashed lines show estimated distributions (assuming strong locality as in A-C red, weak locality orange, where weak locality is defined by bond weights that decay as beyond the cutoff ). Insets: geometric mean of marginal stiffness distribution versus . Error bars are given by the geometric standard deviations divided by the square root of the number of data points; red curves shows the macroscopic shear modulus versus the macroscopic average of the number of local fibers.

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 power-law decay beyond a short-range 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. S4A-C and Fig. S5A-C, 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. S4D-F, Fig. S5D-F). We therefore estimated the local stiffness distribution for all types of networks by transforming the numerically-obtained 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: Numerically-simulated and theoretically-estimated 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

Figure S6: (A and B) Schematic illustration of the procedure used to compute the elastic moduli used in the continuum theory. (A) For the shear modulus, SVD is used to calculate the deformation response to shear forces applied to nodes within buffer zones at the top and bottom of a cubical region of network. The shear modulus is defined as , 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 least-squares linear fit of the deformation profile in the direction of the applied force. (B) For the bulk modulus, SVD is used to simulate the deformation response to compressive forces applied to nodes within a buffer zone at the boundary of a spherical region of network. The bulk modulus is defined as 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. (C) Schematic illustration of the continuum model for mean stiffness of a network. A local stiffness probe is modeled by two force-monopoles, separated by a distance , that apply a contractile force to an isotropic continuum elastic material with moduli calculated from the network. The local stiffness in the continuum model is defined as the sum of the average deformation within spheres of radius centered at the points of force application, where is a short-range cutoff equal to the mesh size, defined for the experimental and RGG networks as the radius of a sphere containing the average volume per node, and for the disordered lattice networks as the lattice spacing.

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 elastic-material model.

To gain insight into the average mechanical response of a disordered fiber network to a local force-dipole, 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 force-balance 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 least-squares 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 force-monopoles applied to nodes at the edge of the network (see Fig. S6B). The deformations of the vertices are obtained as before by solving the force-balance 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 stiffness-measuring force-dipole model.

The local stiffness probe consists of two force-monopoles separated by a distance which exert a contractile force of magnitude (see Figure S6C):

(S45)

The deformation produced by each force-monopole 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 force-monopoles 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 force-monopoles 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 coarse-grained 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 force-monopoles, which add in a linear fashion. Since the self-deformation 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 non-interacting monopoles.

To study the dynamic range of local stiffness as a function of probe length, we consider the fold-stiffening, defined as the ratio of the geometric mean stiffness to the geometric mean stiffness for asymptotically long probes (i.e. ):

(S53)

This fold-stiffening predicted from continuum theory is roughly consistent with the fold-stiffening 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 fold-stiffening 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. S7A-C). 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 force-monopole self-deformation predicted by continuum theory and the average deformation produced by a force-monopole applied to the discrete networks.

s4.3 The effect of long fibers.

Figure S7: (A-D) Local stiffness versus probe length for (A) experimental collagen network, (B) experimental fibrin network, (C) RGG network, and (D) FCC network, showing geometric mean and geometric standard deviation of local stiffness (red); solid curves from conventional continuum elastic theory (blue) and gradient elastic theory (green).

Below around five times the mesh size, the fold-stiffening measured for the experimental and RGG networks begins to deviate rapidly from the continuum prediction (Fig. S7A-C). In contrast, the continuum theory seems to capture the fold-stiffening 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 fold-stiffening 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 fold-stiffening 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 fold-stiffening, 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 fold-stiffening 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 fold-stiffening at short probe lengths. Since only the experimental and RGG networks have this large fold-stiffening 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).

Figure S8: Fold-stiffening versus probe length for RGG network over a range of average coordination number .

To predict the probe length for which we expect the crossover to short-range 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 non-coaxial 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 cross-linking (i.e. non-coaxial) 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 fold-stiffening to deviate from the continuum prediction away from the elastic transitions.

Figure S9: (A-D) Distribution of fiber lengths for (A) experimental fibrin network, (B) experimental collagen network, (C) RGG network, and (D) FCC network.

In contrast to the relatively monodisperse fiber-length 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 fiber-length distribution decays roughly exponentially with a length scale approximately equal to the mesh size . Thus, we expect the long fibers to modify the fold-stiffening out to intermediate probe lengths. Could this effect account for the larger fold-stiffening we observed for the experimental and RGG networks?

The long-ranged 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 higher-order 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 squared-sum of the diagonal components of the strain tensor or the squared-sum 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 non-locality. 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 non-local 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 fold-stiffening 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 fold-stiffening 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 non-locality 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 two-point geometric correlation function for local stiffness measurements and , defined as follows:

(S61)

We calculated the two-point 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 two-point stiffness correlation length, and is the zero-distance 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 zero-distance 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 two-point 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: