Compressive instabilities enable cell-induced extreme densification patterns in the fibrous extracellular matrix: Discrete model predictions

We present a new model and extensive computations that explain the dramatic remodelling undergone by a fibrous collagen extracellular matrix (ECM), when subjected to contractile mechanical forces from embedded cells or cell clusters. This remodelling creates complex patterns, comprising multiple narrow localised bands of severe densification and fiber alignment, extending far into the ECM, often joining distant cells or cell clusters (such as tumours). Most previous models cannot capture this behaviour, as they assume stable mechanical fiber response with stress an increasing function of fiber stretch, and a restriction to small displacements. Our fully nonlinear network model distinguishes between two types of single-fiber nonlinearity: fibers that undergo stable (supercritical) buckling (as in previous work) versus fibers that suffer unstable (subcritical) buckling collapse. The model allows unrestricted, arbitrarily large displacements (geometric nonlinearity). Our assumptions on single-fiber instability are supported by recent simulations and experiments on buckling of individual beams with a hierarchical microstructure, such as collagen fibers. We use simple scenarios to illustrate, for the first time, two distinct compressive-instability mechanisms at work in our model: unstable buckling collapse of single fibers, and snap-through of multiple-fiber groups. The latter is possible even when single fibers are stable. Through simulations of large fiber networks, we show how these instabilities lead to spatially extended patterns of densification, fiber alignment and ECM remodelling induced by cell contraction. Our model is simple, but describes a very complex, multi-stable energy landscape, using sophisticated numerical optimisation methods that overcome the difficulties caused by instabilities in large systems. Our work opens up new ways of understanding the unique biomechanics of fibrous-network ECM, by fully accounting for nonlinearity and associated loss of stability in fiber networks. Our results provide new insights on tumour invasion and metastasis.


Introduction
Cellular processes constitute a fundamental system of complex cascades of intracellular signalling pathways and biomechanical interactions between cells and their environment.Cells continually remodel the extracellular matrix (ECM) through chemical and mechanical signals [1,2].The interplay between cells and the ECM is of great importance as it regulates a number of cellular processes [3,4], such as cell motility and migration in physiological [5] and pathological [6,7,8,9] conditions, stem cell differentiation [10,11], as well as cell and tissue morphology [12,13,14].The intrinsic actomyosin machinery enables cells to contract, thereby exerting tractions to the fibrous ECM, which results in the generation of spatial patterns of localized deformation [13,15,16,17,18].Although the mechanism for their formation has not been clarified yet, there is sufficient evidence for their implication in intercellular mechanical communication.These patterns are characterized by fiber alignment and severe material densification, localized within tethers that join neighbouring cell assemblies [12,16], such as tumors.The tethers are densified bands where matrix density can be three to five times higher compared to the rest of the matrix.Cells were spotted to leave their cluster and advance along the axis of the tether towards the neighbouring cluster [12,16].In addition, individual cells in collagen induced fiber alignment and network densification along lines connecting them [15], while experiments with isolated fibroblasts reported that cells grew protrusions along the generated tether and towards each other [17].Moreover, cell-induced ECM remodelling underlies additional facets in tissue biology.Cancer studies have demonstrated the preference of tumor cells to invade along densified regions of ECM [8,19,20], while aligned collagen matrix serves as a highway-path which they use in order to migrate [7,8].
In essence, every tissue component-cells and fibers-is a biomaterial with unique mechanical properties that responds accordingly to physical cues.The mechanical behavior of the ECM is thus attributed to that of its individual fibers.Previous studies [21,22] revealed the nonlinear elastic behavior of fibers, which explains why cell-induced deformations extend substantially far and are not confined to the cell boundary [17,23].This nonlinearity is manifested by strain stiffening in tension [22,24,25,26,27] and buckling in compression [28,29,17,30].These effects have been extensively investigated before [17,31,32,33,34,35,36,37,38,39], wherein fibers are modelled as homogeneous beams with stable (monotonically increasing) force-stretch responses.This would seem to suggest that the whole matrix (fiber network) would exhibit stable mechanical behavior.It might come as a surprise that this is actually not true even in the case of linear elastic fibers: in [40] it is shown that the energy of such a network is a nonconvex, multi-well function of nodal displacements, which implies instability.In our recent work [18] instability is also encountered in a continuum model, obtained from orientational averaging of individual fiber response, even when the latter is stable.This model predicts localized densification patterns that bear a strong resemblance with experiments [13,16].
Rather than homogeneous rods, ECM fibers have a bundlelike morphology characterized by a complex hierarchical structure [27,41].This assembly gives rise to unexpected mechanical effects.Subjected to large strains, fibers can be extremely extensible without breaking [41] and they stiffen with increasing tension [25,26].Especially interesting is the case when fibers are subjected to compressive forces under which they buckle, losing stiffness and eventually collapsing.Experiments with elastic fiber networks [42] and later on with fibrin [43] revealed multiple regimes of stress-strain response marked by a non-linear softening instability in compression coupled with network densification.Additionally, uniaxial compression experiments on hierarchical beams [44] revealed a transition from hardening (positive slope in stressstrain response) to softening (negative slope) with increasing compression, while they were reversible upon load release.These studies highlight that fiber mechanics is by far more complicated, exhibiting unstable behavior that deviates from what has been addressed so far towards understanding ECM-related deformations.
Considering the above, our objective here is to develop a discrete network model that accounts for the unique intrinsic features of fiber morphology and mechanics and investigate the role of compression instability in deformation patterns associated with ECM densification, tether formation and fiber alignment.We implement two different families of fiber constitutive relations, with distinct nonlinearity and stability features.Family 1 displays a positive but decreasing stiffness with increasing compression, while Family 2 entails a stretch instability phase, where stiffness becomes negative at extreme compressions.Family 1 represents the traditional view of post-buckling behavior.Family 2 is a more radical model, incorporating recent experimental observa-tions on buckling of hierarchical beams [44].Parts of this work have appeared in [45].
We perform extensive simulations of a fiber network containing one or more contracting circular cells, in the spirit of [17].Simulations with our Family-2 models predict ECM densification, the formation of experimentally observed densified tethers between pairs of contracting cells and enhanced fiber alignment localized within the tethers.We show that ECM densification and fiber alignment are simultaneous consequences of fiber compression instability inherent in Family 2. These phenomena are stronger and highly localized within more sharply defined zones compared to predictions of previous models with stable stretch response [17,32,33,34,35,36,37,38,39].Surprisingly, densified zones with aligned fibers also appear in Family-1 simulations, albeit with important differences.First, they require much higher levels of cell contraction and second, they reveal a different instability mechanism, namely snap-through buckling of larger fiber groups, such as triangular elements.This distinct compression instability mechanism is also evident in simulations with the linear model and has not been addressed by previous studies involving discrete networks.These observations show that the presence of compression instability is critical and essential for localized densification and fiber alignment.

Methods
We have developed a 2D discrete model of a fiber network representing the natural ECM to explore its mechanical behavior under cell induced loading.In particular, we partition a circular domain into triangular elements.Each of the three sides of an element represents an individual fiber (Fig. 1a).Triangle vertices are the nodes of the network, where fibers terminate, so that fiber length corresponds to the length of the segment between two nodes.The nodes are modeled as movable frictionless hinges with two degrees of freedom.In general, nodal displacements change the length of fibers, which are modelled as nonlinear elastic springs.As a result, the total energy of the network (sum of individual spring energies) is a function of all nodal displacements.This leads to the problem of minimizing the total energy over nodal displacements.This problem is solved numerically using advanced optimization techniques; see Subsection 2.3 below.

Mechanical properties
For a single fiber, we introduce the effective stretch λ, which equals the distance between its endpoints divided by its undeformed, or relaxed, length.The energy of a single fiber can be written as W (λ) as a function of effective stretch λ.When the fiber is in tension, it is straight and λ equals the actual stretch (strain +1), while W (λ) equals the elastic energy due to stretching of the fiber.When the fiber is in compression, it may be buckled, in which case the elastic (mostly bending) energy of the fiber can still be expressed as a function W (λ) of the distance between its endpoints, hence of the effective stretch λ.See Fig. 1b.In that case W (λ) is chosen to embody the postbuckling response of the fiber.If the deformed-position vectors of the fiber end points (nodes) are x i and x j and the undeformed fiber length is l ij , the energy of the fiber is the quantity within parentheses above being the effective stretch λ of the fiber between nodes i and j.
We explore various models of the mechanical behavior of fibers, characterized by the force-stretch relation of a single fiber S = S(λ), where is the fiber force, nondimensionalized after dividing by a coefficient with dimensions of force.These models were designed to capture fiber stiffening in tension (λ > 1) [25,26], but also softening in compression (0 < λ < 1) due to buckling [17,31,18].For direct observation of buckled fibers in collagen networks see [28,30].We introduce two families of models.Family 1, shown in Fig. 1c, which includes the linear case (k=1) and Family 2 shown in Fig. 1d: While for all models there is stiffening in tension (except for the linear one S 11 ), the difference between the two constitutive families is in compression.In all nonlinear models in Family 1, force S and stiffness dS/dλ both increase monotonically with increasing stretch λ (Fig. 1c).However, stiffness decreases monotonically with compression (as λ decreases to zero) until it vanishes in the crushing limit λ → 0, whereas force reaches a plateau where it remains approximately constant as λ → 0. Models differ in how abrupt the loss of stiffness is and at which level of compressive stretch it occurs.Thus Family-1 fibers can sustain a limited amount of compressive force even after buckling.
See also [18] where similar fiber behavior is used to derive a continuum model.
In contrast [17] use a bilinear model with piecewise constant stiffness that is lower in compression.Loss of stiffness in Family-1 type force response is intended to model the post buckling behavior (load versus distance between ends) of nonlinear elastic bars in compression by axial loads in the context of large deformations.Family-1 behaviour is consistent with experiments and simulations of post-buckling in certain homogeneous nonlinear elastic beams [46].Actual ECM fibers are far from homogeneous, but exhibit a bundle-like morphology with a complex hierarchical structure [27,41].This gives rise to unexpected mechanical effects.In a recent study [44] with hierarchical beams, uniaxial compression experiments revealed a post-buckling response where the force-stretch relation changes from positive to negative stiffness (slope) for high enough compression.This also occurs in beams composed of metamaterials [46].Family-2 models were designed to capture this instability (Fig. 1d).Stiffness decreases monotonically with decreasing λ < 1 as the fiber initially resists the compressive load, after which stiffness becomes negative with further compression, entering a compression instability regime (negative stiffness) up to final collapse as λ → 0.
Suppose one of the models from (3) or (4) has been chosen.The corresponding elastic energy of a fiber is then given by with a minimum at λ = 1 when the fiber is unstretched.Let λ j be the stretch of fiber j, where j = 1,2,..., F and F is the total number of fibers in the network.Therefore, the total fiber network strain energy is equal to Here F is the number of fibers, N the number of nodes and k ij = 1 if there is a fiber joining nodes i and j, 0 otherwise.

Nonconvexity and Instability
Our initial attempt is to impose displacements or applied forces on boundary nodes of the network and minimize the total network energy E(x 1 , . . ., x N ) with respect to all positions of interior nodes x i .However, allowing for large contractile deformations can result in nonphysical solutions due to interpenetration of matter.This happens when triangular elements fold over and snap through to the other side, as there is no resistance against fibers crossing through each other in the model.Examples are shown in Supplementary SFig.1I.This is an instance of the well-known snapthrough instability of structural mechanics, e.g., [47] and shows that the energy E is nonconvex and likely to have multiple local minima, as well as unstable extrema, even in the case of the linear fiber model S 11 (λ)= λ − 1 as schematically shown in Fig. 2a,b.This and other instabilities also occur in a regular (lattice) network of linear springs [40].
Solutions with snapthrough involve interpenetration of matter and orientation reversal, which are both physically unacceptable.In order to exclude such solutions, an energy penalty term is introduced that resists any two fibers with a node in common from crushing into each other.That would be equivalent to the oriented area of the associated triangular element going to zero, then becoming negative with orientation reversal.Interpolating the deformation from nodes to the entire domain in a piecewise affine way (continuous overall and linear in each triangle), we define J to be the Jacobian determinant of the deformation.Hence J is piecewise constant and equal to the ratio of deformed to undeformed oriented triangle area.Letting x and z be two undeformed vector sides of a triangle, and x, z be the deformed sides, we have in terms of the vector product, where k is the out-of-plane vector.Negative J denotes orientation reversal of the respective elements, i.e. folding over and interpenetration of matter.The penalty term is chosen as a function of J, where b > 0 is a small constant and Q > 0 is large (Supplementary SFig.1II, we usually choose Q=50 and b=1/4).Thus Φ(J) is very small for J > b and becomes very large for J < 0. Thus it serves in maintaining positive orientation in the network, as negative orientation (J < 0) is costly in energy.For elements with positive area ratio, Φ(J) is very small, thereby having essentially no contribution to the network's total energy.Physically it corresponds to fibers (which actually have nonzero thickness) resisting being crushed together when network elements are close to collapsing.The modified network potential energy has the form: where F is the total number of fibers in the network, K the total number of elements, W (λ j ) the potential energy of an individual fiber, A k the reference area that element k occupies and J k its oriented area ratio.Even after modifying the energy by adding the penalty term, instabilities due to nonconvexity are still present.We identify some instability modes next:

Nonconvexity due to Large Rotations
Even if the single fiber energy W (λ) from ( 5) is strictly convex with a minimum at λ = 1, the corresponding energy in (1) is a nonconvex function of nodal positions x i because of rotational invariance (Fig. 4 in [40]).This is a source of nonconvexity of the total energy Ê, that is typically entirely missed when small rotations are assumed [40].

Element Collapse Instability
Before an element undergoes snapthrough (its oriented area changes signs) it buckles, or collapses, when a node touches the opposite side.This is actually an unstable equilibrium of the triangle energy which is thus a nonconvex function of element oriented area ratio J. Compressing an element triangle along its height keeping its base fixed (Fig. 2a), we observe that the energy of the triangle as a function W (J) is minimal and vanishes at J = 1 and at J = −1 (Fig. 2b) after the triangle has snapped through to its mirror image.Since W (J) is odd it must be nonconvex with an unstable equilibrium at J = 0.In our model, such total collapse is prevented by the penalty term Φ(J), but bistability and noncovexity of the penalized energy W (J) + Φ(J) remains (Fig. 2c), with an additional, highly compressed solution for some values of compressive force (Fig. 2d).This occurs in both model Families, even in the linear model; it is an example of the well-known snapthrough instability of structural mechanics, e.g., [47].In order to identify this instability mode in our simulations we define the densification ratio of each triangular element to be the ratio of deformed to reference density of a hypothetical continuum deforming as the triangle.This gives

Fiber Collapse Instability
In Family-2 networks there is an additional instability: when a fiber is compressed past the point where the slope of the S(λ) curve becomes negative (Fig. 1d), it enters an unstable regime, tending to collapse to zero effective stretch.For example, the unstable regime for S(λ) = λ 7 −λ 5 is 0 < λ ≤ 0.85.Clearly, fiber collapse would imply area collapse of any element (triangle) with this fiber as a side (Fig. 2e).Eventually the penalty term ( 8),( 9) restabilizes the element against total collapse.An unstable regime remains in general, rendering fiber compression response essentially biphasic, similar to Fig. 2d.
To summarize, all fiber networks are susceptible to element collapse (triangle buckling) instability.Family-2 networks suffer from an additional fiber collapse instability brought about by total loss of strength due to buckling of hierarchical fibers.Simple geometry shows that fiber collapse implies element collapse (Fig. 2e), but not vice versa (Fig. 2f).

Formulation, Software and statistical analysis
By expressing stretches λ j and Jacobians J k in terms of variable nodal positions x i , i = 1, ..., N , we express the total energy Ê in (9) as a function Ê(x 1 , ..., x N ) of nodal positions.See (6) for the first term.We then perform energy minimization on Ê.
The discrete model has been implemented in Python [48].The triangulation has been implemented in FEniCS [49] and the optimization method (for finding energy minima) is explained in [18].The statistical analysis has been done in R [50] and for multi-group comparisons we used one-way analysis of variance (ANOVA).

Model Geometry
Cells are modelled as circular cavities of radius r c within the domain of outside radius R. Thus, the domain containing the ECM network is an annulus with r c < r < R, where r = |x| is the radial distance from the domain center and x is the position vector.In particular, we model contractile cells.We simulate contraction by prescribing an inward radial displacement of cell boundary nodes given by: u Here λ is defined as the ratio of deformed to reference (l 0 ) distance of a fiber's endpoints.From left to right: a relaxed fiber with length l 0 , a fiber under tension (λ > 1) and a buckled fiber under compression (λ < 1).The cyan arrows represent the applied loads at the fiber's endpoints.(c,d) Force-stretch (λ) curves of individual fibers.(c) Family 1:

Results
3.1.Single-cell simulations 3.1.1.Severe localized densification patterns are observed in Family-2 models, moderate ones for Family 1 We simulate a single cell contracting within a fibrous network for each one of the models introduced in (3), ( 4).Simulations at 50% cell contraction exhibit patterns of highly localized, severe densification shown in Fig. 3a-c, SFig.2a-c.These patterns take the form of bands, emanating from the periphery of the contracting cell into the surrounding matrix.Plotting the densification ratio of each element versus distance from the cell shows that in Family-2 models, highly densified elements have densification ratio ≈ 3 and reach up to six cell radii into the ECM (Fig. 3c).In contrast, Family-1 densified triangles are confined within two cell radii (Fig. 3a,b), with densification ratio at most 2.
The distribution of fiber stretches within the deformed networks illustrates similarities and differences between models (Fig. 3d-k, SFig.2).Fibers under tension (λ > 1) align roughly with the radial direction, forming continual paths that propagate a few cell diameters into the matrix (Fig. 3d-f, SFig.2d-f).This happens regardless of the model, though in Family-2 simulations the paths extend further into the matrix (Fig. 3f, SFig.2f).When it comes to compressed fibers, things differ significantly between models (Fig. 3g-k, SFig.2g-k).Fibers under compression (λ < 1) are oriented close to the angular direction, forming loops around the cell (Fig. 3g-k).Within each of these loops, and close to the cell boundary, the stretch is nearly uniform for Family-1 models (Fig. 3g,h).Similar behavior is seen in [39, Fig. 6].
Simulations with Family 2 exhibit two differences: the distribution of compressive stretch around the cell is strongly inhomogeneous (Fig. 3k), and the maximum compression is up to twice as high as in Family-1 simulations, 60% compressive strain (or stretch λ ≈ 0.4) compared to 30% (λ ≈ 0.7) for Family 1 (colorbars in Fig. 3g-k and SFig.2g-k).Compressed Family-2 fibers are still roughly in the angular direction.The most compressed fibers occur within narrow bands emanating radially from the cell and reaching as far as 6 deformed cell radii into the matrix (Fig. 3k,m, SFig.2k,n).Furthermore, in Family 2, network triangles comprising the densified bands are excessively compressed (Fig. 3m), as they contain fibers that have nearly collapsed.Fibers under tension are aligned along the axis of densification bands, roughly perpendicular to fibers under compression (Fig. 3n).When the densification ratio of the networks in Fig. 3c is compared to the compressed fiber distribution of Fig. 3k, it becomes clear that regions of localized excessive densification ( ≈ 3) coincide with the bands containing severely compressed fibers (Fig. 3c,k-n, SFig.2c,k).
In Family-1 simulations, severe compressive stretch is not observed at 50% contraction level, with λ remaining above 0.7, compared to 0.4 for Family 2. Densified zones are much shorter and confined to the immediate vicinity of the cell (Fig. 3a,b, SFig.2a,b) with triangles less compressed ( at most 2 compared to 3 for Family-2).

Comparing critical levels for densification pattern formation in the two
model families Since the main difference between the two model families is the unstable compression regime in Family-2 models, where the S(λ) curve has negative slope (Fig. 1d), the results of 3.1.1suggest that fiber collapse instability is responsible for the sudden growth of localized densification bands.In order to test this hypothesis, we simulate network response under gradual cell contraction, to study the onset of densified band formation.
Fig. 4a-e shows a Family-2 network with a cell contracting at five consecutive levels, from 20% to 40% reduction of cell initial radius.Initially, as contraction progresses, the densification ratio of essentially the same few triangles proximal to the cell increases linearly with cell contraction (Fig. 4,a-d) up to 35%.Remarkably, at the next level of (40%) contraction, densified bands around the contracting cell have appeared, extending noticeably further into the matrix (Fig. 4e).Below each plot of Fig. 4a-e, in a "tree diagram", we plot the stretch of each individual fiber (abscissa) versus distance from the cell center (ordinate) for each contraction level; color indicates orientation relative to the radial direction.The evident asymmetry near the base of each tree at larger contractions shows the difference of compressive versus tensile stretches.Tensile stretches λ > 1 grow gradually with increasing contraction.In fibers under compression (λ < 1), the stretch first decreases slowly, with only a few fibers in the unstable regime λ < 0.85 (red dotted line), all of whom are close to the cell up to 30% contraction.At 35% there is a steep increase in the number of fibers below the threshold, with stretches down to 0.4 and reaching more than 6 cell radii into the network by 40% contraction.Going back to the respective densification ratio configurations, we observe that the jump in fiber compressive stretch and the abrupt appearance of densified bands occur at the same contraction level between 35% and 40%.
This trend in densification localization is reflected in plots of the maximum over the network of the densification ratio inverse 1/ max , and the minimum stretch λ min , at each contraction level (Fig. 4f).We note that 1/ max = J min , the minimum area ratio, corresponding to the most compressed triangular element.The densification ratio first increases slowly with contraction, then there is a steep rise between 35 and 40% contraction, the level at which extensive localized densification is spotted (Fig. 4e).The minimum fiber stretch follows exactly the same behavior as 1/ max , the two curves in Fig. 4f being nearly identical.Initially, λ min decreases approximately linearly with contraction, namely the minimum stretch occurs at the cell boundary as dictated by the boundary conditions with γ the fractional cell diameter decrease.Then there is a sudden drop in stretch magnitude at 35 − 40% contraction, exactly the level of sudden 1/ max drop in Fig. 4f and band growth in Fig. 4e.This shows that element densification is driven by fiber collapse as explained in Fig. 2e.We recall also Fig. 3m showing a collapsed red fiber within each densified green triangle (See Methods 2.2, Fiber collapse instability) The behaviour of Family-1 networks is different (Fig. 4g, SFig.3,SFig.4).The minimal stretch λ min follows (12) all the way up to the largest simulated contraction (red line in Fig. 4g), occurs on the cell boundary, and is equal to cell boundary contraction prescribed by boundary conditions.This shows that fiber collapse is not observed, as expected.In contrast, the maximal densification ratio does undergo a sudden leap (1/ max drop in Fig. 4g) as in Family-2 models, albeit at a higher contraction level of γ ≈ 45% − 50%.This is evidence of an element collapse instability (See Methods 2.2) that is weaker and requires higher cell contraction than the fiber collapse instability of Family-2 models.Notably, this instability occurs in the linear fiber model S(λ) = λ − 1 (SFig.3)as well as the nonlinear ones (SFig.4).Top: densification ratio color plot at each indicated contraction step.Bottom: tree diagrams, fiber distance from cell center versus fiber stretch for all fibers in the network at each contraction step, x axis: fiber stretch, y axis: fiber distance from cell center.(f,g) Maximum densification ratio inverse 1/ max and minimum stretch value λ min over the network at each contraction level, for Family-2 S(λ) = λ 7 − λ 5 and Family-1 S(λ) = λ − 1 respectively.Note that 1/ max = J min .Red solid line: cell boundary stretch imposed by boundary conditions.Colorbars: Densification ratio and fiber radial orientation (in degrees).

Intercellular tether formation in two-cell simulations
We report on simulations involving a pair of cells contracting at 50% of their initial radius, separated by either 6r c or 4r c , where r c is the cell radius (Fig. 5, SFig.5, SFig.6).What distinguishes these from singe-cell simulations is the spontaneous appearance of intercellular tethers, composed of thin, roughly parallel bands of high densification and fiber alignment, that connect the two cells (Fig. 5c, SFig.5c).When cells are separated by a larger distance, 6r c , tethers are generated only with Family-2 models.Additional densified bands emanate radially from each cell (Fig. 5c, SFig.5c) as before.In contrast, in Family-1 simulations, matrix densification is limited close to the cell boundary and cells remain isolated and disconnected (Fig. 5a,b, SFig.5a,b).When cells are closer together, tethers are generated by all models, even the linear one (SFig.6).In this case, we observe that they are substantially stronger in Family-2 simulations, as they extend from one cell to the other and are noticeably wider compared to Family-1 tethers (SFig.6a-c).
When tethers form, we observe a fraction of fibers, located almost entirely in the intercellular region, to be highly stretched (SFig.6d-f).These fibers are densely packed and aligned with the horizontal direction passing through the cell centers, generating straight paths of fibers connecting the two cells.These paths comprise the tether.At the same time, fibers under extreme compression occupy the same region as the tensile ones, but their orientation is nearly perpendicular to the paths of the tensed and aligned fibers (SFig.6g-k).This is true for all models, though fiber compression magnitude is almost twice as large with Family 2, reaching approximately 70% compression (SFig.6k).This indicates that in Family-2 tethers, compressed fibers are well within the regime of the fiber collapse instability.In addition, we observe highly compressed fibers within the densified bands that emanate radially from the cell periphery (SFig.6k).
When cells are separated by a greater distance, 6r c , and tethers are generated only with Family-2 models, fiber stretches highlight a significant difference between families (Fig. 5d-k, SFig.5d-k).In Family-2 simulations, fiber distributions and orientations within the tether are the same as for shorter distances (Fig. 5f,m, SFig.5f,k).For Family-1 models, this is no longer true, as the fiber paths are disrupted and tensile stretches are distributed in a broader region between cells, without the strong alignment we have with Family-2 models (Fig. 5d,e, SFig.5d,e).This is reflected in angle distributions of the tensile fiber orientation, which are substantially different across models within the intercellular region (SFig.5m).This distribution is more localised for Family-2 models, consistent with greater alignment.
Excessively tensed fiber angles within the densified region are narrowly distributed about zero (horizontal direction through cell centers) (Fig. 5m).Compressed fibers are distributed about 80 − 90 • within the densified region, compared to a uniform distribution in non-densified regions (Fig. 5m).On the contrary, compressive stretches in Family-1 models are confined to concentric loops around each individual cell instead of the region between cells, and oriented in the circumferential direction (Fig. 5g,h, SFig.5g,h).
The previous findings hold for 50% contraction.When cells contract more, tethers are eventually generated for Family-1 models as well.In Fig. 6 we present the case of Family-1 model S(λ) = λ 5 − 1 (eq.( 3), Fig. 1c) with two cells separated by 6r c at four contraction levels 45%, 55%, 65% and 75%.We observe that densification between the two cells progressively strengthens.Fiber compression in the cell-cell vicinity is ever-increasing with contraction level, resulting finally in a solid tether at 75% contraction (Fig. 6d).
Working in the same manner for each model separately, we have tested different contraction levels ranging from 5% to 80% decrease in cell radii, for multiple distances separating the two cells.Results are summarized in Fig. 6e.In particular, for each model we obtain a curve that indicates the minimum contraction cells should undergo to produce a tether, expressed as a function of cell distance.That is, above each curve a tether is predicted to form for the respective model.Clearly, Family-2 models are able to sustain tether formation for moderate contraction levels ≤ 50%, and for relatively large cell-cell distances (up to 11r c ).On the contrary, regarding Family-1 models for the same contraction levels ≤ 50%, a full tether is formed when cells have a distance at most 5r c .
What are the mechanisms responsible for the differences between tethers in the two Family models?Extremely compressed fibers occur in Family-2 tethers (λ min ≈ 0.3) but not in Family 1, where λ min ≈ 0.6 (Fig. 7a,b).In Family 2, fiber collapse (extreme fiber compression, sudden λ min drop in Fig. 7b) occurs at the same time as extreme densification (sudden 1/ drop , Fig. 7b).On the contrary, in Family 1, we see extreme densification without fiber collapse (Fig. 7a).In Family 2, most collapsed triangles within the tether contain a highly compressed fiber, oriented within 45°of the vertical (Fig. 7d as in Fig2e).In contrast, in Family-1, collapsed triangles have nearly horizontal bases, while the other two sides are under moderate compression, and are closer to horizontal than vertical after collapse (Fig. 7c as in Fig. 2f).These findings indicate that fiber collapse instability is the main player in Family-2 tethers, whereas the dominant role in Family-1 tethers is played by element (triangle) collapse instability.

Discussion
Extracellular Matrix (ECM) mechanical remodelling by cellular forces brings about unique deformation patterns of excessive matrix densification and fiber alignment, both playing a central role in intercellular communication [12,13,15,17] and in cell motility and invasion [8,20].In order to explore this type of cell-induced ECM deformations, we develop a discrete model that accounts for individual fibers and their intrinsic mechanics.The discrete fiber network model featured in this study complements our prior work [18] where we demonstrated through continuum modelling and experiments that the aforementioned phenomena are the result of a material instability.The latter is caused by a special nonlinearity due to fiber buckling under compression.To understand how discreteness affects these results, here we implement two different families of fiber constitutive relations, with distinct nonlinearity and stability features.Family 1 displays a positive but decreasing stiffness with increasing compression, while Family 2 entails a stretch instability phase, where stiffness becomes negative at extreme compressions.Family 1 represents the traditional view of post-buckling behavior.Family 2 is a more radical model, incorporating recent experimental observations on buckling of hierarchical beams [44].
Our simulations results have revealed instability mechanisms that have not been identified in previous work.In particular, all nonlinear spring networks are susceptible to element collapse (triangle buckling) instability; see Methods.This includes linear spring networks provided large rotations are accounted for [47].In Family-2 networks, an additional fiber collapse instability occurs because of total loss of strength due to buckling of fibers in compression (Methods).Simple geometry shows that fiber collapse implies element collapse (Fig. 2e), but not vice versa.Thus both instabilities can occur in Family-2 networks.In contrast, Family-1 networks can undergo element collapse, despite the fact that individual fibers exhibit stable behavior.The distinction of these different instabilities is an effect of discreteness and is not captured by continuum models, including ones that allow instability [18].
In Family-2 models we observe a sudden increase in densification simultaneously with abrupt fiber collapse, both in single and two-cell simulations (Fig. 4f and Fig. 7b).The majority of the elements in the densified regions contains a severely compressed fiber (red fibers in Fig. 7d).This is strong evidence that the mechanism behind densification in Family-2 networks is fiber collapse instability (Methods and Fig. 2e), driven by cell compression.Densification also occurs in Family-1 networks, but requires higher levels of cell compression.Fiber collapse is not encountered in this case (Fig. 4g and Fig. 7a).Instead, densification is due to triangular elements collapsing (Fig. 7c).Element collapse instability (Methods and Fig. 2f) is the main player behind the appearance of densified zones in Family-1 networks.These observations apply both to intercellular tethers and densified zones around single cells.
One important finding concerns fiber alignment.Regardless of the model, simulations show that fiber excessive alignment occurs simultaneously with densification and at the same locations in the ECM, both in single and twocell cases (Fig. 3, Fig. 5).In particular, in Family 1 we see that severe element compression forces all sides of the densified triangles to align with each other (Fig. 7c), a direct effect of element collapse instability (Fig. 2f).Especially in two-cell cases, before tethers form we observe a moderate tendency of fibers to align (Fig. 6a,b).After tethers form, fibers within the tethers are aligned almost perfectly (Fig. 6d) and the elements they belong to are the ones that show extreme densification.In Family-2 tethers and bands, we see the stretched sides of the densified elements aligned with each other while the highly compressed one is roughly perpendicular to them (Fig. 3m, Fig. 7d), evidence that fiber collapse instability brings about the alignment of stretched fibers within the densified regions.These significant observations indicate that ECM compression instabilities are responsible for both matrix densification and fiber alignment.
In general it is much easier to form a tether with Family-2 models than with Family-1.For the same distance between two cells, a Family-1 tether requires a much higher compression in order to form.For example, for a distance 6r c , where r c is cell radius, the linear model requires 80% compression whereas 25% is sufficient with λ 7 − λ 5 (Fig. 6e).Given a level of compression, say 50%, a Family-1 tether is formed when cells are very close, less than 5r c .On the contrary, a Family-2 tether can form when cells are more than twice as far away from each other (Fig. 6e).
Our most prominent prediction is the formation of densified tethers between two cells, and densified radial bands emanating from each cell.Zones of high densification (tethers) have been observed experimentally joining two clusters of contractile cells [12] [13] [16] [51] while thinner bands were seen to emerge from each cluster [12] [13], extending and gradually diminishing within the matrix.Densified tethers and radial bands were also observed in [18] where contracting active particles were used in place of living cells, thereby excluding non-mechanical causes behind densification.Fibers within the tethers are highly aligned along the tether axis.Individual cells from each explant [12] [13] or acinus [16] start migrating along the tether in an attempt to reach the neighbouring cluster.Moreover, isolated fibroblasts grew protrusions towards each other, along the tether that formed following their contraction [17].These studies illustrate the significance of tethers and radial bands in cell migration, motility and intercellular communication.Our simulations identify the formation of these densified zones as a direct consequence of compression instabilities.An additional effect of these instabilities is the close alignment of fibers within the densified zones.Alignment and densification are therefore seen to be part of the same mechanism.
The most essential application of this work is cancer invasion and metastasis.Tumor explants cultured in an initially randomly organized matrix aligned the collagen fibers around them by contracting.This allowed individual cancer cells to use the tracts of aligned fibers as highway paths to invade the ECM [7] [8].In fact, both densification [52] and fiber alignment [7] are considered prognostic biomarkers for breast carcinoma [53], specifically, "bundles of straightened and aligned collagen fibers that are oriented perpendicular to the tumor boundary" [53].Contractility is a necessary ingredient in their formation [8].All these observations apply to both tethers and radial densified bands.A different phenomenon is seen in expanding tumors where the densified layer surrounding them consists of fibers parallel (not perpendicular) to their boundary [53].Remarkably, our simulations of an expanding cell confirm this, see Supplementary Material, SFig.7.Additionally, aligned collagen fibers offer biochemical molecules transportation between cells [54].Focusing on the elevated fiber alignment within the tethers, our predictions highlight the contribution of compression instability to cancer related ECM mechanisms.
Family-2 models (unstable stretch responses) result in well-defined tethers with highly localised densification and fiber alignment very close to the tether axis.These results are in qualitative agreement with experiments [13,16,18].On the contrary, in models with stable stretch response (Family 1) tethers are diffused and not localized, while fibers under tension distribute in the broader intercellular region.In addition, tethers with Family 2 occur under experimentally observed physiological levels of cell contraction ≈ 50% [16,17,18], in contrast to Family 1 where extreme contraction is required.Considering the above, Family-2 models are preferable over Family-1 as they describe experimental observations better.The unstable response of Family-2 fibers is justified by recent work [44] on post-buckling behavior of beams with hierarchical structure, such as ECM fibers [27] (see Methods).
There have been considerable efforts in ECM modelling [17,31,32,33,34,35,36,37,39,38].In these studies, ECM fibers are usually modelled as Timoshenko beams [32,36] or as elements with asymmetric elastic responses to extension and compression, obeying either piecewise linear stress-strain curves or combining strain-stiffening with compression softening intended to model fiber buckling [17,31,34,37,39,35,38].Even though these approaches have explored nonlinear aspects of fiber behavior, they are limited by stable stretch responses (monotonic) and sometimes small deformations.As far as the stretch response is concerned, all these models are similar in spirit to our Family-1 models.None of these works address instability.Here, we recognise that instability plays a central role in the appearance of ECM densification and fiber alignment.Our work shows that instability can occur even in models with stable stretch response (Family 1), but it does so at unreasonably high cell contraction levels.This prompts us to introduce models (Family 2) whose stretch response becomes unstable in compression.This particular instability allows tether formation and fiber alignment under experimentally observed cell contraction levels.
Fiber alignment has been explored by previous studies [37], [39], yet as a separate ECM mechanism involved together with compression buckling in intercellular force transmission [37] or matrix elastic anisotropy [39].Sopher et al., [37] suggest that elevated tension in the intercellular region obliges fibers to stretch and align.In our simulations densification and fiber alignment suddenly jump to much higher values compared to aforementioned works.This transition occurs when a compression instability is activated.For models with stable stretch response (Family 1), this would occur in much higher contraction levels than considered in previous works [17,31,34,37,39,35,38], which explains the lower levels of densification and alignment observed there.Family-2 models not only require moderate cell contraction, but also densification is much stronger, tethers are solid and substantially wider and fibers almost perfectly aligned to the tether axis.
Deformation induced anisotropy has been studied [39] as a possible mechanism for long-range cell communication.We remark that the compressive instabilities studied here create strong anisotropy in the densified state of the network, because they create a strongly aligned and dense uniaxial distribution of fibers from an originally roughly isotropic random fiber distribution.
Our results show that there is a single unifying mechanism behind densification, fiber alignment and matrix anisotropy, and this is compression instability.These phenomena occur simultaneously within the same localized zones as soon as compression instability is triggered either due to fiber buckling (collapse) or element collapse (snapthrough buckling).
Compression instability due to buckling of elastic fibers in networks was first identified in [42] as a mechanism of localized densification.The continuum model of [55,18], obtained from a fiber model through orientational averaging, also exhibits a compression instability, and predicts highly localized densified tethers and radial bands.Their morphology is qualitatively similar to the ones reported here.In general, it is practically impossible to obtain an explicit continuum constitutive law for a random network through a rigorous limit as the fiber length approaches zero [56,40].The advantage of our model is that it captures the discrete nature of the actual fibrous network.As a result it can distinguish between different types of compression instabilities (buckling of fibers versus fiber elements) and it clarifies the close connection of instability with fiber alignment and densification.
Possible extensions of our work include three-dimensional network models, more sophisticated modelling of joints and crosslinking between fibers, viscoelastic effects in fiber behavior, and the possibility of larger scale instabilities in analogy with periodic discrete networks [40,57].
Our models highlight compression instability due to buckling as a crucial nonlinear mechanism underlying the mechanical behavior of fibrous ECM and give rise to new insights in exploring the nature of cell-induced deformations that underlie matrix densification and fiber alignment and their implications in intercellular biomechanical interaction, cancer metastasis and cell motility., where J is ratio of deformed to undeformed oriented triangle area.Q > 0 is large and b > 0 is small constant.As a result, negative values of J have high energy cost, whereas positive values have negligible contribution to the network's total energy.Bottom: Simulations of a cell contracting by 50%, either with or without the penalty term for the area ratio J. Without penalizing J, the optimizer finds solutions that are physically unacceptable, as J < 0 corresponds to elements (red) that changed orientation.Colorbar: J values.

Figure 2 :
Figure 2: Instability Mechanisms.(a) Element collapse Instability (snapthrough of triangular elements) under compressive force S (cyan) (b) Energy of a triangular element as a function of its oriented area ratio J.Note nonconvexity and two-well structure.(c) Dotted line: as in (b) for J > 0. Solid line: energy with penalty Φ(J) added.(d) Penalized energy has two stable equilibria J 0 and J 1 under suitable compressive force (equal to the slope of the red straight line) (e) Fiber collapse (red fiber) causes triangular element area collapse.(f ) The converse is not true.Triangular element collapse can happen without fiber collapse.

Figure 3 :
Figure 3: Fiber collapse instability and severe localized densification.Simulations with a single cell at 50% contraction with Family-1 models λ − 1 and λ 5 − 1 and Family-2 model λ 7 − λ 5 (a-c) Densification ratio of triangular elements (color plot) in deformed networks (d-f ) tensile stretches and (g-k) compressive stretches in deformed fibers (m) stretch of deformed fibers and (n) radial orientation distribution of fibers within the densified bands in Family-2 case (c).Colorbars: (a-c) densification ratio ρ of the deformed networks, (d-k) fiber stretch λ.

Figure 4 :
Figure 4: Progressive cell contraction and densification localization.(a-e) Simulations with Family-2 model S(λ) = λ 7 − λ 5 of a cell contracting in the range 5% − 80%.Top: densification ratio color plot at each indicated contraction step.Bottom: tree diagrams, fiber distance from cell center versus fiber stretch for all fibers in the network at each contraction step, x axis: fiber stretch, y axis: fiber distance from cell center.(f,g) Maximum densification ratio inverse 1/ max and minimum stretch value λ min over the network at each contraction level, for Family-2 S(λ) = λ 7 − λ 5 and Family-1 S(λ) = λ − 1 respectively.Note that 1/ max = J min .Red solid line: cell boundary stretch imposed by boundary conditions.Colorbars: Densification ratio and fiber radial orientation (in degrees).

Figure 5 :
Figure 5: Intercellular tether formation.Simulations with two cells contracting at 50%, for three different models (three columns in a-k).Cell centers are separated by 6r c , where r c is the undeformed cell radius.(a-c) densification ratio of triangular elements (color plot) in deformed networks (d-f ) tensile stretches and (g-k) compressive stretches of deformed fibers (m) orientation distribution of deformed fibers within the densified zones (tether and radial bands) in Family-2 case (c) and within the highlighted non-densified zones.Horizontal direction is the one parallel to the axis connecting the cell centres.Radial direction is the one passing through the cell centre.Colorbars: (a-c) densification ratio of the deformed networks, (d-k) fiber stretch.

Figure 6 :
Figure 6: (a-d) Tether formation in Family-1 networks.Simulations with two cells contracting at different levels in a Family-1 network (model S(λ) = λ 5 − 1).Cell centers have distance 6r c , where r c is the undeformed cell radius.Densification ratio color plot of triangular elements (up), tensile (middle) and compressive (bottom) stretches of deformed fibers at each contraction step.(e) Contraction versus cell-cell distance required for tether formation in various models.Simulations with two cells contracting in the range 5% − 80% decrease in cell radius (y axis).Cells are separated by a distance proportional to cell undeformed radius r (x axis).Each curve corresponds to a different model and depicts the minimum contraction level required for a solid tether joining the cells as a function of distance separating them.

Figure 7 :
Figure 7: Mechanisms of densification within tethers in the two Families.(ab) Maximum densification ratio inverse 1/ max and minimum stretch value λ i min over the network at each contraction level excluding fibers on the cell boundary, for Family-1: S(λ) = λ − 1 and Family-2: S(λ) = λ 7 − λ 5 respectively.Note that 1/ max = J min .Red solid line: cell boundary stretch imposed by boundary conditions.(a) Element area collapse (yellow dots) occurs without fiber collapse (blue dots), indicating element collapse instability.(b) Densification (yellow dots) occurs simultaneously with fiber collapse (blue dots), suggesting fiber collapse instability.(c-d) Stretch of fibers (red: compression; green: tension) at 80% contraction within a Family-1 tether (c) and Family-2 tether (d).Note scarcity of compressed fibers despite presence of collapsed triangles in (c), indicating element collapse instability.In contrast, in (d) most collapsed triangles contain a highly compressed fiber, pointing towards fiber collapse instability.

SFigure 1 :
Interpenetration of matter and penalty term.(I) Various presentations of triangulated rectangular truss elements.Each edge in the structures represents a linear spring.Dirichlet boundary conditions were applied on the upper boundary nodes by imposing a displacement u = (h, 0.0), h being the scale to x direction.Deformed structures contain triangles that have changed orientation, resulting in interpenetration of matter.(II) Top: Penalty term Φ(J) = exp(−Q(J−b))

SFigure 4 :
Progressive contraction and densification.(a-b) As contraction level rises, densification strengthens in the close proximity of the cell for Family-1 models.The bands consisting of densified elements do not propagate far from the cell boundary, reaching as far as 3 deformed cell radii at 80%.On the contrary, in Family-2 simulations (c-d) densification is evident at much lower contraction levels, 35%.With increased contraction, more densified bands are generated and extend substantially further into the matrix.Colorbar: densification ratio of deformed networks.