Cells exploit a phase transition to mechanically remodel the fibrous extracellular matrix

TL;DR

This paper reveals that cells mechanically remodel the extracellular matrix through a densification phase transition caused by fiber buckling, which explains the formation of tethers and their role in cell migration and invasion.

Contribution

It uncovers that tether formation results from a phase transition driven by fiber buckling, combining modeling, simulation, and experiments to explain this biological phenomenon.

Findings

Tether formation is a densification phase transition of the ECM.

Simulations predict strain discontinuities localized in tethers.

Experimental active particles mimic cell-induced patterns.

Abstract

Through mechanical forces, biological cells remodel the surrounding collagen network, generating striking deformation patterns. Tethers-tracts of high densification and fiber alignment-form between cells, thinner bands emanate from cell clusters. While tethers facilitate cell migration and communication, how they form is unclear. Combining modeling, simulation and experiment, we show that tether formation is a densification phase transition of the extracellular matrix, caused by buckling instability of network fibers under cell-induced compression, featuring unexpected similarities with martensitic microstructures. Multiscale averaging yields a two-phase, bistable continuum energy landscape for fibrous collagen, with a densified/aligned second phase. Simulations predict strain discontinuities between the undensified and densified phase, which localizes within tethers as experimentally observed. In our experiments, active particles induce similar localized patterns as cells. This shows how cells exploit an instability to mechanically remodel the extracellular matrix simply by contracting, thereby facilitating mechanosensing, invasion and metastasis.

Tables

Table 1. Table S1: Simulations parameters for the Monostable and Bistable Elastic Energy.

ECM Model Parameters
Elastic Energy	Equation (5.2)	Equation (5.3)	Equation (2.4)
Monostable	$a_m=0$	$A=80$, $b=0.22$	$\epsilon =0$ except Fig. 1d see Fig. caption
Bistable	$a_m=0.4696$	$A=80$, $b=0.11$	$\epsilon =0$

Table 2. Table S2: Experimental data. Distance unit: μ ​ m 𝜇 𝑚 \mu m .

Experiment data shown in Fig. 4a.
Temperature	${26}^{\circ }$C	${39}^{\circ }$C
Radius of left particle	75.8	52
Radius of right particle	72.6	49.6
Distance between centers	227.7	219.1
Experimental data shown in Fig. 5.
Temperature	${26}^{\circ }$C	${39}^{\circ }$C
Radius	Fig 5a : 50.1	Fig 5b : 25.3
Radius	Fig 5e : 38.1	Fig 5d : 21.7

Equations

$$\hat{W}({\lambda }_1,{\lambda }_2)=\frac{1}{2\pi }{\int }_0^{2\pi }w\left(\sqrt{({\lambda }_1\cos \theta {)}^2+({\lambda }_2\sin \theta {)}^2}\right)d\theta ,$$

$$W=\bar{W}({\lambda }_1,{\lambda }_2)=\hat{W}({\lambda }_1,{\lambda }_2)+\Phi ({\lambda }_1{\lambda }_2),$$

$$E\{\mathbf{y}\}={\int }_{\Omega }W(\nabla \mathbf{y}(\mathbf{x}))d\mathbf{x}$$

$$G_{\epsilon }\{\mathbf{y}\}=\frac{{\epsilon }^2}{2}{\int }_{\Omega }|\nabla \nabla \mathbf{y}(\mathbf{x}){|}^{2}d\mathbf{x}.$$

$$C\{{\bf y};{\bf z}_{1},\ldots,{\bf z}_{n}\}=\sum_{i=1}^{n}\int_{\Gamma_{i}}{k\over 2}\Bigl{[}|{\bf y}({\bf x})-{\bf z}_{i}|-(R-u_{0})\Bigr{]}^{2}ds_{\bf x},$$

$$\Phi \{\mathbf{y};{\mathbf{z}}_1,\dots ,{\mathbf{z}}_n\}=E\{\mathbf{y}\}+G_{\epsilon }\{\mathbf{y}\}+C\{\mathbf{y};{\mathbf{z}}_1,\dots ,{\mathbf{z}}_n\}.$$

$$S(\lambda )={\lambda }^5-(a_m+1){\lambda }^3+a_m\lambda ,$$

$$w(\lambda )={\lambda }^6/6-(a_m+1){\lambda }^4/4+a_m{\lambda }^2/2+1/12-a_m/4$$

$$\hat{W}({\lambda }_1,{\lambda }_2)=\frac{1}{2\pi }{\int }_0^{2\pi }w(\lambda (\theta ))d\theta =\frac{1}{2\pi }{\int }_0^{2\pi }w\left(\sqrt{({\lambda }_1\cos \theta {)}^2+({\lambda }_2\sin \theta {)}^2}\right)d\theta ,$$

$$\tilde{W}(\mathbf{F})=\frac{1}{96}[5I_1^3-12I_1{J}^2-(1+a_m)(9I_1^2-12J^2)+24a_m{I}_1+8].$$

$$\Phi (J)=\exp [A(b-J)]$$

$$W^{\ast }(\mathbf{F})=\tilde{W}(\mathbf{F})+\Phi (\det \mathbf{F}).$$

$$C\{\mathbf{y};{\mathbf{z}}_1,\dots ,{\mathbf{z}}_n\}=\sum _{i=1}^n{\int }_{{\Gamma }_i}\frac{k}{2}{\left|\mathbf{y}(\mathbf{x})-{\mathbf{z}}_i-(1-u_0/R)(\mathbf{x}-{\bar{\mathbf{z}}}_i)\right|}^{2}d{s}_{\mathbf{x}}$$

$${\Psi }_h[{\mathbf{y}}_h]$$

$$\displaystyle-\sum_{e\in E_{h}}\Big{[}\int_{e}\{\kern-3.125pt\{\nabla\nabla{\bf y}_{h}\}\kern-3.125pt\}\cdot\llbracket\nabla{\bf y}_{h}\otimes n_{e}\rrbracket+\frac{60}{h_{e}}\int_{e}|\llbracket\nabla{\bf y}_{h}\rrbracket|^{2}\Big{]}\Bigg{)}$$

$$\{\{\nabla \nabla {\mathbf{y}}_h\}\}$$

$$[[\nabla {\mathbf{y}}_h\otimes n_e]]$$

$$[[\nabla {\mathbf{y}}_h]]$$

$${\lambda }_2\ge 1,{\lambda }_1\le 1$$

$$g(\phi )=\frac{(1/2\pi )}{({\lambda }_2/{\lambda }_1){\cos }^2\phi +({\lambda }_1/{\lambda }_2){\sin }^2\phi }$$

Full text

Cells exploit a phase transition to mechanically remodel the fibrous extracellular matrix

Georgios Grekas,1,∗ Maria Proestaki,2 Phoebus Rosakis,3,4,∗ Jacob Notbohm,2

Charalambos Makridakis3,4,5 & Guruswami Ravichandran6

1Aerospace Engineering and Mechanics, University of Minnesota, Minneapolis, USA

2Department of Engineering Physics, University of Wisconsin-Madison, USA

3Department of Mathematics and Applied Mathematics, University of Crete, Greece

4Institute of Applied & Computational Mathematics, Foundation for Research & Technology-Hellas, Heraklion, Greece

5 Department of Mathematics, MPS, University of Sussex, Brighton, United Kingdom

6Division of Engineering and Applied Science, California Institute of Technology

Through mechanical forces, biological cells remodel the surrounding collagen network, generating striking deformation patterns. Tethers—tracts of high densification and fiber alignment—form between cells, thinner bands emanate from cell clusters. While tethers facilitate cell migration and communication, how they form is unclear. Combining modelling, simulation and experiment, we show that tether formation is a densification phase transition of the extracellular matrix, caused by buckling instability of network fibers under cell-induced compression, featuring unexpected similarities with martensitic microstructures. Multiscale averaging yields a two-phase, bistable continuum energy landscape for fibrous collagen, with a densified/aligned second phase. Simulations predict strain discontinuities between the undensified and densified phase, which localises within tethers as experimentally observed. In our experiments, active particles induce similar localised patterns as cells. This shows how cells exploit an instability to mechanically remodel the extracellular matrix simply by contracting, thereby facilitating mechanosensing, invasion and metastasis.

1 Introduction

Biological cells remodel the extracellular matrix (ECM) either biochemically, producing and degrading collagen fibers, or mechanically, by exerting mostly contractile forces [1]. Mechanical remodelling of fibrous collagen ECM by cellular forces results in distinctive deformation patterns extending over long distances. If the ECM were a typical elastic material, deformations due to cells contracting would decay with distance within a few cellular diameters [2, 3]. Instead, the observations of Weiss [4] and Harris & Stopak [5, 6] led the way [7, 8, 9] in showing dramatic spatial patterns of densification and fiber alignment, localised within tether-like bands joining distant cell clusters (Fig. 1a). Tethers were recently observed between individual cells [8, 3] as well (Fig. 1b). In addition, radial hair-like emanations from each cluster of cells (Figs 1a, 1c) were observed.These patterns are a central feature of mechanical ECM remodelling and its role in cell motility, invasion and intercellular communication: individual cells leave their cluster and move along a tether to neighbouring clusters [5, 6, 8, 9]; single fibroblasts joined by a tether grow appendages along it toward each other [3]. Cancer cells preferentially invade along densified regions of the ECM with high fiber alignment [10, 11, 12, 13, 14, 15, 16, 17]. Specific patterns of ECM density and alignment have been identified as biomarkers for breast cancer [13, 15, 16, 17, 12].

Along the axis of each tether, the ECM is stretched by as much as 40%, but also compressed in the transverse direction to half or even a quarter of its original thickness [9, 18]. Density can be 3 to 5 times higher within tethers than without. These deformations are not only severe, but also spatially localised within tethers [5], and along radial hairlike bands that issue from individual cell clusters [6] (Fig. 1c). What is the mechanism underlying the formation of these patterns? Many studies focus on the behaviour of collagen in tension, and argue that alignment and densification are induced by tensile strain [8], or are due to a stiffening nonlinearity of the stress-strain behaviour of the ECM in tension [19, 20]. At the same time, density within tethers exceeds the undeformed value by almost an order of magnitude [9, 18], which implies compressive strains at least twice as large as tensile ones in magnitude (Supplementary Materials: Compressive Stretch Estimate). While cells, by contracting, apply radial tensile forces to the ECM they adhere to, they also decrease their perimeter, thereby inducing compression in the circumferential direction [3, 21]. The ECM ligament between two cells, where a tether forms, is thus under axial tension and transverse compression, but the observed compressive strains are unexpectedly large [9]. The degree of localisation of densification is also unexpected. Centimeter-scale tethers induced by multi-cell tissue explants have well defined boundaries, across which density is virtually discontinuous [5, 6] (Fig. 1a). Numerous thinner radial bands issue from each explant. At scales comparable to fiber length, density localisation is gradual but still pronounced [3], whereas fewer radial bands emerge from each cell (Fig. 1b). Solitary cell clusters also produce radial bands (Fig. 1c) [5, 6]. Deformations are strongly inhomogeneous in the angular coordinates and radial symmetry is broken. What is surprising in these observations is the persistent appearance of localised, inhomogeneous deformations, not only at the scale of fiber network inhomogeneity [3], but at macroscopic scales as well [6]. Current nonlinear continuum models do not seem to explain these patterns of localisation [18, 20, 21, 22].

Our model and simulations demonstrate that these phenomena are possible because of a material instability that stems from a special nonlinearity in mechanical behaviour of individual network fibers. Using active particles, our experiments establish that these severe, localised deformations can be brought about simply by mechanical forces, such as the ones due to cell contraction, without involvement of biochemical factors, and they remain after forces are removed. This indicates that they are in fact a form of mechanical ECM remodelling [14, 11, 23].

The connection with instability was first established through experiments and modelling in fiber networks and open-cell foams [24] and later for fibrin [25]. Rather like rubber bands, individual fibers support tensile forces and may stiffen with increasing tension [26, 27], but buckle under compression, losing stiffness and eventually collapsing. Larger polyhedral groups of fibers buckle and collapse under compression [24]. These microscopic buckling instabilities cause the appearance of macroscopic bands of intense compressive deformation and high density, within which fibers are mostly buckled and compacted, alternating with regions of normal density and low compressive strain, where fibers are largely unbent and loosely arranged. These two region types (high- and low-density) are separated by sharp interfaces, across which strain and density jump discontinuously at the macroscopic scale. This behaviour is suggestive of coexistence of a densified phase and an undensified phase. To understand this, we start from the buckling behaviour of a single fiber, and use multiscale averaging to obtain a higher dimensional energy landscape for general deformations. Surprisingly, this energy exhibits bistability in strain. The theory of mechanical, diffusionless, isothermal phase transitions predicts the spontaneous appearance of the distinctive patterns of localised densification and alignment that feature in ECM remodelling and serve as biomarkers for tumour invasion in fibrous ECM.

2 Modelling the Energy Landscape of the Fibrous ECM

We model the nonlinear elastic energy landscape of fibrous collagen ECM, starting from the relation $S=S({\lambda})$ between force $S$ and effective stretch ${\lambda}=d/l_{0}$ for a single flexible elastic fiber; here $d$ is the distance between endpoints, $l_{0}$ is the relaxed (reference) length. See Methods: Model and Supplemental Fig. S1a for details. We assume $S({\lambda})$ stiffens in tension (its slope increases with ${\lambda}>1$) as is common for biopolymers [26, 28, 27]. Aside from direct observation of buckling [29, 30, 31], very little is known about the post-buckling behaviour of individual collagen or fibrin fibers in compression [32], as ${\lambda}$ decreases from $1$ toward [math] (total collapse). This depends on their bending behaviour, which is difficult to characterise, because of their inhomogeneous, hierarchical structure [28] as bundles of loosely connected fibrils. We choose $S({\lambda})$ as in Fig. 2a so that there is stiffening in tension (${\lambda}>1$) but loss of stiffness in compression (${\lambda}<1$) because of buckling. An example is $S({\lambda})=\mu({\lambda}^{3}-1)$ (Fig. 2a). The corresponding elastic energy of the fiber is then $w({\lambda})=\int_{1}^{\lambda}S({\lambda}^{\prime})d{\lambda}^{\prime}$. Suppose the undeformed 2D network has a uniform distribution of fiber orientations. When subjected to a 2D deformation, the stretch ${\lambda}={\lambda}(\theta)$ of a fiber making an angle $\theta$ with the principal axes of stretch in the reference configuration is ${\lambda}(\theta)=\sqrt{({\lambda}_{1}\cos\theta)^{2}+({\lambda}_{2}\sin\theta)^{2}}$, where ${\lambda}_{1}$, ${\lambda}_{2}$ are the principal stretches of the deformation (Supplementary Materials: Principal Stretches). Summing over all orientation angles, following Treloar’s approach [33] and the virtual internal bond model [34, 35], we find a continuum elastic 2D energy density of the network:

[TABLE]

see Methods: Model for more details and an explicit form Eq. (5.2) of the 2D energy.

Instability in Uniaxial Compression.

Surprisingly, the 2D energy density $\hat{W}$ has an instability, evident in Fig. 2b, the stress-strain relation in uniaxial compression $S_{1}({\lambda}_{1})=\partial{\bar{W}}({\lambda}_{1},1)/\partial{\lambda}_{1}$. This relation is not monotonic, but its slope becomes negative for ${\lambda}_{1}<0.45$ (Fig. 2b). This is unexpected, since the original single fiber stress $S({\lambda})$ is monotone increasing. Nonetheless, with increasing compression as ${\lambda}_{1}$ decreases toward [math], the network loses stiffness because of buckling of more fibers. Also, fibers reorient due to compression, increasing the angle they make with the compression axis, so they support less load along the compression axis. As a result, the stiffness eventually becomes negative (strain softening), triggering a compression instability. For a detailed explanation, see Supplementary Materials: Stiffness Loss.

Physically, we expect that increasing compression (${\lambda}_{1}\to 0$) leads to densification, and fibers getting squeezed together, thereby resisting further compression and eventually restoring stability. To account for resistance to crushing (due to fiber volume), we add a volume penalty term to the energy (2.1) that resists collapse to zero volume in the crushing limit as the volume ratio $J={\lambda}_{1}{\lambda}_{2}\to 0$:

[TABLE]

where the penalty function $\Phi(J)$ increases rapidly as the volume ratio $J={\lambda}_{1}{\lambda}_{2}\to 0$ as shown in Supplemental Fig. S1b. We use the specific form given by Eq. (5.3). This restores positive stiffness at extreme compression as ${\lambda}$ approaches zero (Fig. 2c). Thus the decreasing branch with negative slope in the $S_{1}$-${\lambda}$ curve in Fig. 2c, which is unstable, is now sandwiched between two increasing branches. These correspond to two stable phases: a low compression phase and a high compression, densified phase. At a suitable compressible stress $S_{0}<0$ (red line in Fig. 2c and Supplemental Fig. S2b), there are three corresponding stretch states (red and blue dots in Figs 2c, S2b). The middle one is unstable. The states ${\lambda}_{d}$ and ${\lambda}_{u}$ (red dots) are stable. These two states of compressive stretch ${\lambda}_{d}$ (densified) and ${\lambda}_{u}$ (undensified) can coexist side by side in neighbouring parallel bands in the ECM, under the same compressive stress. Such banding deformations with alternating zones of high and low densification, separated by sharp interfaces normal to the compression axis, are observed in open-cell foams and fibrin [24, 25]. This bistable behaviour has much in common with that of shape-memory materials [36, 37].

A Bistable Energy Density.

In the model just described, the densified phase is only stable under compressive stress and disappears upon removal of the latter; this is observed in fully crosslinked collagen [8]. In uncrosslinked collagen, after compression is removed, part of the densified phase remains in the ECM [8, 29, 38, 39, 18]. Various factors contribute to this, e.g., adhesion and new crosslinking of fibers coming into contact due to densification, subject to van der Waals, or other noncovalent types of attraction [40, 38]. To model this, we add a short range attractive adhesion potential $w_{a}({\lambda})$ (blue curve in Fig. 2d) to the single fiber buckling potential $w({\lambda})$ (green curve in Fig. 2d). This represents minimal modelling of new crosslinks developing during fiber proximity [40, 41, 42] that occurs during collapse caused by buckling. Crosslinks can be broken by sufficient stretching, but can be reformed by subsequent compression-induced fiber proximity [43]. This is modelled here by the short-range adhesion potential $w_{a}({\lambda})$. It renders fiber collapse energetically favourable at short distances of the endpoints (${\lambda}$ near zero). The resulting single fiber potential $w_{*}({\lambda})=w({\lambda})+w_{a}({\lambda})$ is now a two-well potential (red curve in Fig. 2d) with a new minimum corresponding to a collapsed fiber state (${\lambda}=0$) that is stable under zero load. Replacing $w({\lambda})$ by $w_{*}({\lambda})$ in Eq. (2.1), we obtain a new bistable energy $W^{*}$ that has multiple global minima (Supplemental Figs S2a, S3a) compared to the single minimum of the monostable energy $W$ (Supplemental Figs S2b, S3b). The new minimum, ${\lambda}_{d}$ in Fig. S2a, corresponds to a crushed (densified) state that is now stable in the absence of compressive stress, modelling the situation where densely packed fibers (collapsed due to buckling) are held together by new crosslinks that develop between them, without the need of external compressive stress.

Nonconvexity in 2D

Next we consider the energy landscape in 2D. What sets both energy densities $W$ and $W^{*}$ apart from previous continuum models of fibrous ECM [18, 21, 20] is nonconvexity (Supplemental Figs S2, S3) and the instability associated with it. This they share with non-convex nonlinear elastic models developed for austenite-martensite phase transformations, twinning, and the shape memory effect [44, 45, 46, 47]. Physically, instability occurs because of stiffness loss and collapse in compression caused by fiber buckling at the microscopic level. Mathematically, this instability is reflected in the fact that both $W$ and $W^{*}$ suffer a loss of a property known as strong ellipticity [48, 49, 44] at some strains, whereby some higher dimensional measure of stiffness becomes negative [50]. See Supplementary Materials: Ellipticity Loss, also Supplemental Figs S3d, S3e.

The bistable energy density $W^{*}$ is a multi-well potential (Supplemental Figs S2a, S3a). This allows the densified phase to be stable at zero stress, as has been observed in uncrosslinked collagen [8, 29, 38, 39, 18]. Although the monostable energy density $W$ has a single minimum (Supplemental Figs S2b, S3b), the associated Gibbs free energy $W({\bf F})-{\bf S}_{0}\cdot{\bf F}$ is a double well potential for a special value of the stress ${\bf S}_{0}$, with minima similar to those of $W^{*}$ (Supplemental Fig. S3c). Thus under suitable compressive stress, the monostable energy exhibits bistable behaviour and coexistence of phases, but predicts the disappearance of the densified phase under zero stress, consistent with crosslinked collagen behaviour [8].

For standard values of model parameters used in our simulations (Supplemental Table S1), the minima (energy wells) of the 2D energy density $W^{*}({\lambda}_{1},{\lambda}_{2})$ occur at the undeformed state $({\lambda}_{1},{\lambda}_{2})=(1,1)$, and at the densified state $({\lambda}^{*}_{1},{\lambda}^{*}_{2})=(0.2,1.06)$. The latter is a severe compression with strain $\varepsilon_{1}={\lambda}^{*}_{1}-1=-0.8$, combined with a small extension $\varepsilon_{2}={\lambda}^{*}_{2}-1=0.06$ in the perpendicular direction. These two energy wells correspond to the undensified phase and the densified phase of the material, respectively. The densified phase involves a five-fold increase in density with ratio $1/({\lambda}^{*}_{1}{\lambda}^{*}_{2})\approx 4.7$.

The energy minima just described satisfy geometric compatibility conditions [45] allowing the two corresponding strain states to coexist in the material, separated by a coherent phase boundary of strain discontinuity [46]. These conditions are described in Supplementary Materials: Compatibility. This means that bands of the densified state can coexist in equilibrium, side by side within the undensified state. An additional minimum at $({\lambda}_{1},{\lambda}_{2})=(0.45,0.45)$ is not compatible with the undeformed state $(1,1)$; see Supplementary Materials: Compatibility. This explains why it is never encountered in our simulations.

Additional Energy Contributions.

The elastic energy of a deformation ${\bf y}({\bf x})$ is

[TABLE]

where $\Omega$ is the undeformed region occupied by the ECM. Going from the discrete energy of a random fiber network of characteristic fiber size $\varepsilon$ to the continuum energy in an asymptotic expansion [51, 34] as $\varepsilon\to 0$, one obtains Eq. (2.3) as the first term, followed by a higher gradient term quadratic in $\varepsilon\nabla\nabla{\bf y}({\bf x})$. We choose a simple form of isotropic higher gradient energy [52, 53]

[TABLE]

To model contractile multi-cell clusters (explants [5, 6], acini [9, 18]), we let $\Omega$ contain initially circular cavities. At the boundary of each cavity, forces can be exerted onto the ECM. In contrast to previous studies [3, 20, 18, 54, 55], cellular clusters are not assumed in our model to remain circular after deformation, and their centres are allowed to move, since they are deformable, and attached to a deforming ECM. Accordingly, each cellular cluster is represented by a distribution of linear springs, connecting each cavity boundary point to a central point which is free to move (Supplemental Fig. S1c). This contributes to the energy a term

[TABLE]

where $\Gamma_{i}$ is the boundary of the $i^{\rm th}$ cavity ($i=1,\ldots,n$), $k$ is the spring stiffness, ${\bf z}_{i}$ the variable centre position of the $i^{\rm th}$ cluster, $R$ the undeformed radius, and $u_{0}$ the cellular spring contraction. Our soft model for clusters allows them to contract and exert radial forces to the ECM but change shape as well. The model mimics the radial contracting actomyosin network of contractile cells, e.g., [56]. Simulations yield deformed cluster shapes consistent with observations: egg-shaped and pointed towards the tether (observation Fig. 1a reproduced from [6]; simulations Fig. 1e). Our active particles are an order of magnitude stiffer than typical cellular clusters [57, 58]. They are modelled as stiff spheres, so that they do not deviate much from sphericity; see Methods: Model for Active Particles, Eq. (5.5), Supplemental Fig. S1d, and Supplementary Materials: Simulation Parameters.

The total energy to be minimised is the sum of Eqs (2.3), (2.4) and (2.5):

[TABLE]

3 Experiments and Simulations of ECM Deformations Induced by Active Particles

How can we ensure that the densification observed in experiments [8, 9, 18] is not due to biochemical factors (e.g., cells depositing more fibers) but entirely mechanical? In these experiments, cells are concentrated within well defined clusters while the densification takes place, and they only migrate away from these clusters after densification has occurred. So there are no cells (that could deposit new fibers) in the densified zones until the latter have fully developed. That the densification is associated with compressive strain can be directly deduced from [​[9] supplemental video sm14], where the evolving compression is recorded and the compressive strain can be estimated from the deformation of gridlines deposited on the ECM.

To further ascertain that tether formation and concomitant densification and strain localisation can occur solely due to mechanical forces, not to other biochemical factors, we embedded active hydrogel microspheres [31] into Collagen I instead of cells. These PNIPAAm particles undergo a temperature-induced phase transition, causing their radius to contract by as much as 60% when heated above $32^{\circ}$C (Methods: Experiments with Active Particles). Advantages of this are that we can control the amount of contraction via temperature, and even cause reverse re-expansion by cooling, without recourse to chemical means. See Methods: PNIPAAm Particle Generation for more details. As we report below, the densification patterns caused by active particle compression are consistent with those caused by cell clusters [8, 9, 18], providing strong evidence for their mechanical origin.

After developing a finite element scheme that can handle deformation gradient discontinuities (Methods: Numerical Method) we minimise the total energy, Eq. (2.6), with respect to the deformation field ${\bf y}({\bf x})$ and the particle/cluster center positions ${\bf z}_{i}$. The simulated ECM domain contains one or more initially circular cavities of radius $R$, representing cell clusters, or the active contracting particles in our experiments. Choosing $u_{0}>0$ in Eq. (2.5) contracts the natural length of the springs comprising the cell cluster model from $R$ to $R-u_{0}$, thus exerting contractile centripetal forces onto the cavity boundaries. For simulations of active particles we use Eq. (5.5) in the energy (2.6).

Tethers.

The most striking feature of our numerical solutions involving two contracting cavities is the spontaneous appearance of a tether, a zone of high density joining the two clusters [5, 6, 9, 18], and thinner hairlike bands emerging from each cluster in the radial direction [5, 6], tapering off and terminating within the domain. This is the ubiquitous morphology reported previously (Figs 1a-1c) and encountered in our simulations for a wide range of model parameters (Figs 1d-1f). In these simulations, within each tether and radial band, stretches are in the densified phase; outside they take values in the undensified phase. Density is discontinuous across the boundary of tethers and radial bands; the ratio of densities outside and inside the tether is in the range 3-5. Within each tether, there is tension along the tether axis and compression in the transverse direction. The compressive stretch is discontinuous across the tether boundary and as low as $20\%$ (compressive strain is as high as $80\%$); the tensile stretch is much smoother although it is higher within tethers, with tensile strains as high as $30\%$. This is consistent with stretch values observed in experiments [9, 18]. The large discontinuous jump in compressive stretch is related to the fact that energy bistability occurs in compression. Compatibility of discontinuous strains plays an important role. Essentially, compatibility restricts possible discontinuities of strains (or stretches) across a phase boundary to ensure that the displacement remains continuous across it. It allows the compressive stretch normal to the tether boundary to be discontinuous across it, but the tangential tensile stretch is restricted to be continuous, in order to ensure displacement continuity across a phase boundary (see Supplementary Materials: Compatibility). This is consistent with experiments [[18] Figs S2, S5] and is observed in our simulations (Figs 3a, 3b).

How close do two active particles have to be (Fig. 4), and how much must they contract in order to form a tether joining them? We performed multiple simulations of a pair of particles. We independently varied particle contraction and distance between particles. The simulations provided a separatrix curve of average particle radial strain versus distance between particles (blue curve in Fig. 4g). Above this curve, our model predicts that a tether forms joining the two particles; below the curve no tether will form. Data from our experiments agreed with this prediction: blue points in Fig. 4g are data points from particle pairs with a tether observed joining them, red points correspond to pairs without a tether between them.

Shape Change and Particle Motion.

Departing from common practice [3, 20, 18, 54, 55], we allow clusters to change shape and move during ECM deformation in our model. See Eq. (2.5), also Methods: Model for Active Particles and Supplemental Fig. S1d. In experiments [5, 6], isolated explants remain roughly circular (Fig. 1c) after contraction. In contrast, neighbouring clusters connected by a tether lose circularity [5, 6] and become egg-shaped with pointed ends toward each other (Fig. 1a) due to tether tension. The parameter controlling shape change for clusters and particles is the stiffness $k$ in (2.5), (5.5). Active particles are an order of magnitude stiffer than cell clusters [57, 58]. For low values of cellular spring stiffness, $k=1$ in Eq. (2.5), our model predicts similar shape changes with pointed end where the tether makes contact with the cluster (compare 2nd and 3d panels in Fig. 1e to shape of explants in Fig. 1a). In order for simulated explants to remain circular after contraction as typically assumed, we must choose an artificially high stiffness ($k=100$ in Eq. (2.5), 1st panel in Fig. 1e). See Supplementary Materials: Simulation Parameters for stiffness values.

Contracting active particles in our experiments (Figs 4a, 4d) and simulations (Figs 4b, 4e) move towards each other by as much as 10% of their original distance (Fig. 6g), while remaining nearly spherical. These particles are roughly an order of magnitude stiffer than cell clusters. To simulate them, we use a high spring stiffness value $k=10$ (compared to $k=1$ for cell clusters) in (5.5).

Tether Splitting.

In our experiments, we observe differences in morphology of tethers between active particles. A tether is sometimes separated into thinner parallel bands (Fig. 4a). Other tethers are in the form of uniform bands making full contact with particles (Fig. 4d). Mammary acini [8, 9, 18] are typically joined by uniform tethers. The tether continues into a layer of the densified phase enveloping each acinus.

Why this difference in morphologies? Unexpectedly the answer comes from our understanding of austenite-martensite transformations, where the Bain strain at the martensitic energy minimum [36] is incompatible with zero-strain austenite [45, 46] (see Supplementary Materials: Compatibility). Namely, these two minimal-energy strains cannot occur on either side of a phase boundary (strain discontinuity) without causing a mismatch in displacements. This forces splitting and tapering of twin bands in a crystal near an incompatible boundary [59]. An example of such split martensitic twins is shown in Fig. 4c. Here also, energy minimisation compels strains not to stray far from energy-density minima. The active particles in Fig. 4a contracted by $u_{0}/R=38$%. The azimuthal stretch ${\lambda}_{\theta}=1-u_{0}/R=0.62$ imposed at the particle boundary by contraction is incompatible with the stretch ${\lambda}^{*}_{1}=0.2$ corresponding to the densified-phase energy well. To avoid this mismatch while maintaining displacement continuity, the tethers splits into narrow bands to minimise contact with the particle boundary (Experiment:Fig. 4a, simulation 4b). The particles of Fig. 4d contracted more ($52\%$). The resulting tether does not split, but is in full contact with the left particle (torn fibers near the right particle keep the tether from extending to the left particle boundary). A corresponding simulation (Fig. 4e) confirms the absence of tether splitting near the particles for this level of contraction. For even higher contraction levels of $80\%$, possible for acini [8, 9, 18] but not for active particles, the azimuthal stretch ${\lambda}_{\theta}\approx{\lambda}^{*}_{1}$ (is compatible with the energy minimal stretch ${\lambda}^{*}_{1}$). In hypothetical simulations of this contraction level, the densified phase is in full contact and envelopes the particle because of enhanced compatibility (Fig. 4f). These observations support the hypothesis that tether morphology is decided by strain compatibility as in other types of coherent phase transitions (martensitic).

Contracting vs Expanding Particles.

What causes the densified phase to split into thin radial bands around particles, such as “hairs” emanating from clusters in Fig. 1c (reproduced from [5, 6]) and in Fig. 5b? A contracting inclusion induces radial tension and circumferential compression. Since the energy is bistable in compression, phase change tends to occur along the direction of compression, with phase boundaries normal to it, roughly along radial lines (Fig. 5b). Instead, an expanding particle would create radial compression, hence a circumferential phase boundary and densified ring enveloping the particle (Fig. 5e). The stark contrast between the ECM’s response to contracting and expanding particles is seen in our experiments (contracting particle Fig. 5b vs expanding Fig. 5e) and for the first time captured by a continuum model simulation (contracting Fig. 5c vs expanding particle Fig. 5f). See also the densified ring surrounding a re-expanded particle in experiment Fig. 6c, simulation 6f). Remarkably, this explains the sharply defined circumferential layer of densification typically observed surrounding expanding tumour spheroids; see, e.g., [​ [14] Figs 1e,f].

Splitting Hairs. Numerical Mesh Dependence.

As noted above, for contracting particles, the radial orientation of interfaces forces the layer of densification around each particle to split into roughly radial bands. Another factor that promotes even finer splitting is compatibility of strains, in order to lower the cost of matching displacements at the particle boundary. This is the same phenomenon as splitting of a tether into finer bands (see Tether Splitting above). These effects are responsible for the radial hair morphology around each particle (Figs 1, 4). As a result, each contracting particle is surrounded by a mixture of the undensified and densified phases, the latter occurring within thin, roughly radial, hairlike bands. These spatially fine phase mixtures establish a further connection with the nonlinear elastic theory of phase transitions in crystals [45, 46]. In multi-well energy minimisation problems for martensitic transitions, it is known that the energy cannot reach a global minimum because of strain incompatibility [45, 60]. Decreasing the energy is facilitated by refinement: increasing the number of bands in alternating phases, while decreasing their thickness, in effect decreases the mismatch in displacements that costs energy to overcome. In theory, the energy approaches an infimum only in the limit of infinite refinement. This explains observed finely twinned martensitic microstructures [45, 61]. In computations [62, 63] this phenomenon manifests itself as an increase in the number and spatial frequency of twin boundaries with increasing computational finite element mesh resolution. Similarly, when we decrease the computational mesh size in our simulations, the number of radial bands issuing from each cluster increases and their thickness decreases (Fig. 1f).

To show that this finite element mesh dependence is not an artifact of the numerical method, we add the higher-gradient term Eq. (2.4) to the energy. This introduces a length scale proportional to $\varepsilon$, an additional material parameter related to characteristic fiber length, bending stiffness and other fiber network parameters [34]. For larger $\varepsilon$, there are fewer, thicker radial bands (Fig. 1d). Above a critical value of $\varepsilon$, they disappear, but the tether persists. The presence of this term has a smoothening effect as it penalises high strain gradients; strain discontinuities are replaced by transition layers with thickness of order $\varepsilon$ [52, 53]. This eliminates mesh dependence at numerical mesh sizes below $\varepsilon$, which controls the scale of the phase mixture. In practice we expect $\varepsilon$ to be of the order of the free fiber length in the network. Accordingly, many hairlike bands are experimentally observed (Fig. 1a, 1c) issuing from millimetre-size explants [6], (where $\varepsilon$ is very small compared to explant size), while only a few (Fig. 1b) from single micron-scale fibroblasts [3].

Residual Densified Zones and Inelasticity.

Do the densification patterns in the ECM persist after cell-exerted contractile forces cease [8, 18]? In our experiments, active particles (Fig. 6a) contracted upon heating to $39^{\circ}$C, causing densified tethers and radial hairs to appear (Fig. 6b), then expanded to their original radius upon cooling back to $26^{\circ}$C (Fig. 6c) thus providing a controllable method of performing one—or several—loading/unloading cycles. When contraction was reversed in our experiments, some residual tethers remained (Fig. 6c), but they were thinner and less prominent than ones appearing during particle contraction (Fig. 6b). Many radial bands issuing from particles largely disappeared.

We performed contraction-expansion cycle simulations (Figs 6d-6f) of an active particle pair undergoing quasistatic, gradual contraction, followed by gradual expansion to original size, using the bistable energy density (we matched initial diameters, distance and contractile strains of both particles; see Table S1). Experimental observations (Figs 6a-6c) confirmed by our model simulations (Figs 6d-6f) included: (i) a residual tether remaining upon re-expansion, which is weaker, thinner and disjointed compared to the tether appearing during contraction (ii) most radial bands disappearing, (iii) a new circumferential layer of densification appearing upon re-expansion on each particle boundary (compare simulation Fig. 6f to experimental Fig. 6c).

Is plasticity of collagen necessary for tether/densification pattern formation [8, 18, 22]? To answer this, we also performed cycle simulations using the monostable energy density $W$. Typical tethers and radial hairs appeared during particle contraction, but disappeared upon re-expansion to original size. This shows that the appearance of localised densification bands is due to the elastic microbuckling instability, which is accounted for in the monostable energy density and does not require inelasticity. In our model, whether these patterns persist once contractile forces are removed, depends on whether the densified phase is stable under zero stress. This is the case for the bistable energy density, with its equally stable minima, but not for the monostable energy density.

This difference is consistent with experiments in fully crosslinked collagen, where the densified phase is only observed during application of compressive stress [8], versus uncrosslinked collagen, where part of the densified phase remains after compression is removed [8, 29, 38, 39, 18]. Analogously, some austenite-martensite transitions occur under stress only, with martensite disappearing upon its removal; this behaviour is called pseudoelasticity [36, 37], described by monostable-type nonconvex elastic energy models. Under different circumstances, such as temperature, martensite persists after load removal, as captured by a bistable energy [36, 61, 64]. The description of both the pseudoelastic and residual transitions, and the present densification transition, does not require plasticity type models, as hysteresis and residual deformations are accounted for by phase transition models [36, 61, 64]. This is demonstrated by our simulations Figs 6d-6f.

Tether-Crack Interaction and Tether Networks.

Our model captures complex experiments of Shi et al [9], where a cut (crack) is made between two acini in order to interfere with tether formation (experiment: dotted line in Fig. 7c, reproduced from [9]). The original tether disappears; instead tethers form that bypass the crack by going around its corners (arrows in Fig. 7c). Our simulation (Fig. 7d) captures this phenomenon, with predicted tethers bypassing the opened crack.

Multiple contractile cell clusters or acini generate a network of tethers, which brings about extensive remodelling of the ECM. See., e.g., [​​ [9], Fig. 1D]. A simulation of our model with a hypothetical distribution of contractile acini is shown in Fig. 7. In Fig. 7b the deformed positions of acini are superposed on the reference positions to demonstrate considerable motion of particles due to matrix deformation. This simulation suggests that multiple acini contracting undergo more dramatic motion and distortion than isolated pairs of acini.

4 Discussion

The continuum model we have proposed differs from previous ones, because it is based on a nonlinear elastic, nonconvex strain energy density, that features a regime of instability in compression. This unstable strain regime separates two distinct stable strain regimes, the undensified and densified phases of the material. The macroscopic instability stems from the microbuckling behaviour of individual fibers. It is the central player in the prediction of the densification patterns due to contractile cells and clusters. While various previous models (including our own) have incorporated fiber buckling in different ways [3, 21, 55], they have restricted attention to small deformations, and have not taken full account of kinematic nonlinearity. Other nonlinear elastic models allow large deformations but assume a phenomenological behaviour that is stable to begin with [20]. Some models have focused on tensile stiffening behaviour of network fibers [19]. In most cases, the unstable strain regime involving compression/densification is not probed by continuum or discrete models so far, e.g., by discrete models studying stress amplification and force chains [54, 55]. Other models investigate nonaffinity in network displacements due to factors such as network randomness, usually restricted to small deformations [65]. The present model averages over this type of nonaffinity by coarse-graining from a discrete network to a continuum description, but identifies a more severe type of nonaffine deformations: strain localisation due to a mechanical instability. These deformations are unusually large but observed in tethers. This requires us to retain full nonlinear kinematics and to consider arbitrarily severe strains and large rotations.

Plasticity has been proposed as a mechanism for tether formation, e.g., [18, 22]. On the other hand, based on experiments, Vader et al. [8] conclude that “the reversibility of fiber alignment and gel densification, seen in crosslinked collagen samples, show that these effects are primarily elastic”. Our simulations with the monostable energy model agree with this conclusion for crosslinked collagen. It seems that this behaviour is in fact pseudoelastic [36, 37], namely due not to plasticity but to an elastic instability. The macroscopic instability is caused by microbuckling and reflected by nonconvexity of the elastic energy landscape.

For uncrosslinked networks, residual tether-like deformations are predicted by a plasticity model incorporating crosslinking [18], but also by the bistable version of our model. Nonetheless, there are various characteristics observed in our experiments and in those of [18], that are not captured by plasticity models. These include virtually discontinuous compressive stretch (but smoothly varying tensile stretch) in the densified zones, and the split tether and split radial hairlike band microstructures seen in our experiments. Our phase transition model does predict these numerically and explains them theoretically as arising from geometric incompatibility of energy minimising strains with particle imposed deformations. Our bistable energy model does account for new crosslink formation and captures the residual densification after force removal observed in our experiments. This is due to stability of the densified phase under zero stress, analogous to residual martensitic transformations associated with the shape memory effect that are not due to plasticity [61].

Our model simulations tend to overestimate the length of radial hairlike bands (Figs 4-6) and the change in distance between particles, Fig. 6g, caused by deformation from our experiments. On the other hand, the model only involves five constitutive parameters (Table S1), less than half the number employed in plasticity based models.

How is tether formation involved in aiding and abetting cell migration, invasion and intercellular communication? After acini contract causing a tether to form, individual cells from each acinus start migrating along the tether, towards the acinus at the other end [5, 6, 8, 9]. Isolated fibroblasts grow appendages toward each other along the tether that forms (Fig. 1b) as a result of their contraction [3]. Given this observed tendency of cells to approach each other, the advantage offered by the biphasic behaviour of the fibrous ECM becomes clear: to detect its neighbours, all a cell has to do is contract, uniformly in all directions, without prior cues as to the direction of neighbours. The automatic response of the ECM is to form sharply defined, densified paths (tethers) leading directly to nearby cells or clusters. Our model identifies this as a special feature of the instability-driven multiphase behaviour of the ECM, not possible in stable, single-phase elastic materials.

Two ECM-related factors identified as breast cancer risk biomarkers are density [66, 17] and fiber alignment [12, 15, 16, 13]. Both arise in our model and experiments, as a result of phase change brought about by buckling fiber collapse. The densified phase in our model and experiments, also in [9, 18], is characterised by a 3- to 5-fold increase in density, but also strong fiber alignment toward the direction of maximum stretch, as compression squeezes an originally uniform angular distribution of fibers away from the compressive direction and toward the orthogonal tensile direction (see Supplementary Materials: Fiber Alignment in the Densified Phase). Hence, within a tether formed between contracting particles (Fig. 6e) the fiber angle distribution we predict, Eq. (S2), is as in Fig. 3e, with a peak along the tether axis, normal to the contracting spheroid boundary. This compares well with Fig. 3f (reproduced from [15]) which shows a typical fiber distribution associated with TACS-3 [15, 16, 13, 12] or Tumour-Associated Collagen Signature 3 (aligned fibers normal to contractile tumour spheroid boundary) a biomarker for invasion of tumour cells along tracts of aligned fibers. Tether formation exhibits the defining features of TAC3; accordingly, when contractility is inhibited, tethers typical of TACS-3 are suppressed [16]. At earlier stages when a tumour is growing and pressing against the surrounding ECM, our simulation Figs 5f,3c,3d of expanding particles (Fig. 5c) predicts a fiber distribution Fig. 3g peaked in the direction tangent to the (expanding) boundary. This agrees qualitatively with collagen signature TACS-2 [15, 16, 13, 12] (fibers compacted and aligned parallel to a growing tumour boundary) with measured distribution Fig. 3h (reproduced from [15]).

Our work provides strong evidence that two different TACS in ECM associated with tumour growth and invasion are distinct manifestations of a phase change in the ECM, caused by tumour growth (TACS-2 [15, 13]) or tumour cell contraction (TACS-3) [16].

Other than serving as a track for cell migration, the aligned and densified collagen fibers within tethers can enhance anisotropic transport of macromolecules and extracellular vesicles. This can enhance cancer-stroma communication and facilitate tumour progression [67, 68].

Taken together, our experiments, model and simulations show that the densification patterns caused by contractile cells and cellular clusters remodelling fibrous ECM exhibit many salient features of stress-induced phase transformations in solids, captured by nonlinear models featuring a nonconvex strain energy [44, 45, 46, 47, 48, 53, 60, 61, 62, 63, 64, 59]. Remarkably, we find that such nonconvex energies arise naturally in modelling the mechanics of fibrous biomaterials, once the microbuckling instability mechanism is accounted for. Our model is minimal, with a handful of parameters, and demonstrates the necessity of a paradigm shift: material instability has to be taken seriously if we are to understand the behaviour of fibrous networks, such as collagen ECM. This is what enables the model to capture the distinctive morphology characteristic of mechanical remodelling of the fibrous ECM for the first time, in the form of cell-induced, two-phase, pattern-forming deformations. We expect it will provide new insights into the role played by these singular deformations in cell migration, communication, tumour cell invasion, and alteration of mechanical ECM properties caused by cell-induced remodelling.

5 Methods

Experiments with Active Particles

PNIPAAm Particle Generation.

PNIPAAm particles were generated using two different recipes, one giving larger particles (diameter 150 $\mu$m) and the other giving smaller particles (diameter 75 $\mu$m). Although these sizes are greater than the size of a cell, they could match the size of a cell cluster or mammary acinus. More importantly, the contraction of the particles create controllable localised forces that bring about the densified zones of interest to this study. To generate the larger particles, kerosene with 3.5% Span 80 (Tokyo Chemical Industries) was degassed under vacuum for 1 hour. It was then maintained under nitrogen for 10 minutes before stirring at 350 rpm at 22∘C. A solution containing 1 g N-isopolyacrylamide (Sigma), 7.5 mL of 2% bis-acrylamide solution (Bio-Rad), 0.05 g Ammonium Persulfate (Bio-Rad) and 1.5 mL of 1 $\times$ tris-buffered saline was prepared. 10 $\mu$l TEMED (Bio-Rad) was added, and the solution was immediately added to the kerosene. To generate the smaller particles, cyclohexane, rather than kerosene, with 3.5% Span 80 was used as the solvent, and stirring occurred at 450 rpm instead of 350 rpm. A solution of 1 g N-isopolyacrylamide (Sigma), 3.75 mL of 2% bis-acrylamide solution (Bio-Rad), 0.05 g Ammonium Persulfate (Bio-Rad), 1.5 mL of 1 $\times$ tris-buffered saline and 3.75 mL deionised water was prepared. Again, 10 $\mu$l TEMED (Bio-Rad) was added, and the solution was immediately combined with the cyclohexane. For both large and small particles, the solutions were stirred for 30 min, and then the PNIPAAm particles were allowed to settle overnight. The particles were then washed with hexane and again allowed to settle. Washes were repeated with isopropyl alcohol, ethanol, and finally deionised water. The particles were then filtered with a cell strainer to keep particles of diameter $>$ 40 $\mu$m. The particles were resuspended in 1 $\times$ PBS.

The PNIPAAm particles were embedded in networks of rat tail collagen I (Corning), which was fluorescently labeled with Alexa Fluor 488 (Thermo Fisher Scientific) [31]. To adhere the particles to the collagen, the particles were first treated with sulfo-SANPAH (Proteochem) as in our previous studies [31]. Next, the pH of the collagen was neutralised using 100 mM HEPES buffer to a final concentration of 3 mg/mL and mixed with the PNIPAAm particles. The collagen then polymerised at 27∘C for 1 hour. After polymerisation, the collagen formed a network of fibers that branch into multiple fibers at discrete nodes. The average distance along a fiber between nodes, which we refer to as the fiber length, was measured as in our previous work [39]. Briefly, images of the fibers were segmented [69], giving segmented areas of the void space between fibers. The average of the segmented areas was 3.9 $\mu$m2. Hence, the mesh size, which is proportional the square root of the segmented areas, was 1.97 $\mu$m. Next, we accounted for the fact that some fibers in the image cross without connecting [70], by selecting representative fibers and manually counting the number of unconnected crossings. The product of the mesh size and number of unconnected crossings gave the fiber length, of 19.7 $\pm$ 2 $\mu$m (mean $\pm$ standard deviation).

Microscopy and Imaging.

Images of PNIPAAm particles embedded in collagen networks were collected using a spinning disk confocal microscope (Yokogawa CSU-X1) with a Nikon Ti-E base and a 20 $\times$ 0.75 NA air objective (Nikon) using a Zyla sCMOS camera (Andor). The particles diameters and the distances between particle centers were measured manually using ImageJ.

The temperature was controlled with an H301 incubator (Okolab) mounted on the microscope stage and controlled with an UNO controller (Okolab). To verify that the samples reached the desired temperature, we used a digital thermometer having accuracy of 0.1∘C (Fisherbrand Traceable) with its probe inside a dish with water placed next to the samples. Images were collected when samples were at 26∘C (undeformed state) and 39∘C (contracted state).

Model

Energy Density.

For a single fiber, we introduce the effective stretch ${\lambda}$, 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({\lambda})$ as a function of effective stretch ${\lambda}$. When the fiber is in tension, it is straight and ${\lambda}$ equals the actual stretch (strain $+1$), while $w({\lambda})$ equals the elastic energy due to stretching of the fiber. When it is in compression, it may be buckled, in which case the elastic bending energy of the fiber can still be expressed as a function $w({\lambda})$ of the distance between its endpoints, hence of the effective stretch ${\lambda}$. See Supplemental Fig. S1a. In order to model a 1D two well energy $w({\lambda})$ for a single fiber like the red curve in Fig. 1d, we start with the derivative $S({\lambda})=dw({\lambda})/d{\lambda}$, which represents force as a function of stretch. We choose a polynomial that vanishes at [math], $1$ and an intermediate value,

[TABLE]

to qualitatively represent the curve in Fig. 2e. Here $a_{m}$ is a parameter that controls the relative height of the two wells (minima) of $w({\lambda})$. Integrating this with respect to ${\lambda}$ gives an energy

[TABLE]

We model the ECM as a 2D nonlinear elastic continuum undergoing possibly large deformations ${\bf y}({\bf x})$ where a particle with position vector ${\bf x}$ in the undeformed state is mapped to deformed position ${\bf y}={\bf y}({\bf x})$. The elastic strain energy density of the material can be written as a function $W({\bf F})$ of the deformation gradient ${\bf F}=\nabla{\bf y}$, which is a $2\times 2$ matrix. We model the random fiber network as an isotropic material, which means that $W$ depends on ${\bf F}$ only through the principal stretches ${\lambda}_{1}$, ${\lambda}_{2}$, whose squares are the eigenvalues of the Cauchy-Green deformation matrix ${\bf F}^{T}{\bf F}$. Equivalently $W$ is a function of the two deformation invariants $I_{1}({\bf F})=\hbox{tr}({\bf F}^{T}{\bf F})={\lambda}_{1}^{2}+{\lambda}_{2}^{2}$ and $J({\bf F})=\det{\bf F}={\lambda}_{1}{\lambda}_{2}$. Here $J$ is the Jacobian determinant of the deformation and the ratio of deformed density $\rho$ to undeformed density $\rho_{0}$ satisfies $\rho/\rho_{0}=1/J$.

To connect the single fiber energy with the 2D strain energy density $W$ we follow [34, 35]. We assume the network contains a uniform distribution of fibers. A homogeneous deformation (with constant strain) is equivalent to a biaxial stretch in two orthogonal directions with stretches ${\lambda}_{1}$, ${\lambda}_{2}$. A fiber of undeformed length $l$ that makes an angle $\theta$ with the principal stretch axes in the undeformed state, will have endpoints at $(0,0)$ and $(l\cos\theta,l\sin\theta)$. After deformation the latter will become $({\lambda}_{1}l\cos\theta,{\lambda}_{2}l\sin\theta)$. As a result, the stretch ratio of the fiber will be ${\lambda}(\theta)=\sqrt{({\lambda}_{1}\cos\theta)^{2}+({\lambda}_{2}\sin\theta)^{2}}$, and its energy will be $w({\lambda}(\theta))$. Summing over all fiber orientation angles $\theta$, we find the elastic energy density of the network

[TABLE]

which gives Eq. (2.1). It turns out that for $w({\lambda})$ given by Eq. (5.1), the integral in Eq. (2.1) can be evaluated explicitly:

[TABLE]

Here $I_{1}=\hbox{tr}{{\bf F}^{\rm T}{\bf F}}={\lambda}_{1}^{2}+{\lambda}_{2}^{2}$ and $J=\det{\bf F}={\lambda}_{1}{\lambda}_{2}$ are the 2D deformation invariants. We then add a fiber volume penalty term to the energy to account for resistance of densified fibers to complete crushing by virtue of their nonzero volume. This term increases the energy abruptly when the Jacobian (volume ratio) $J=\det{\bf F}$ becomes less than a small positive constant $b<<1$, while it becomes negligible as $J$ increases from $b$. Such a function is given by

[TABLE]

(Supplemental Fig. S1b) where $A$ is a large positive constant. The total energy is

[TABLE]

Model for Active Particles.

We consider two models, one for a cell cluster (e.g., acinus) and another for active particles. The model for acini is given by Eq. (2.5); see also Supplemental Fig. S1c. Active particles are modelled as spheres whose radius can undergo a stress-free transformation strain due to temperature (Methods: Experiments with Active Particles), but they can also deform in response to boundary forces. This deformation is modelled minimally using linear springs on the boundary of the sphere. Each point on the sphere circumference is connected to the matrix by a linear spring of zero natural length.

[TABLE]

where $\Gamma_{i}$ is the boundary of the $i^{\rm th}$ cavity ($i=1,\ldots,n$), $k$ is the spring stiffness, ${\bf z}_{i}$ the variable (deformed) centre position of the $i^{\rm th}$ cluster, $\bar{\bf z}_{i}$ is its undeformed position, $R$ the undeformed radius, and $u_{0}$ the particle radius contraction, with percentage referring to $u_{0}/R$. See Supplemental Fig. S1d.

Numerical Method

We employ the finite element method, based on a triangulation of the domain $\Omega$. Our approximations for the deformation are sought in the space of continuous piecewise polynomial functions of degree two. Consequently, our computational methods are based on a discrete minimisation problem on the finite element space. In order to capture areas of high densification accurately a local mesh refinement strategy close to the areas of phase transition is adopted.

Challenges include the subtle non-linear character of the problem, and the nearly singular behaviour of solutions in areas where phase transitions take place. It is known that numerical algorithms can become quite subtle exactly at these areas, and thus special care should be given to the reliable resolution of the interfaces.

Incorporating higher gradients into the energy functional, $\varepsilon>0$ in Eq. (2.4), introduces additional challenges, because the finite element spaces based on piecewise continuous polynomials have reduced smoothness and are not consistent with the standard energy setting of the model. Furthermore, it is very desirable to have a computational model that works seamlessly when introducing regularisation by higher gradients. We thus use the same discrete spaces in all models considered in our study ($\varepsilon=0$ and $\varepsilon>0$). To this end, we adapt to our problem an approach based on the discontinuous Galerkin methodology. Our approximations are still sought on the same spaces of piecewise continuous polynomial functions over a triangulation of the domain, however the energy functional is modified to account for possible discontinuities of normal derivatives at the element faces. We introduce a novel discrete energy functional, which includes terms accounting for the jumps of the higher gradients at the element interfaces, as well as penalty terms which enforce weak continuity of the higher gradients and coercivity. To be more specific, let ${\bf y}_{h}$ denote a function of the discrete finite element space of piecewise polynomials ($C^{0}$-conforming), then the discretised functional for has the form:

[TABLE]

where $E_{h}$ is the set of the interior facets of the triangulation, $h_{e}$ is the length of the edge $e$ and the average and jump operators $\{\kern-3.125pt\{\cdot\}\kern-3.125pt\}$, $\llbracket\cdot\rrbracket$ are defined as follows

[TABLE]

here the superscripts $+$ and $-$ indicate functions evaluation on opposite sides of an edge $e$, $n_{e}^{+}$, $n_{e}^{-}$ are the corresponding outward normal to the edge and $\otimes$ denotes the tensor product.

The discretisation has been implemented in FEniCs [71]. For the minimisation of the discrete energy functional a parallelised nonlinear conjugate gradient method [72] has been developed. The reliability of the computational experiments is guaranteed by a separate detailed mathematical study [73], which demonstrates the convergence of the discrete numerical solution to the solutions predicted by the model.

The energy functional is minimised using a parallelised version of the nonlinear conjugate gradient method [72], a successful iterative method for large scale nonlinear optimisation problems. When ellipticity of the corresponding Euler-Lagrange equation for the unregularised problem holds, Newton’s method is an efficient nonlinear minimisation technique. However, because ellipticity fails in our model (and phase transition occurs) the Hessian matrix of second derivatives of the energy is not positive definite, thus the nonlinear conjugate gradient is preferred, for both the unregularized and regularised energy functionals.

Author Contributions

PR, GR and JN planned the research. GG and PR developed the model. GG and CM developed the numerical method. GG performed the simulations. MP and JN designed the experiments. MP performed the experiments; MP and JN analyzed the data. All authors discussed and analyzed the results and contributed to writing the manuscript. Conflict of Interest. The authors declare no conflict of interest.

Funding

The work of PR, CM and GG was partially supported by the EU Horizon 2020 Research and Innovation Programme under the Marie Sklodowska-Curie project ModCompShock (modcompshock.eu) agreement No 642768. The research of GG was also supported by a Vannevar Bush Faculty Fellowship. The work of JN and MP and was partially supported by National Science Foundation grant number CMMI-1749400. GR acknowledges the support of the National Science Foundation (DMR No. 0520565) through the Center for Science and Engineering of Materials at the California Institute of Technology.

Data Accessibility

The numerical code used in the simulations reported herein as well as experimental data are available at https://github.com/ggrekas/rsif-2020-0823.

Acknowledgments

We thank Brian Burkel for assistance in microscopy.

SUPPLEMENTARY INFORMATION FOR:

Cells exploit a phase transition to mechanically remodel the fibrous extracellular matrix,    J. R. Soc. Interface 20200823

Georgios Grekas, Maria Proestaki Phoebus Rosakis, Jacob Notbohm,

Charalambos Makridakis & Guruswami Ravichandran

Supplemental Notes

Supplemental Figures S1-S3

Tables S1, S2

Supplemental Notes

Compressive Stretch Estimate.

To estimate compressive strains we note that in the tether the density ratio $\rho/\rho_{0}$ of deformed to undeformed density is typically in the range [18] 3-5 also observed in our experiments. From mass balance we have that $\rho{\lambda}_{1}{\lambda}_{2}=\rho_{0}$. Thus the inequality $\rho/\rho_{0}\geq 3$ and the fact that the tether is in longitudinal tension so that ${\lambda}_{2}\geq 1$ give the inequality ${\lambda}_{1}<1/3$. This is roughly 70% compressive strain $\varepsilon_{1}={\lambda}_{1}-1\approx-0.67$. This is consistent with the compressive stretch we estimated from [​[9] supplemental video sm14], where the evolving compression is recorded and the compressive strain can be obtained from the deformation of gridlines deposited on the ECM.

Stiffness Loss.

For a network of fibers with force-stretch relation as in Fig. 2a, we show that the orientation averaged stress goes to zero in the crushing limit. Consider compression along the $y$-axis with stretch $0<{\lambda}_{c}<1$ combined with tension along the $x$-axis with stretch ${\lambda}_{t}\geq 1$. Then a fiber that makes an angle $\theta$ with the compression ($y$)-axis in the undeformed state, will make an angle $\theta_{\ast}$ such that $\tan\theta_{\ast}=({\lambda}_{t}/{\lambda}_{c})\tan\theta$. Since ${\lambda}_{t}/{\lambda}_{c}>1$ we have $\theta_{\ast}>\theta$ and the angle of the fiber with the direction of compression increases and approaches $90^{\circ}$ in the crushing limit as ${\lambda}_{c}\to 0$. This reorientation tends to decrease the component of the load in the compression direction a fiber can sustain. For example, consider uniaxial compressive stretch with ${\lambda}_{t}\geq 1$ fixed and ${\lambda}_{c}\to 0$. The fiber stretch is $\sqrt{({\lambda}_{c}\cos\theta)^{2}+({\lambda}_{t}\sin\theta)^{2}}$ The component of fiber axial force in the $y$ direction of compression is $T=S\cos\theta_{\ast}$. Now $S$ is bounded in compression and as ${\lambda}_{c}\to 0$, $\theta_{\ast}\to\pi/2$ so $\cos\theta_{\ast}\to 0$. As a result the component of the fiber force along the compression direction $T\to 0$. The compressive stress is $\int Td\theta=\int S\cos\theta_{\ast}d\theta$. Now $S\cos\theta_{\ast}$ is bounded for compression and $\cos\theta_{\ast}$ converges (almost everywhere) to zero (except when $\theta=0$) so bounded convergence implies that $\int Td\theta\to 0$ as ${\lambda}_{c}\to 0$. As a result the material fails to sustain uniaxial compressive stress in the crushing limit, and the qualitative behaviour is necessarily as in Fig. 2b.

Single Fiber behaviour.

Modeling the force-stretch relation of a single fiber depends on the post-buckling behaviour of a flexible beam, which can be quite complicated as buckling is a bifurcation/instability phenomenon. There are typically three main types of post-buckling response $S({\lambda})$, depending on whether stability is maintained or lost on the bifurcating branch of buckled states [32], and whether it is later regained if lost: (i) Stable post-buckling, where stiffness is diminished but load carrying capacity in compression is not lost (as ${\lambda}$ decreases below $1$). An example is $S({\lambda})=k({\lambda}^{3}-1)$ with $k>0$ a constant (Fig. 2a). (ii) Unstable (collapse) post-buckling, where the stiffness becomes negative (stability is lost) with increasing compression, e.g., $S({\lambda})=k({\lambda}^{5}-{\lambda})$. (of the form shown in Fig. 2b). (iii) Snap-through post-buckling, where stability is first lost but eventually regained at higher compression [32], for example $S({\lambda})=k\left({\lambda}^{5}-3{\lambda}^{3}/2+{\lambda}/2\right)$ (qualitatively as in Fig. 2e).

Ellipticity Loss.

The strong ellipticity condition [49] plays a central role in coherent phase transformations in solids. It is closely related to the rank-one convexity condition [48] and implies local stability of the material. In one dimension it is equivalent to positive slope (tangent modulus) of the (nonlinear) stress-strain curve, so it would fail for some intermediate strains where the slope is negative in Figs 2c, 2f. In higher dimensions its meaning is more complex [50]. An elastic energy density function that is globally strongly elliptic cannot sustain strain discontinuities and phase transitions [44]. In order to determine for what values of the principal stretches strong ellipticity fails in our model energy density functions $W$ and $W^{*}$, we use the criteria of Knowles & Sternberg [49]. The result is shown in Supplemental Figs S3d, S3e.

Compatibility.

The compatibility condition for strain discontinuities states that for two different values of the deformation gradient ${\bf F}_{1}$ and ${\bf F}_{2}$ to occur on either side of a strain discontinuity surface, across which the displacement is continuous (such as a coherent phase boundary or twin boundary), they must be rank-one connected, namely satisfy ${\bf F}_{2}-{\bf F}_{1}={\bf a}\otimes{\bf n}$ for some vector ${\bf a}$, where ${\bf n}$ is the unit normal to the surface [44]. In case ${\bf F}_{1}={\bf 1}$ (undeformed state) then the principal stretches of ${\bf F}_{2}$ must then satisfy the inequalities [45, 46]

[TABLE]

It turns out that in our model the minimal energy principal stretches $({\lambda}^{*}_{1},{\lambda}^{*}_{2})$ at the energy well corresponding to the densified phase always satisfy these inequalities. In the specific example of the bistable energy $W^{*}$ used in our simulations, $({\lambda}^{*}_{1},{\lambda}^{*}_{2})=(0.2,1.06)$ (Fig. S3a). Since the energy is a symmetric function of $({\lambda}_{1},{\lambda}_{2})$, the minimum at $({\lambda}^{*}_{2},{\lambda}^{*}_{1})$ corresponds to the same state as $({\lambda}^{*}_{1},{\lambda}^{*}_{2})$. The additional minimum at $({\lambda}_{1},{\lambda}_{2})=(0.45,0.45)$ is not compatible with the undeformed state $(1,1)$, since it violates Eq. (S1). The principal stretches in our simulations are never observed to take values at or near this state, This further illustrates the role of compatibility.

In Fig. 3 the phase boundary between the tether and the rest of the ECM is roughly parallel to a principal direction of stretch. In that case compatibility requires that the stretch in the direction parallel to the phase boundary has to be continuous across it (in this case the tensile stretch ${\lambda}_{2}$, Fig. 3b, whereas the compressive one ${\lambda}_{1}$ (Fig. 3b) may be discontinuous as seen in the simulation.

Principal Stretches.

An isotropic elastic energy density function $W({\bf F})$ in 2D, normally expressed in terms of the deformation gradient matrix ${\bf F}$, can be written as a function of the principal stretches ${\lambda}_{2}\geq{\lambda}_{1}>0$ which are the eigenvalues of the right stretch tensor ${\bf U}=({\bf F}^{T}{\bf F})^{1/2}$. Typically in a tether between similar clusters, there is moderate tension along its axis and sever compression normal to it, thus the principal axes of stretch roughly correspond to these two directions. Because of compatibility, the compressive stretch ${\lambda}_{1}$ can jump across the tether boundary, but the tensile one cannot, as it stretches the material roughly parallel to the boundary. This is encountered in experiments [18] and captured in our simulations as seen in Fig. 3

Fiber Alignment in the Densified Phase.

By considering how a homogeneous deformation changes fiber directions, it is easy to show that an undeformed uniform distribution of fibers becomes a distribution

[TABLE]

where ${\lambda}_{2}\geq{\lambda}_{1}>0$ are the principal stretches and $\varphi$ is the angle between the deformed fiber and the maximum stretch direction in the deformed state, or eigenvector of the left Cauchy-Green tensor ${\bf F}{\bf F}^{T}$ with eigenvalue to ${\lambda}_{2}^{2}$. This distribution becomes uniform when ${\lambda}_{1}={\lambda}_{2}$, whereas for ${\lambda}_{2}>{\lambda}_{1}$ it has a peak at $\varphi=0$, the direction of highest tension, Fig. 3. To derive this, consider a biaxial stretch with principal stretches ${\lambda}_{2}\geq{\lambda}_{1}>0$ along orthogonal directions. Let $\varphi$ be the angle between the deformed fiber and the deformed direction of largest stretch. Then $\varphi=f(\theta)=\arctan({\lambda}_{2}\sin\theta/{\lambda}_{1}\cos\theta)$ where $\theta$ is the corresponding undeformed angle. The referential distribution is then $g_{0}(\theta)=df(\theta)/d\theta$ and the deformed distribution $g(\varphi)$ satisfies $\int g(\varphi)d\varphi=\int g_{0}(\theta)d\theta$, which after changing variables from $\theta$ to $\varphi=f(\theta)$ and some algebra gives Eq. (S2), with a constant factor ensuring that $\int_{0}^{2\pi}g(\varphi)d\varphi=1$.

In the densified phase, stretches are close to the energy minimum at $({\lambda}^{*}_{1},{\lambda}^{*}_{2})=(0.2,1.06)$. For these values the deformed fiber distribution is as in Fig. 3e. This situation occurs within tethers. Here the angle is plotted relative to the particle boundary so $\varphi=90$ is along the tether axis.

Simulation Parameters.

The model parameters used in all simulations of our own experiments are listed in Supplemental Table S1. In (5.1) we chose the smallest even powers that will allow the S-shaped form shown in Fig. 2e. There are 4 model parameters. The value of $b$ in (5.3) was chosen to control the density ratio of the densified phase to roughly equal 5, which is within the range 3-6 reported by [18] and deduced from [[9] supplemental video sm14]. Parameter $A$ was given a high value to ensure a steep enough barrier against total volume collapse that at the same time does not affect the energy for larger volume ratios. Parameter $a_{m}$ in (5.1) was chosen so as to result in equal heights of the energy wells (densified and undensified phases) of the bistable energy $W^{*}$ in (5.2), (5.4).

The Young’s Modulus of the PNIPAAm active particles with the same concentration of N-isopropylacrylamide as used here, but a lower concentration of crosslinker, was measured to be $>20$ kPa [57]. Also the modulus of acini such as the ones used by [9] is reported in the range 1–4 kPa [58], so that active particles are roughly an order of magnitude stiffer than cell clusters. Accordingly, in our simulations we take the bulk modulus of active particles to be ten times larger that that of acini. In terms of nondimensional spring stiffness, this translates to $k=10$ for active particles in (5.5), whereas for acini we take $k=1$ in (2.5). In comparison, the shear modulus of the bistable energy $W^{*}$ for small deformations is $0.1$. This is consistent with measured values of the shear modulus of 3 mg/mL collagen in the range 80–100 Pa [74].

The values reported in Table S1 were used in all simulations of our own experiments. In simulations shown in Figs 4-6 corresponding to our own experiments, we used the experimentally measured values of initial particle radii, initial center distances and particle contractions. These are listed in Table S2.

In Figs 1d-1f, simulations are not intended to match specific experimental data, but to qualitatively show that the basic tether/hair morphology persists for different parameter choices although its details are affected. We use the monostable energy; see Table S1. All other parameters are as indicated in the caption.