Nonlinear mechanosensation in fiber networks

In a diversity of physiological contexts, eukaryotic cells adhere to an extracellular matrix (ECM), a disordered network with complex nonlinear mechanics. Such cells can perform mechanosensation: using local force probing they can measure and respond to their substrate's mechanical properties. It remains unclear, however, how the mechanical complexity of the ECM at the cellular scale impacts mechanosensation. Here, we investigate the physical limits of mechanosensation imposed by the inherent structural disorder and nonlinear elastic response of the ECM. Using a theoretical framework for disordered fiber networks, we find that the extreme mechanical heterogeneity that cells can locally sense with small probing forces is strongly reduced with increasing force. Specifically, we predict that the accuracy of mechanosensation dramatically improves with force, following a universal power law insensitive to constitutive details, which we quantitatively confirm using microrheology experiments in collagen and fibrin gels. We provide conceptual insights into this behavior by introducing a general model for nonlinear mechanosensation, based on the idea of an emergent nonlinear length-scale associated with fiber buckling. This force-dependent length-scale enhances the range over which local mechanical measurements are performed, thereby averaging the response of a disordered network over an enlarged region. We show with an example how a cell can use this nonlinear mechanosensation to infer the macroscopic mechanical properties of a disordered ECM using local measurements. Together, our results demonstrate that cells can take advantage of the inherent nonlinearity of fibrous networks to robustly sense, control, and respond to their mechanical environment.

Cell behavior is steered by various cues from their extracellular environment.Such cues include chemical, electrical, and topographic signals that regulate key cell functions such as migration [1] and thereby impact processes ranging from embryonic development [2][3][4] and tissue maintenance to disease progression [5,6].In particular, there is growing evidence for mechanosensation: cells sense and respond to the mechanical properties of their environment [5,[7][8][9][10].The stiffness of the cell's substrate can guide developmental processes in vivo, such as axonal growth [12].In vitro model systems further revealed that cells mechanically probe their substrate and subsequently modify behaviors such as differentiation [13], gene expression [14], and motility [11,15,16].It remains unclear, however, what mechanical information cells can perceive inside the complex environments they encounter naturally [17,18].
In vivo, many cell types mechanically interact with the extracellular matrix (ECM) by adhering to network fibers and exerting local forces.The polymerization and gelation processes through which these collageneous matrices form result in an inherently disordered fiber network with large structural variations at the cell scale [19,20].Consequently, the mechanical properties cells can locally perceive depend strongly on network location [21][22][23][24][25][26].
This implies that cells face a highly heterogeneous mechanical environment in which the cell-scale linear stiffness, measured locally at different locations in a single network, exhibits relative variations as large as between the macroscopic stiffness of tissues as distinct as brain (1 kPa) and bone (100 kPa) [22].Thus, even if cellular mechanosensors are ideal and can perfectly measure the local linear mechanical response, mechanosensation remains limited by matrix disorder.Importantly however, cells can exert forces of up to few nanonewtons [27][28][29][30] to probe their environment, easily exceeding the linear response regime of the ECM [31][32][33].Indeed, collagen networks exhibit a pronounced nonlinear response at relatively small stresses or deformations.It is unclear how such nonlinearities impact cellular mechanosensation inside the disordered ECM.
Here, we investigate nonlinear mechanosensation from the perspective of a cell as an ideal mechanosensor inside a disordered matrix.Using a network model, we study the relation between bulk and local network properties.We derive a simple relation describing how the fiber constituent stiffening and density control the macroscopic nonlinear mechanics, which explains the unique density-independence of the nonlinear stiffness of collagen networks [42].Surprisingly however, we find that these macroscopic nonlinear properties are largely irrelevant microscopically.At the scale of a cell, the local mechanical response is strongly modulated by variations in fiber density.Using both theory and experiments, we discover a generic power law decay of the variability of local stiffness measurements with force.For a range of fiber constitutive nonlinearities, the local stiffness becomes increasingly insensitive to network disorder.To provide conceptual insights we develop a model for this nonlinear mechanosensation.Our model shows how large probing forces applied by a mechanosensor induce fiber buckling over an extended range, thereby effectively enhancing the length-scale over which a mechanical measurement is averaged in a disordered network.Thus, we here find that elastic nonlinearities can be exploited by cells to overcome the inherent disorder of their environment such that local measurements can be employed to accurately infer the macroscopic mechanical properties of the ECM.

MODEL FOR DISORDERED NONLINEAR FIBER NETWORKS
To investigate the consequences of structural heterogeneity in the ECM on the local mechanical environment cells can perceive, we build on a broadly used minimal model for a disordered fibrous matrix [36,63].In this model (Materials and Methods), we introduce structural disorder by randomly depleting bonds on a regular lattice.The lattice fibers are represented by these bonds that are present with a probability p, setting the fiber density.Fiber bonds resist both bending and longitudinal deformations.Here, we describe the longitudinal fiber response with a nonlinear force-extension constitutive law (CL) τ = f ( ), with the fiber's deformation and τ its tension.This CL is chosen to be asymmetric in compression and tension.Indeed, fibers buckle and soften under compression ( < 0) and stiffen beyond a characteristic tension, with a power law increase k ∝ τ x of their differential stiffness k = dτ /d , where the exponent x characterizes the stiffening mechanism.Of particular interest is the case x = 1, corresponding to an exponential CL (EXP, Figs.1A-C), reflecting the empirically established stress-strain relationship of tendon and reconstituted collagen networks [42,64].To assess the effects of the nonlinear fiber micromechanics on cellular mechanosensing, we also consider three other CLs, Figs.1A-C.These CLs exhibit buckling-induced fiber softening and several distinct tensile responses described by a power law k ∝ τ x : x = 0 (LIN) corresponds to linear non-stiffening springs, x = 1/2 (QUAD) to a quadratic force-extension relation, and x = 3/2 (WLC) describes a divergent entropic force-extension relation of the worm-like chain model [34,45].Throughout this article, we use the mechanical equilibrium response to global and local loading of these model networks to study the ability of cells to glean information on the stiffness of their surrounding heterogeneous environment.FIG. 2. Stress-and fiber density-dependence of the bulk modulus.Macroscopic differential elastic modulus (K) versus bulk stress (σ) of 3D depleted networks with various fiber densities set by the bond occupation parameter p (low to high value indicated by a dark to bright color) and constituted of A: EXP-fibers, B: LIN-fibers, C: QUAD-fibers and D: WLCfibers, the inset highlights the p−dependence in the large stress regime.Tension distribution in 2D EXP-fiber networks (same color code as Fig. 1C) of E: randomly depleted networks with, from left to right, p = 0.5, 0.6 and 0.8 at a fixed stress in the asymptotic large stress regime.F: p-dependence of the large-stress modulus of 3D depleted networks.The normalization constant K1 is the modulus for the p = 1 homogeneous network.G: same as E for regular networks of decreasing fiber density.H: n-dependence of the large-stress modulus of the regular networks.The normalization constant K1 is the modulus for the n = 1 denser network.I: K versus σ of regular networks with low to high fiber densities for the different CLs, J: corresponding rescaled response, with A0 and σ0 CL-dependent constants.
to some regions of the network being denser than others, thereby modulating the local mechanical properties [21][22][23][24][25][26].To understand how, we first investigate how fiber density affects the bulk stiffness of the network.While this dependency is well-understood in the linear, lowstress regime [63,65], this is not the case in the nonlinear stiffening regime arising at larger stress, which is relevant for cell-ECM mechanical interactions [31][32][33].
We simulate the response to a dilatation strain of networks with varying fiber density of EXP-fibers representing collagen (Fig. 1D, Fig. 2A).The network stiffness is quantified by the differential bulk modulus K = dσ/dγ, with σ and γ the macroscopic stress and strain (Materials and Methods).At low stress, the response is linear: network stiffness is stress-independent and, intuitively, increases with fiber density.In contrast, after a cross-over at intermediate stresses, the modulus increases with stress as a power law K ∝ σ that reflects the fiber constitutive nonlinearity.Strikingly, in this regime the macroscopic elastic responses converge to a stresscontrolled value, insensitive to fiber density, consistent with macrorheology experiments on collagen gels [42].
To elucidate the stress-and density-dependence of the nonlinear macroscopic response of collagen, we propose a differentially affine model.Indeed, we observe that at large stresses a tense subnetwork that carries most of the stress emerges and remains stable under further loading (Fig. 2E, Supplementary Movie S1).In this regime, we assume that the stress is evenly distributed among the bonds of this load-bearing subnetwork.These bonds have a density n and a tension τ , resulting in a macroscopic stress In our model, we further assume that an increase δγ of the macroscopic strain results in an equal stretch δ = δγ of the load-bearing fibers, i.e. that the system is differentially affine.This implies that the macroscopic differential modulus K at large stress directly reflects the microscopic differential stiffness k = dτ /d of the loadbearing fibers: Importantly, Eqs. ( 1) and (2) imply strong constraints connecting the stress-and density-dependence of the bulk modulus.Indeed, for fibers with power-law stiffening, k ∝ τ x , we find at large stress that This equation implies that a single exponent x of the fiber-level CL controls both the n-and σ-dependence of the macroscopic modulus.Strikingly, for collagen-like fibers with x = 1, we find K ∝ σ, independently of n.Our simple differentially affine model thus recapitulates the observations for collagen at large stress.While an alternate explanation involving normal stresses was previously proposed under shear [42], our model proposes a simple and general rationalization of the densityindependence of the nonlinear elastic modulus of collagen.
Our differentially affine model (Eq.( 3)) also makes predictions for other CLs.For linear elements (x = 0), we recover a stress-independent modulus proportional to n.For 0 < x < 1, K(σ) increases with the load-bearing bond density.Remarkably, if x > 1, we predict that, K(σ) decreases with n.This startling behavior can be understood by considering the loading of a set of two bonds in parallel.If one of these segments is cut, the load is transferred to the remaining bond, doubling its load.If x = 1, however, this would also double the bond stiffness, thus leaving the rigidity of the system unchanged.For x < 1, depleting the network leads to a reduced modulus.By contrast, if x > 1 the stiffness of the remaining bond more than doubles, clarifying our counter-intuitive prediction that at constant stress, network depletion leads to stiffening.
We further confirm these predictions by simulating the fiber density-and stress-dependence of K for networks with various fiber CLs (Fig. 1).Considering first the simple case of regular networks of variable mesh size (Fig. 2G), where all fibers are load-bearing, we recover precisely the scaling behavior predicted by Eq.(3) (Fig. 2H-J) and, in particular, the n-dependence of the large-stress modulus (Fig. 2H).For networks with random depletion (Fig. 2E), the fiber density is controlled by the depletion parameter p, and the connection with the load-bearing fiber density n is less evident.Qualitatively however, we observe that the influence of p is consistent with our prediction (Fig. 2F).For all four CLs considered here, our model adequately captures the stress scaling of the differential modulus K ∝ σ x (Fig. 2A-D).
The scaling of K we observe can be compared with macrorheology experiments that report a stiffening exponent 3/2 for F-actin [45], fibrin [47], vimentin and neurofilaments networks [46] and biomimetic hydrogels [48], whereas Zn 2+ -modified fibrin networks exhibit an exponent 1/2 [66].Our minimal model adequately captures the density dependencies observed for the different stiffening exponents (Fig. 2H).In particular, for x > 1, as for WLC-fibers, the differential modulus decreases for denser networks (Figs.2D,F), in agreement with earlier experiments [45].By contrast, for 0 ≤ x < 1 denser networks display an increased modulus (Figs.2C,F), as for the Zn 2+ -modified fibrin networks [66].Taken together, our results show how at large stress the fiber density controls the mechanical response of soft heterogeneous networks in a way that depends sensitively on the nonlinear micromechanics of the constituents.Interestingly, collagen networks stand out by uniquely displaying a stresscontrolled mechanical response independently of network fiber density.

FORCE ENSURES ROBUST LOCAL RESPONSE
To determine the macroscopic mechanical information a cell can obtain by performing local mechanical measurements inside a disordered network, we assume that cells are ideal mechanosensors: they probe the network by actively exerting a force at the scale of the network mesh-size and measure the network's local compliance without error.To conceptually understand how cells can locally perceive their mechanical environment in a simple way, we study the response of fiber networks to point-like force monopoles.
Experimentally, we perform active microrheology experiments on reconstituted collagen gels.We embed beads that are large enough to be trapped in the collagen mesh (Fig. 3A).Using optical tweezers, we apply controlled forces on these beads to locally measure the mechanical response at various locations in the network (Fig. 3A, see Materials and Methods and [67]).This setup is informative for mechanosensing as it allows us not only to probe the network at the cellular scale, but also to apply forces that are large enough to locally trigger the nonlinear response of a collagen matrix, as observed in the vicinity of cells embedded in fibrous matrices [31][32][33].The force-displacement curves obtained at various locations in the same network (Fig. 3B) reveal two key features: i) the curves exhibit large variability, with a broad distribution of displacements at any force level, and ii) they are nonlinear and exhibit a pronounced stiffening response.For each force-displacement curve F (u) we measure locally, we determine the differential stiffness k = dF/du as a function of force (Fig. 3C).Indeed, while it has been suggested that cells could be sensitive to several quantities such as the strain energy [68], viscoelastic properties [69], and stiffness with extensive evidence [13-16, 33, 70, 71], a complete determination of the mechanical variables cells respond to is still lacking, especially in nonlinear environments.However, as cells have been shown to adapt to the stiffness in collagen matrices, including the local differential stiffness increased by the forces exerted by the cell [33], we here characterize the cell-scale mechanical response in terms of this differential mechanical quantity.
At low forces (F < 0.01pN ), k is independent of F and this linear response is highly heterogeneous.Thus, at low forces cells can only acquire unreliable information about the mechanical landscape of their environment, as previ- a probe (red sphere) in the center of a numerically generated disordered network, inset: tension distribution in fibers located behind the probe (F = 10, same colorcode as in Fig. 1D).Small random sample of measurements performed, E-F: same as B-C for the simulated probings.G: Signal-to-noise ratio (r k ) of the stiffness measurements versus applied force, normalized by the force F * at the onset of the power law regime, for simulated networks with different CLs (Fig. 1) and microrheology measurements on collagen and fibrin gels.
ously observed [21][22][23][24][25].As the force is increased, however, the network stiffens, with locally softer networks stiffening at lower probe forces.Remarkably, at large forces (F ≥ 0.1pN ), the stiffness no longer strongly varies relative to the mean (Fig. 3C).In this nonlinear regime, the local stiffness that cells could measure inside a fibrous matrix becomes reliable: the probing force sets the measured stiffness, as observed for cells in collagen [31,33], and this measured stiffness becomes increasingly robust to local fluctuations in fiber density.
To theoretically understand this robustness of nonlinear stiffness measurements, we employ the model introduced in Fig. 1.We simulate local loading induced by a point-force monopole in a large spherical network with fixed boundary conditions (Fig. 3D, Materials and Methods).To avoid boundary effects and correlations between individual measurements, we perform a local mechanical probe in the center of independently sampled network configurations.For collagen-like fibers with EXP CL, we measure many statistically independent forcedisplacement curves (Fig. 3E) and the resulting differential stiffness (Fig. 3F).We find that these simulations exhibit the same general trend as we observed experimentally: the differential stiffness is highly fluctuating in the linear response regime where soft bending deformations dominate the response [22], but tends towards a single k(F ) curve for a large force.Thus, both our experimental and simulated behaviors reveal that the local mechanical response of a collagen matrix becomes progressively reliable and insensitive to inherent structural disorder in the network when the force triggers the system's nonlinear response.
This may appear unsurprising in light of the fiber density-independence of the nonlinear macroscopic modulus of a collagen network at fixed stress [42] (Fig. 2A, Eq.3).However, to show that the observed macroscopic density-independence of collagen networks does not explain the convergence of the microscopic k(F ) curves, which become independent of local fiber density heterogeneity, we perform microrheological simulations on networks with other fiber CLs (Figs. 1A-C) whose nonlinear bulk stiffness is not independent of average fiber density.Surprisingly, we observe the same features for all CLs (Fig. S5): the local stiffness strongly fluctuates in the linear response regime, while the different differential measurements robustly tend to a single master curve at large forces.
To further quantify this increased robustness, we characterize the ensemble of independent stiffness values at a given force in terms of the signal-to-noise ratio r k = k /std(k) (Fig. 3G).As r k increases with force, the mechanical signal becomes stronger relative to the stiffness heterogeneity.Interestingly, we find that the growth of r k with increasing force can be approximated as a power law in the nonlinear regime: We measure α ≈ 0.6 for both simulated and experimentally measured mechanical responses.This exponent is ) for an applied probe force F = 0.001.E: Radius of the buckling region R * , normalized by the system radius R, versus F applied by a point-force probe.F: Comparison of the measured (dots and diamonds) and expected (dashed lines) signal-to-noise ratio, determined using Eqs.( 4) and ( 6), as a function the size of the probe taken as R * .
also independent of network dimensionality (Fig. S1).In addition, α does not vary substantially when changing the CL: convergence of measurements is observed in all cases according to this universal power law.We experimentally confirm this predicted constituentindependence of the exponent α with microrheology measurements that we perform on a reconstituted fibrin network (Fig. S11), with nonlinear bulk properties that are distinct from collagen [72].Indeed, in fibrin we observe a marked increase of r k in the nonlinear response regime, with a power law increase similar to that observed in collagen (Fig. 3G).Therefore, the origins of this microrheological robustness cannot be the same as for the macrorheological convergence of nonlinear bulk modulus, which is specific to collagen.
In summary, at low force a single mechanical measurement is a poor estimator of the network's average mechanical properties: mechanosensing is strongly limited by structural heterogeneities.By contrast, our experiments and simulations indicate that the local mechanical response of fiber networks becomes largely insensitive to structural disorder at large force.

NONLINEAR MECHANOSENSING MODEL
As the increased robustness of local micromechanical measurements at large force is generically observed for a range of CLs, its physical origin must lie in the fibrous structure of the network, rather than in the specific micromechanical properties of its constituents.This hypothesis is further supported by the modest influence of the CL on differential stiffness in the nonlinear regime: while at the macroscopic scale the bulk modulus followed K ∼ σ x at large stress, the exponent x does not set the microscopic stiffening response [67] (Fig. S5G).Building on recent theoretical results showing that the nonlinear response of fiber networks to force dipoles results in an effective increase of the dipole size [59], we propose that this robustness can instead be understood in terms of an effective increase of the size of the probed region.As network heterogeneities average out on larger scales, this would imply that, as the force increases, the probe becomes less sensitive to local density fluctuations.
To explore this idea, we first examine the linear local response and investigate how the size of the probe affects relative stiffness fluctuations of both 2D and 3D networks with EXP-fibers.While a measurement integrates mechanical contributions at all scales, the rapid decay of linear elastic deformations [73] leads to a dominant contribution of density fluctuations in a small volume in the vicinity of the probe [22].Consequently, individual measurements strongly depend on their location and the corresponding stiffness values display a large variability.A larger probe, however, samples the local mechanics of a bigger region in the vicinity of the network, thereby averaging more effectively over structural heterogeneity.Larger probes are thus more informative about the system's macroscopic response.
To demonstrate that increasing probe size indeed leads to more robust measurements, we perform simulations of a circular rigid body of radius R 0 applying a small force monopole (F = 0.001) to the network.(Fig. 4A).For each probe radius R 0 , we compute the signal-to-noise ratio of linear stiffnesses, revealing a power-law increase with β 0.5 in 2D and β 1.1 in 3D.This power law increase of r k with probe size confirms that the responses depend less on local fiber density fluctuations.
We map this increased robustness of local sensing with probe size to the increased robustness we observe in the nonlinear response for large local probing forces (Fig. 3G).Thus, we argue that a sufficiently large applied force triggers the response of the network over an effectively larger region than in the linear response regime.To determine the force-dependent length-scale that sets the local response, we note that at the onset of the nonlinear response, fibers start to buckle near the probe (Fig. 4B,  left).This buckling spreads to a larger region in the network as the applied force increases [59] (Fig. 4B, right).This leads to the emergence of a 'buckling zone' of growing size R * with a large density of buckled fibers and the formation of tensed rope-like structures.Consequently, within the buckling zone the network strain stiffens as the elastic response is dominated by stretching of the ropes, which is a much stiffer mode of deformation than the fiber bending modes that govern the linear response regime.A probe force deforms both the network inside the buckling zone and the surrounding network beyond R * .These two network sections thus effectively act in series.Because the buckling zone strain stiffens, however, it becomes much stiffer than the network section beyond R * that is still dominated by soft bending modes (Fig. 4C).To first approximation, the buckling zone therefore becomes effectively rigid, and the compliance in response to the probe is dominated by the linear response of the network beyond R * .Thus, elastic nonlinearity renders the network disorder irrelevant inside the stiffened buckling zone, and local stiffness fluctuations are instead determined by network disorder outside the nonlinear zone.Put simply, the emergent length scale R * renormalizes the size of the local probe in a force-dependent way.Local probes with large enough forces thus effectively probe the local linear mechanical properties of the network over a larger length scale making the response less sensitive to local disorder.
To understand the force-dependence of the buckling zone radius induced by monopole probing, we perform an analysis similar to previous work on dipole-induced buckling [57,59,60].Away from the probe, the stress decays as σ(r) ∼ F/r D−1 due to force conservation.From the buckling condition, here written in terms of stress as σ ∼ σ b , with σ b the buckling stress, we expect σ(R * ) ∼ σ b .Therefore, buckling occurs over a region of size with ζ = 1/(D − 1).Indeed, we measure ζ 0.45 in 3D and 0.95 in 2D in our simulations (Fig. 4E, Figs.S7&S8).
To complete our nonlinear mechanosensing model, we now quantitatively connect the increase of r k with probe size in the linear response regime (Eq.( 5)) to the power law increase of r k with the applied point force (Eq.( 4)).
To do so, we identify the effective probe size R 0 induced by nonlinear effects with the buckling length-scale R * .Using Eq.( 4) & (6), our nonlinear mechanosensing model predicts a power law increase of r k with the buckling range R * : r k ∼ R * α/ζ .Importantly, for both 2D and 3D simulated responses the expected power law increase is consistent with Eq.( 5) (Fig. 4F, Supplementary Table 1).Thus, our scaling model establishes that locally probing the network in the nonlinear regime can be conceptualized as a linear probe with a renormalized probing radius associated to the radius of the buckling zone.This model quantitatively explains the CL-independent increase of mechanosensation reliability at large forces in fiber networks.

NONLINEAR MECHANOSENSATION IS RELIABLE
Can cells make use of nonlinear mechanosensation to reliably infer the large-scale mechanical properties of the surrounding matrix from local stiffness measurements?To address this question, we consider local stiffness measurements performed in networks with varying mean fiber density (Fig. 5A), which sensitively tunes the network macroscopic modulus [42, 45-48, 63, 65, 66].We find that in the low force regime, there is a large overlap of the local stiffness measurements obtained on networks with different p values: the variation of individual measurements exceeds the difference in mean stiffness of networks with different fiber densities.By contrast, at large forces the stiffness measured on networks with different fiber densities become clearly separated: by probing the mechanical response with large force, a cell can thus robustly infer the average stiffness of its surroundings and discriminate between the mechanical properties of networks of varying fiber density.
To quantify the ability to infer the macroscopic properties of the network via nonlinear mechanosensing, we compare the local stiffness measured at different forces and the linear bulk modulus K 0 of various networks.While the average local stiffness increases with the bulk modulus, the large scatter of local measurements at low probing forces would prevent a reliable local estimation of the bulk stiffness (Fig. 5B, inset).By contrast, at large forces, the relative mechanical heterogeneity of local measurements diminishes: it becomes possible to discriminate networks with distinct bulk moduli using local stiffness measurements.Indeed, we find that the local stiffness and bulk modulus are increasingly correlated with increasing probe force (Fig. 5B).The Pearson correlation coefficient approaches 1 at large forces: nonlinear mechanosensing enables accurate local inference of bulk network properties.

DISCUSSION
To understand the physical limits of nonlinear mechanosensing we studied the response of disordered networks to local force probes at the cell scale.Using a fiber network model and microrheology experiments, we demonstrated that a local probe can reliably determine both the local and macroscopic mechanical response of the disordered network by triggering either the linear response or elastic nonlinearities.This nonlinear mechanosensation becomes progressively insensitive to network disorder with increasing force.Nonlinear mechanosensing thus offers cells a reliable strategy to both locally determine and control the mechanical properties of a disordered ECM.
We showed that macroscopic and microscopic nonlinear responses of disordered fiber networks are set by distinct mechanisms.Macroscopically, the differential bulk stiffness is a power law of the applied stress in the nonlinear regime (Figs.2A-D), with an exponent x set by the constitutive stiffening of a single fiber (Fig. 1C).We capture this macroscopic behavior with a differentially affine model.This model predicts that the load-bearing fiber density dependence of the differential bulk modulus is set by an exponent 1 − x (Eq.( 3)), which thus solely depends on the constitutive fiber stiffening (Fig. 2H).Remarkably, this implies that the differential modulus of collagen (x = 1) becomes insensitive to fiber density at large stress (Fig. 2A), as observed experimentally [42].Our model further offers insights into experiments with different stress and density-dependencies on various reconstituted networks with distinct stiffening behaviors [42,[45][46][47][48]66]. Microscopically, the local differential stiffness also stiffens with force as a power law (Fig. S5), but with an exponent that is not determined by the constitutive fiber stiffening [67].Instead, this local stiffening is controlled by force-induced buckling and network nonlinearity in the form of a bending-stretching stiffening transition, giving rise to a stiffened buckling zone embedded in a linearly responding network (Fig. 4C).Thus, local force-stiffening is caused by the effective probing of the linear network stiffness at increasingly larger scales set by the buckling zone radius.Our results on the difference between microscopic and macroscopic stiffening mechanisms could be used to further develop accurate approaches for 3D traction force inference or stress inference around cells in 3D matrices [31,33].
The emergence of the stiffened buckling zone explains the mechanical robustness to local force probes.A sufficiently large force effectively probes the linear response averaged over the structural disorder of an enlarged region of the network, with dimension set by the forcedependent buckling length-scale.This elastic regime is relevant for cell-ECM mechanical interactions.Indeed, several cell types apply traction forces of the order of tens of nanonewtons [27][28][29][30]74].These forces can trigger the nonlinear response of reconstituted networks, as shown by our microrheology measurements in collagen and fibrin (Fig. 3B, Figs.S11).This nonlinear response is consistent with observations of contractile cells in a matrix: fibers buckle [31] and displacements are enhanced, decreasing more slowly than predicted by linear elasticity [33,58,75].Furthermore, the network stiffens in the wake of a cell applying traction [31,33,75], and our nonlinear mechanosening model can be used to determine the stiffness such cells would perceive.The magnitude of the force cells need to exert on their substrate to employ nonlinear mechanosensation in vivo likely depends on context, where softer environments typically require smaller forces to trigger nonlinearity.Thus, even though neurons apply smaller forces [76] than e.g., fibroblasts or cardiomyocites [74], differences in stiffness and nonlinear force thresholds of the natural surroundings could still allow such different cells to employ nonlinear mechanosensation.
The formation of the force-controlled buckling zone is qualitatively independent of the local structure of the probe forces.In particular, the characteristic length-scale emerging in response to dipole loading is well characterized [59][60][61].We thus expect a similar robustness increase as a power law for dipole mechanosensors, but with modified exponents.Since cells are mechanically better described as force dipoles [77], this anticipated power law increase in response to dipole loading further supports nonlinear mechanosensing as a cellular strategy.
Many aspects of cellular mechanosensation are still debated, including the internal cellular machinery and pro-cesses that are involved, as well as the mechanical variable that can be sensed by cells [78,79].Here, we characterized the reliability of the mechanical response in terms of stiffness, a mechanical property that is experimentally shown to influence cellular behavior [13-16, 33, 70, 71].Yet, cells are also found to respond to other mechanical variables, such as the substrate strain energy [68].Therefore, one could consider other mechanical quantities to assess the limits of nonlinear mechanosensing, such as strain and elastic energies, which we characterized as alternatives (Fig. S6).In these cases, the signal-to-noise ratio also generically increase with force for networks with a range of fiber constituents.This is understandable since local sensing becomes nonlocal due to the emergence of the buckling zone that facilitates disorder averaging over an increased length scale, regardless of the precise mechanical variable that is considered.Therefore, we argue that enhanced nonlinear mechanosensing is a general characteristic of disordered fibrous networks.
Finally, our work suggests that cells could employ nonlinear mechanosensation as a strategy to reliably sense and respond to the stiffness of their environment.In the linear regime, cells face a highly heterogeneous mechanical landscape [21][22][23][24][25].If this linear regime dominated cell-ECM interactions, we would expect the mechanical heterogeneity perceived by cells to lead to erratic stiffness-dependent cell behaviors following the large local stiffness fluctuations in the ECM.In contrast, if cells trigger the nonlinear response with large forces, then our nonlinear mechanosensation model implies that cells instead face a strikingly homogeneous mechanical landscape that is directly correlated with the matrix' macroscopic modulus, and where stiffness-dependent cell behavior would not depend randomly on the cell's location in the matrix.Nonlinear mechanosensation also has implications for the cell's ability to control the stiffness of their environment.The idea that cells actively stiffen their matrix and adapt in response to the enhanced stiffness has long been introduced [33,72,75], and it was speculated that such a feedback mechanism aims at reaching a specific substrate resistance.We here propose that nonlinear mechanosensation allows cells to exploit this mechanical feedback to accurately control their surrounding stiffness despite the inherent randomness of their local environment, allowing them to robustly perform mechanosensitive cellular functions even in a highly disordered ECM.

Random network generation
Networks are generated by placing straight fibers on an ordered triangular (2D) or Face Centered Cubic (3D) lattice.These networks are randomly depleted with a bond occupation probability p. Unless stated otherwise, we use p = 0.6 in 2D and p = 0.4 in 3D.

Mechanical model
Fibers are discretized with bonds of rest length 0 = 1.Each bond acts as a spring with linear stretching modulus µ = 100 and their nonlinear longitudinal response f ( ij ) is described by a CL displayed in Fig. 1.Fibers also resist transverse deflection with a bending rigidity κ = 1 that penalizes deflections of angle θ.Fibers are connect by freely deforming hinges at their intersection.We consider a probe located on sites i of position R i and applying a force F i .The Hamiltonian of the system is The boundaries of the network are displaced to impose isotropic dilatation.Our depleted networks have dimensions 30 × 30 × 30 and periodic boundary conditions are imposed.The results are averaged over 3 independent random networks.The bulk stress and bulk modulus are computed from the first and second derivative of the system's energy: where γ is the applied dilatation strain and V the system's volume.

Local probing
A point force F is applied to a vertex at the center of spherical (circular in 2D) depleted networks of radius R = 40 in 3D (R = 100 in 2D) for QUAD-, EXP-and WLC-fibers, and R = 45 for LIN-fibers.Finite size effects on the stiffness statistics are displayed in Fig. S3.The probe loading direction is [0, 1, 0], not following any fiber direction in the undeformed network.The network boundaries are fixed.For each CL, 100 disorder realizations are performed.To determine r k vs. R 0 (Fig. 4D), the same boundary conditions apply, a sphere (disk) of radius R 0 is placed in the system's center (R = 30 in 3D, refer to Fig. S4 for a characterization of finite size effects and R = 200 in 2D) and applies a force F = 10 −3 in the direction [0, 1, 0] to all incident bonds, 100 networks realizations are used.

Numerical resolution
At each applied loading, mechanical equilibrium is obtained by minimizing the total energy using the GNU Scientific Library BFGS implementation of the Broyden-Fletcher-Goldfarb-Shanno algorithm.
from the lowest to the highest force, Fig. S2 shows that the obtained mean stiffness and signal-to-noise ratio as a function of force is similar to that obtained in Fig. 2F of the main text over an ensemble of varying size.
Finite size effect: point-force probing Networks for our simulation are generated on a spherical domain of radius R. We verify the effect of the finite system size of 3D EXP-fibers networks on the statistics of the stiffness measurements (Fig. S3).At large forces, the network radius affects the average local response (Fig. S3A).In addition, we observe a sharp increase of stiffness fluctuations that is R−dependent at these large forces (Fig. S3B).These two contributions result in a saturation of the signal-to-noise ratio at forces that vary with system size (Fig. S3C).The radius R = 40 used in the main text ensures a nonlinear response over a large force range before finite size effects set the mechanical response.We note that finite size effects occur at lower forces for fibers responding linearly in tension (LIN-fibers, Fig. S3D-F).We thus use R = 45 for these networks.

Finite size effect: varying probe sizes
We determine the probe size radius at which the finite size of the network affects the statistics of local stiffness values in our simulations.Fig. S4 reveals that as the probe size R 0 exceeds an R−dependent value, the signalto-noise ratio increases more slowly with the probe size.Therefore, the estimate of the β exponent in Eq.( 5) is performed on the largest system of radius R = 30 and for R 0 ≤ 8.

Local probing of networks of different fiber constituents
To evaluate the role of fiber-level nonlinearity in limiting nonlinear mechanosensation, we study the mechanical response to local point-force probing of networks constituted of fibers with a range of stiffening behavior (Fig. 1).For all the CLs considered, the force-displacement individual measurements display a strong variability and a highly nonlinear behavior (Fig. S5A-C).The corresponding local stiffness also behaves as previously observed for 3D EXP-fibers network: there are strong fluctuations in the linear regime and while the local stiffness converge towards a single curve at large force (Fig. S5D-F).The mean stiffness as a function of the applied force reveals a mild effect of the fiber-level stiffening exponent x: LINfiber networks (x = 0) follow k ∼ F 1/2 , as expected from the stiff body of radius R * ∼ F 1/2 probing the linear response of the surrounding network [67].For other constituents with x > 0, the mean stiffness increases slightly faster as x increases, revealing a contribution of mechanical response of the stiff nonlinear region being probed.

ALTERNATIVE MECHANICAL CHARACTERIZATION: NONLINEAR ROBUSTNESS OF THE ENERGY
As cells could be responsive to several mechanical variable determined by probing the mechanical response of the ECM, we also evaluate the robustness of the mechanical response in terms of strain energy W , defined as the area under the force-displacement curve, as an alternative to local stiffness measurements.Considering the force dependency of the local stiffness: the strain energy scales with force as This force dependence of the strain energy is well recovered numerically for 3D depleted networks of EXP-fibers (Fig. S6A).This sample of strain energy measurements as a function of the applied force again highlights the strong heterogeneity of local mechanical measurements and their non-linearity, as observed by characterizing the mechanical response in terms of stiffness.For EXPfibers network, we measure r E ∼ F α W with α W 0.8 (Fig. S6B).Most importantly here, the power law exponent is independent of the fiber CL considered, consistent with our interpretation of the nonlinear mechanical response being set by the response at large scales.
If one discards the loading history and focuses on the elastic energy E = F × u, then using Eq.( 7), the elastic energy also scales with F as which is well recovered numerically (Fig. S6C).According to our effective probe size interpretation, we can use the scaling R * ∼ F 1/2 and the measured signal-to-noise ratio dependency on the probe size r E ∼ R 1 0 (Fig. S6D) to predict r E ∼ F 1/2 .This power law scaling is confirmed numerically (Fig. S6E).

Determination from the compressive bond force
To determine the spatial extent of the buckling region, characterized by a high proportion of buckled fibers, we compute the radius of gyration of the density of buckled bonds.To do so, for a given disorder realization and a given applied force, we divide the network into squares (cubes) of size δ = 1.5 and compute the density of buckled bonds as ρ = N b /N T , where N b is the number of buckled bonds and N T the number of bonds enclosed in that volume.At a given force, the density is averaged over simulations, leading to the example density map shown for 2D networks on Fig. S7A.We then successively determine the center of mass (red cross) and radius of gyration (red circle) of the density and identify the radius of gyration as a measure of R * .The constitutive law of the fibers accounts for bondlevel buckling of the fibers via a softening response in compression (Fig. 1 of the main text).We set the state of the bond to be buckled if the compressive force on that bond is larger than a threshold force f b .The value of threshold does not affect the exponent of the power law increase of R * with F (Fig. S7B-C,E).As displayed in Fig. S7A, the buckling region is delocalized to the back of the probe loading (gray arrow) as a result of the lack of compressive forces transmission ahead of the probe.The nonlinear region thus extends over a distance D * = R * + y CM from the probe, where c y is the location of the center of mass.This length scale approaches the radius R of the system at large forces (Fig. S7D).

Determination from fiber reorientation
The loss of compressive resistance of buckled bonds is accompanied with large bond reorientations and the formation of the rope-like structure over R * .Thus, the buckling region corresponds to a spatial domain where the initial geometry of the network is progressively replaced by a highly reorganized network (Figs.4B).We use this observation to provide an alternative measurement of R * for 2D networks, although the same method can be extended to 3D systems.The statistics of bonds orientation at a given distance from the probe is shown as a function of force in Fig. S8A-C.At each distance from the probe (low to high indicated by dark to bright curve), we determine the force at which the statistics deviates by more than a threshold amount from its value in the linear regime.Considering the deviation of the signalto-noise ratio of the bond orientation r θ by 2.5%, we measure a power law R * ∼ F 1 (Fig. S8D) in good agreement with that determined from the density of buckled bonds.Again, the threshold selected does not affect the power-law exponent (Fig. S8E) and one can perform this measurement on any of the three quantities r θ , θ and std(θ), as shown in Fig. S8E.

EXPONENTS OF THE SCALING ARGUMENTS
The exponents α, β and ζ measured for depleted 2D and 3D EXP-fiber networks are summarized in Table I, where the values expected from our scaling arguments, leading to ζ = α/β, are also displayed.

EXPERIMENTAL MEASUREMENTS
The microrheology measurements are performed on three distinct reconstituted collagen networks and one fibrin gel.

Collagen networks
We label each ensemble of measurements performed on collagen as Set 1 to 3, where Set 3 is reproduced from [31].The measurements performed in each network are shown in Fig. S9 and the signal-to-noise ratio evolution with force of the different sets are compared in Fig. S10.

Outliers removal
For each set, both the force-displacement responses (Fig. S9A-C) and the corresponding differential stiffness measurements (Fig. S9D-F) greatly vary from one location to another.For each data set, one to a few curves stand out and are identified as outliers (red dashed curves).These behaviors can be the result of several experimental limitations.Indeed, it is challenging to control the vertical location of the bead in the network.Yet, a bead located anomalously close to a boundary will lead to a finite size effect that results in an increased stiffening at low forces.The effect of removing these outlying curves when performing the computation of r k is shown on Fig. S10G-I.Discarding these curves does not substantially affect the exponent of the power law increase of r k but ensures that it takes place over an extended force range.Thus, in the main text and in the next section we do not consider these outlying curves in computing the ensembles statistics.

Replicability
To compare the measurements performed on the distinct networks, we show the signal-to-noise ratio r k as a function of force in Fig. S10.The curves are normalized by the lowest value r 0 and horizontal collapse of the curves is obtained by rescaling the force.The overlapping of the curves reveals the replicability of the experiments and strongly supports our theoretical interpretation.

Fibrin networks
The mechanical response of microrheology measurements on fibrin gels are displayed in Fig. S11.We recover the characteristic variability and nonlinearity observed for local measurements in collagen networks and simulated loadings.

MOVIE DESCRIPTION
Our differentially affine model relies on the hypothesis that at asymptotically large stresses a tense subnetwork that carries most of the stress emerges and remains stable under further loading.The emergence and stability of this subnetwork at large stresses is visible in Movie S1 which displays the ension distribution in a 2D depleted EXP-fiber network in response to macroscopic dilatation strain and the corresponding mechanical response.At each distance, these three quantities are normalized by their value at the lowest applied force.D: buckling region radius versus applied force, where R * is defined as the largest distance from the probe at which r k deviates by more than r th k = 2.5% from r0.E: The value r th k considered does not affect the power law increase of the nonlinear region with force.F: R * can alternatively be determined from the variation of the mean bond orientation or fluctuations from their low force value.

FIG. 3 .
FIG.3.Nonlinear response to local probing in disordered fibrous networks.Microrheology experiments on reconstituted collagen networks, A: confocal reflective image of a gel and a trapped 2 µm diameter bead (white dot), inset: schematic of the local probing using optical tweezers, B: applied force (F ) versus measured displacement u of the beads embedded at various locations, B: corresponding differential stiffness (k) versus F .Local probing of 3D depleted EXP-fibers networks, D: a probe (red sphere) in the center of a numerically generated disordered network, inset: tension distribution in fibers located behind the probe (F = 10, same colorcode as in Fig.1D).Small random sample of measurements performed, E-F: same as B-C for the simulated probings.G: Signal-to-noise ratio (r k ) of the stiffness measurements versus applied force, normalized by the force F * at the onset of the power law regime, for simulated networks with different CLs (Fig.1) and microrheology measurements on collagen and fibrin gels.

FIG. 4 .
FIG.4.Scale-dependent nonlinear mechanosensation.Local probing of depleted EXP-fibers networks.A: Tension distribution in 2D networks loaded by probes (gray circular objects) of two different radius R0 (F = 0.001, same tensions color code as in Fig.1D).B: Highlight of the nonlinear region R * (dashed circle) surrounding the point-force probe (arrow).The buckled fibers are shown in green (only on the right half of the figure).C: Schematics of the stiffened buckling zone embedded in a linearly deforming soft network.D: Signal-to-noise ratio (r k ) of stiffness measurements as a function of the probe radius (R0) for an applied probe force F = 0.001.E: Radius of the buckling region R * , normalized by the system radius R, versus F applied by a point-force probe.F: Comparison of the measured (dots and diamonds) and expected (dashed lines) signal-to-noise ratio, determined using Eqs.(4) and (6), as a function the size of the probe taken as R * .

FIG. 5 .
FIG. 5. Robust nonlinear mechanosensing.A: Mean stiffness as a function of force for different values of p, the shaded areas show the stiffness standard deviation.B: Pearson correlation coefficient between the local stiffness (k) and the linear bulk modulus (K0) as a function of the local probing force.Inset: local stiffness measured at 3 different probing forces versus K0.The normalization k is the mean of the measurements computed over all values k measured at a given force level.

2 FIG. S4 .
FIG. S4.Finite system size effect during local probing.Effect of the system radius R of 3D depleted EXP-fibers networks on the local stiffness signal-to-noise ratio as a function of the probe size R0.Model parameters are as described in Materials and Methods.
FIG. S5.Robust nonlinear mechanical response of networks of various constituents.Effect of the fibers constitutive law on the local mechanical response of 3D depleted networks.Force-displacement curves for a small sample of measurements performed for A: LIN-, B: QUAD-and C: WLC-fibers.D-F: corresponding differential stiffness versus force.G: Mean of the local stiffness measurements as a function of the force for the different CLs considered.
FIG.S8.Determination of R * from bonds reorientation in 2D EXP-fiber networks.Considering all bonds located at a distance [r, r + dr[ from the probe, where r varies from 0 (dark blue) to R − 2 (bright yellow) in increments of dr = 2, A: mean, B: standard deviation and C: signal-to-noise ratio of bonds orientation as a function of the applied force.At each distance, these three quantities are normalized by their value at the lowest applied force.D: buckling region radius versus applied force, where R * is defined as the largest distance from the probe at which r k deviates by more than r th k = 2.5% from r0.E: The value r th k considered does not affect the power law increase of the nonlinear region with force.F: R * can alternatively be determined from the variation of the mean bond orientation or fluctuations from their low force value.

1 FIG. S9 .
FIG. S9.Statistics of stiffness measurements in reconstituted collagen networks.For each network, A-C: forcedisplacement response and D-F: corresponding differential stiffness versus force applied by the beads.The red dashed lines indicate identified outliers.G-I: signal-to-noise ratio of the stiffness ensembles when considering all curves (blue) or discarding outliers (orange).

FIG. S10 .
FIG.S10.Replicability of microrheology measurements on reconstituted collagen networks.Comparison of the ensemble of stiffness measurements obtained on three distinct reconstituted collagen networks.For each set of measurement, the normalized signal-to-noise ratio r k is shown as a function of the rescaled force F .Here, r0 is the lowest value of r k and F * corresponds to the force onset of the power law response.