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

Through mechanical forces, biological cells remodel the surrounding collagen network, generating striking deformation patterns. Tethers—tracts of high densification and fibre 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 fibres 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.


Introduction
Biological cells remodel the extracellular matrix (ECM) either biochemically, producing and degrading collagen fibres, 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 fibre alignment, localized within tether-like bands joining distant cell clusters (figure 1a). Tethers were recently observed between individual cells [3,8] as well (figure 1b). In addition, radial hair-like emanations from each cluster of cells (figure 1a,c) 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 towards each other [3]. Cancer cells preferentially invade along densified regions of the ECM with high fibre alignment [10][11][12][13][14][15][16][17]. Specific patterns of ECM density and alignment have been identified as biomarkers for breast cancer [12,13,[15][16][17].
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 © 2021 The Author(s) Published by the Royal Society. All rights reserved.
original thickness [9,18]. Density can be three to five times higher within tethers than without. These deformations are not only severe, but also spatially localized within tethers [5], and along radial hairlike bands that issue from individual cell clusters [6] (figure 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 (electronic supplementary material, 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 localization of densification is also unexpected. Centimetre-scale tethers induced by multi-cell tissue explants have well defined boundaries, across which density is virtually discontinuous [5,6] (figure 1a). Numerous thinner radial bands issue from each explant. At scales comparable to fibre length, density localization is gradual but still pronounced [3], whereas fewer radial bands emerge from each cell (figure 1b). Solitary cell clusters also produce radial bands (figure 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 localized, inhomogeneous deformations, not only at the scale of fibre network inhomogeneity [3] but at macroscopic scales as well [6]. Current nonlinear continuum models do not seem to explain these patterns of localization [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 fibres. Using active particles, our experiments establish that these severe, localized 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 [11,14,23].
The connection with instability was first established through experiments and modelling in fibre networks and open-cell foams [24] and later for fibrin [25]. Rather like rubber bands, individual fibres support tensile forces and may stiffen with increasing tension [26,27], but buckle under compression, losing stiffness and eventually collapsing. Larger polyhedral groups of fibres 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 fibres are mostly buckled and compacted, alternating with regions of normal density and low compressive strain, where fibres are largely unbent and loosely arranged. These two region types (high-and low-density) are separated by  table S1 for parameters. (a) Tether between millimetre sized explants (reproduced from [6]). (b) Tether between micrometre-sized individual cells (reproduced from [3]). (c) Radial hairs issuing from an isolated explant (reproduced from [6]). (d ) Effect of higher gradient energy, equation (2.4), in controlling number and fineness of hairlike radial bands. (i) ε = 0, (ii) ε = 0.01R, (iii) ε = 0.05R, where R is particle radius. (e) Shape of deformable explants and relative motion (i) k = 100 with fixed centres. (ii), k = 1 with fixed centres. (iii), k = 1 but centres are free, and move from black to green dot. ( f ) Fineness of phase mixture in simulations depends on numerical mesh resolution (increasing from left to right). Colour bar: ratio of deformed to undeformed density.
royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 18: 20200823 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 fibre, 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 localized densification and alignment that feature in ECM remodelling and serve as biomarkers for tumour invasion in fibrous ECM.

Modelling the energy landscape of the fibrous extracellular matrix
We model the nonlinear elastic energy landscape of fibrous collagen ECM, starting from the relation S = S(λ) between force S and effective stretch λ = d/l 0 for a single flexible elastic fibre; here d is the distance between endpoints, l 0 is the relaxed (reference) length. See Methods: Model and electronic supplementary material, figure S1a for details. We assume S(λ) stiffens in tension (its slope increases with λ > 1) as is common for biopolymers [26][27][28]. Aside from direct observation of buckling [29][30][31], very little is known about the post-buckling behaviour of individual collagen or fibrin fibres in compression [32], as λ decreases from 1 toward 0 (total collapse). This depends on their bending behaviour, which is difficult to characterize, because of their inhomogeneous, hierarchical structure [28] as bundles of loosely connected fibrils. We choose S(λ) as in figure 2a so that there is stiffening in tension (λ > 1) but loss of stiffness in compression (λ < 1) because of buckling. An example is S(λ) = μ(λ 3 − 1) (figure 2a). The corresponding elastic energy of the fibre is then w network has a uniform distribution of fibre orientations. When subjected to a 2D deformation, the stretch λ = λ(θ) of a fibre making an angle θ with the principal axes of stretch in the reference configuration is l(u) ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi (l 1 cos u) 2 þ (l 2 sin u) 2 p , where λ 1 , λ 2 are the principal stretches of the deformation (electronic supplementary material, 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 see Methods: Model for more details and an explicit form equation (5.2) of the 2D energy.

Instability in uniaxial compression
Surprisingly, the 2D energy densityŴ has an instability, evident in figure 2b, the stress-strain relation in uniaxial compression S 1 (l 1 ) ¼ @ W(l 1 , 1)=@l 1 . This relation is not monotonic, but its slope becomes negative for λ 1 < 0.45 (figure 2b). This is unexpected, since the original single fibre stress S(λ) is monotone increasing. Nonetheless, with increasing compression as λ 1 decreases towards 0, the network loses stiffness because of buckling of more fibres. Also, fibres 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 electronic supplementary material, Stiffness loss.
Physically, we expect that increasing compression (λ 1 → 0) leads to densification, and fibres getting squeezed together, thereby resisting further compression and eventually restoring stability. To account for resistance to crushing (due to fibre volume), we add a volume penalty term to the energy  Figure 2. One well or two? Designing an energy density for fibrous ECM. (a) Force-stretch curve of a single fibre that loses stiffness in compression and stiffens in tension. (b) Stress-stretch curve for uniaxial compression of 2D (orientation-averaged) energy density W from (2.1) (corresponding to (a)) has a decreasing unstable branch. Here S 1 (λ 1 ) = ∂W(λ 1 , 1)/∂λ 1 . (c) As in (b) but with added energy penalty (2.2) due to fibre volume, which resists extreme compression. At a compressive stress S 0 (red line), there are multiple stretch states, two stable (red dots) and one unstable (blue dot). (d) Adding an attraction potential (blue) to the fibre potential (green) corresponding to (a) produces a two-well bistable potential (red). (e) Force-stretch curve corresponding to the red potential in (d). (f ) Uniaxial compression curve of the corresponding orientation-averaged energy density W* (with volume penalty) shows bistability, with two stable stress-free states (red dots), and an unstable one (blue dot) in between.
royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 18: 20200823 (2.1) that resists collapse to zero volume in the crushing limit as the volume ratio J = λ 1 λ 2 → 0: where the penalty function Φ(J ) increases rapidly as the volume ratio J = λ 1 λ 2 → 0 as shown in electronic supplementary material, figure S1b. We use the specific form given by equation (5.3). This restores positive stiffness at extreme compression as λ approaches zero (figure 2c). Thus the decreasing branch with negative slope in the S 1 (λ 1 ) curve in figure 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 figure 2c and electronic supplementary material, figure S2b), there are three corresponding stretch states (red and blue dots in figure 2c, S2b). The middle one is unstable. The states λ d and λ u (red dots) are stable. These two states of compressive stretch λ d (densified) and λ 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,18,29,38,39]. Various factors contribute to this, e.g. adhesion and new crosslinking of fibres coming into contact due to densification, subject to van der Waals, or other noncovalent types of attraction [38,40]. To model this, we add a short range attractive adhesion potential w a (λ) (blue curve in figure 2d ) to the single fibre buckling potential w(λ) (green curve in figure 2d ). This represents minimal modelling of new crosslinks developing during fibre 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 fibre proximity [43]. This is modelled here by the shortrange adhesion potential w a (λ). It renders fibre collapse energetically favourable at short distances of the endpoints (λ near zero). The resulting single fibre potential w Ã (l) ¼ w(l) þ w a (l) is now a two-well potential (red curve in figure 2d ) with a new minimum corresponding to a collapsed fibre state (λ = 0) that is stable under zero load. Replacing w(λ) by w Ã (l) in equation (

Non-convexity in two dimensions
Next we consider the energy landscape in two dimensions. What sets both energy densities W and W* apart from previous continuum models of fibrous ECM [18,20,21] is nonconvexity (electronic supplementary material, figures 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 fibre 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 [44,48,49] at some strains, whereby some higher dimensional measure of stiffness becomes negative [50]. See electronic supplementary material, Ellipticity loss, also electronic supplementary material, figures S3d, S3e. The bistable energy density W* is a multi-well potential (electronic supplementary material, figures S2a, S3a). This allows the densified phase to be stable at zero stress, as has been observed in uncrosslinked collagen [8,18,29,38,39]. Although the monostable energy density W has a single minimum (electronic supplementary material, figures S2b, S3b), the associated Gibbs free energy W(F) − S 0 · F is a doublewell potential for a special value of the stress S 0 , with minima similar to those of W* (electronic supplementary material, figure 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 (electronic supplementary material, table S1), the minima (energy wells) of the 2D energy density W*(λ 1 , λ 2 ) occur at the undeformed state (λ 1 , λ 2 ) = (1, 1), and at the densified state (l Ã 1 , l Ã 2 ) = (0.2, 1.06). The latter is a severe compression with strain ε 1 = l Ã 1 − 1 = −0.8, combined with a small extension ε 2 = l Ã 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 fivefold increase in density with ratio 1/(l Ã 1 l Ã 2 ) ≈ 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 electronic supplementary material, Compatibility. This means that bands of the densified state can coexist in equilibrium, side by side within the undensified state. An additional minimum at (λ 1 , λ 2 ) = (0.45, 0.45) is not compatible with the undeformed state (1, 1); see electronic supplementary material, Compatibility. This explains why it is never encountered in our simulations.

Additional energy contributions
The elastic energy of a deformation y(x) is where Ω is the undeformed region occupied by the ECM. Going from the discrete energy of a random fibre network of characteristic fibre size ε to the continuum energy in an royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 18: 20200823 asymptotic expansion [34,51] as ε → 0, one obtains equation (2.3) as the first term, followed by a higher gradient term quadratic in 1rry(x). We choose a simple form of isotropic higher gradient energy [52,53] To model contractile multi-cell clusters (explants [5,6], acini [9,18]), we let Ω contain initially circular cavities. At the boundary of each cavity, forces can be exerted onto the ECM. In contrast to previous studies [3,18,20,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 (electronic supplementary material, figure S1c). This contributes to the energy a term where Γ i is the boundary of the ith cavity (i = 1, …, n), k is the spring stiffness, z i the variable centre position of the ith 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: eggshaped and pointed towards the tether (observation figure 1a reproduced from [6]; simulations figure 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, equation (5.5), electronic supplementary material, figure S1d, and electronic supplementary material, Simulation parameters. The total energy to be minimized is the sum of equations (2.3), (2.4) and (2.5):

Experiments and simulations of extracellular matrix 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 fibres) 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 fibres) in the densified zones until the latter have fully developed. That the densification is associated with compressive strain can be directly deduced from [9, electronic supplementary material, video movie_S14], 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 localization 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°C (Methods: Experiments with active particles). Advantages of this are that we can control the amount of contraction via temperature, and even cause reverse reexpansion 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 minimize the total energy, equation (2.6), with respect to the deformation field y(x) and the particle/cluster centre positions 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 equation (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 equation (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 (figure 1a-c) and encountered in our simulations for a wide range of model parameters (figure 1d-f ). 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 electronic supplementary material, Compatibility). This is consistent with experiments [18, figs S2, S5] and is observed in our simulations (figure 3a,b).
How close do two active particles have to be (figure 4), and how much must they contract in order to form a tether royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 18: 20200823 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 figure 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 figure 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,18,20,54,55], we allow clusters to change shape and move during ECM deformation in our model. See equation (2.5), also Methods: Model for active particles and electronic supplementary material, figure S1d. In experiments [5,6], isolated explants remain roughly circular (figure 1c) after contraction. By contrast, neighbouring clusters connected by a tether lose circularity [5,6] and become egg-shaped with pointed ends toward each other (figure 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 equation (2.5), our model predicts similar shape changes with pointed end where the tether makes contact with the cluster (compare panels (ii) and (iii) in figure 1e to shape of explants in figure 1a). In order for simulated explants to remain circular after contraction as typically assumed, we must choose an artificially high stiffness (k = 100 in equation (2.5) Figure 3. Principal stretches and fibre distributions: principal stretches: (a) λ 1 < 1 (compressive) and (b) λ 2 > 1 (tensile) from the simulation in figure 6e. Note sharp change in λ 1 across tether boundary indicating discontinuity of stretch normal to phase boundary, in contrast to gradual change in λ 2 (tangential stretch) as predicted by compatibility; (c) and (d ) same but for expanding particle, figure 5e. Scale bar: principal stretch. (e-h) Deformed fibre distribution at point where λ 1 < 1 and λ 2 > 1, with peak along largest-stretch direction. (e) Predicted distribution at a point within tether in (a,b) (λ 1 , λ 2 ) = (0.2, 1.06) with peak along the tether axis, normal to contractile-particle boundary. (f ) Measured distribution from [15] outside contractile tumour spheroid, where angles 0 and 90 are along and normal to tumour boundary, reproduced from fig. 5 in [15]. (g) Predicted distribution outside simulated expanded particle in (c,d ), also figure 5f , with (λ 1 , λ 2 ) = (0.113, 1.5) and peak at direction parallel to particle boundary. (h) Measured distribution from [15] outside growing tumour spheroid with peak in direction parallel to tumour spheroid boundary, reproduced from fig. 5 in [15].
royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 18: 20200823 Contracting active particles in our experiments (figure 4a, d) and simulations (figure 4b,e) move towards each other by as much as 10% of their original distance (figure 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 (figure 4a). Other tethers are in the form of uniform bands making full contact with particles (figure 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 zerostrain austenite [45,46] (see electronic supplementary material, 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 figure 4c. Here also, energy minimization compels strains not to stray far from energy-density minima. The active particles in figure distance in mean deformed particle radii % mean decrease in particle radius  (d) Two smaller particles at higher 50% contraction with full-contact tether (except torn fibres at right particle) and partially enveloping densification. (e) Simulation of (d ) shows full-contact tether with no splitting. ( f ) Hypothetical simulation of particles of (a) but extremely high contraction of 80% (not attained by active particles) shows uniform tether and densification enveloping particles, due to compatibility of contraction with densified phase. (g) Predicting whether a tether forms between two particles. Blue curve: separatrix constructed from simulations. Axes: % decrease in particle radius versus deformed distance (in deformed particle radii). Above separatrix, tethers are predicted to form between particle pairs. No tether is predicted to form below separatrix. Our experimental data (each particle pair is one point) abide by the prediction: blue points: tether has formed. Red points: no tether has formed. Colour bar: ratio of deformed to undeformed density.
royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 18: 20200823 4a contracted by u 0 /R = 38%. The azimuthal stretch l u ¼ 1 À u 0 =R ¼ 0:62 imposed at the particle boundary by contraction is incompatible with the stretch l Ã 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 minimize contact with the particle boundary (experiment: Fig. 4a, simulation 4b).
The particles of figure 4d contracted more (52%). The resulting tether does not split, but is in full contact with the left particle (torn fibres near the right particle keep the tether from extending to the left particle boundary). A corresponding simulation (figure 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 l u % l Ã 1 (is compatible with the energy minimal stretch l Ã 1 ). In hypothetical simulations of this contraction level, the densified phase is in full contact and envelopes the particle because of enhanced compatibility (figure 4f ). These observations support the hypothesis that tether morphology is decided by strain compatibility as in other types of coherent phase transitions (martensitic).

Contracting versus expanding particles
What causes the densified phase to split into thin radial bands around particles, such as 'hairs' emanating from clusters in figure 1c (reproduced from [5,6]) and in figure 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 (figure 5b). Instead, an expanding particle would create radial compression, hence a circumferential phase boundary and densified ring enveloping the particle (figure 5e). The stark contrast between the ECM's response to contracting and expanding particles is seen in our experiments (contracting particle figure 5b versus expanding figure 5e) and for the first time captured by a continuum model simulation (contracting figure 5c versus expanding particle figure 5f ). See also the densified ring surrounding a re-expanded particle in experiment figure 6c, simulation figure 6f . Remarkably, this explains the sharply defined circumferential layer of densification typically observed surrounding expanding tumour spheroids; see, e.g. [14, fig. 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 (figures 1 and 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 Figure 5. Contracting versus expanding particles. Experiments (green) and simulations (purple) matching initial geometry of particles and contraction level. (a-c) Contracting particle. (a) Undeformed PNIPAAm particle. (b) PNIPAAm particle radially contracted by 50%. (c) Simulation of contracting particle in (b). Note radial 'hairs' and radial symmetry breaking. (d-f ) Expanding particle. (d) Unexpanded (undeformed) PNIPAAm particle. (e) PNIPAAm particle radially expanded by 50%. Note bright circumferential densification layer. (f ) Simulation prediction for radially expanded particle in (e). Note circumferential densified phase layer surrounding particle, circular phase boundary and radial symmetry maintained. Colour bar: ratio of deformed to undeformed density.
royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 18: 20200823 fine phase mixtures establish a further connection with the nonlinear elastic theory of phase transitions in crystals [45,46]. In multi-well energy minimization 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 (figure 1f ).
To show that this finite-element mesh dependence is not an artefact of the numerical method, we add the highergradient term equation (2.4) to the energy. This introduces a length-scale proportional to ε, an additional material parameter related to characteristic fibre length, bending stiffness and other fibre network parameters [34]. For larger ε, there are fewer, thicker radial bands (figure 1d). Above a critical value of ε, they disappear, but the tether persists. The presence of this term has a smoothening effect as it penalizes high strain gradients; strain discontinuities are replaced by transition layers with thickness of order ε [52,53]. This eliminates mesh dependence at numerical mesh sizes below ε, which controls the scale of the phase mixture. In practice we expect ε to be of the order of the free fibre length in the network. Accordingly, many hairlike bands are experimentally observed (figure 1a,c) issuing from millimetre-size explants [6], (where ε is very small compared to explant size), while only a few (figure 1b) from single micrometre-scale fibroblasts [3].

Residual densified zones and inelasticity
Do the densification patterns in the ECM persist after cellexerted contractile forces cease [8,18]? In our experiments, active particles (figure 6a) contracted upon heating to 39°C, causing densified tethers and radial hairs to appear (figure 6b), then expanded to their original radius upon cooling back to 26°C (figure 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 (figure 6c), but they were thinner and less prominent than ones appearing during particle contraction (figure 6b). Many radial bands issuing from particles largely disappeared. We performed contraction-expansion cycle simulations (figure 6d-f ) 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 (figure 6a-c) confirmed by our model simulations (figure 6d-f ) 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 figure 6f to experimental figure 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 localized 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,18,29,38,39]. 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 non-convex 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 figure 6d-f .

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 figure 7c, reproduced from [9]). The original tether disappears; instead tethers form that bypass the crack by going around its corners (arrows in figure 7c). Our simulation (figure 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. E.g. [9, fig. 1D). A simulation of our model with a hypothetical distribution of contractile acini is shown in figure 7. In figure 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.

Discussion
The continuum model we have proposed differs from previous ones, because it is based on a nonlinear elastic, nonconvex strain energy density, which 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 fibres. 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 fibre 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 fibres [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 non-affinity in network displacements due to factors such as network randomness, usually restricted to small deformations [65]. The present model averages over this type of non-affinity by coarse-graining from a discrete network to a continuum description, but identifies a more severe type of non-affine deformations: strain localization 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 royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 18: 20200823 fibre 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 minimizing 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 (figures 4-6) and the change in distance between particles, figure 6g, caused by deformation from our experiments. On the other hand, the model only involves five constitutive parameters (electronic supplementary material, 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? Figure 7. Simulation with multiple acini. Their contraction causes an extensive network of tethers to form which mechanically remodels the ECM (a) deformed configuration. White: contracted acini, red: densified phase, blue: undensified phase. (b) Superposed undeformed configuration (dark blue acini) prior to contraction and deformed configuration (white acini) showing motion and shape distortion as well as contraction. Orange: densified phase, light-blue: undensified phase in the deformed configuration. (c-d ) A crack (dotted line) is cut across a tether between acini [9], the tether disappears. (c) Two new tethers turn around the corners of the crack to bypass it (reproduced from [9]). (d ) Simulation of pair of acini with crack predicts tether bypassing the crack. Colour bar: ratio of deformed to undeformed density.
royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 18: 20200823 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 (figure 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 [17,66] and fibre alignment [12,13,15,16]. Both arise in our model and experiments, as a result of phase change brought about by buckling fibre collapse. The densified phase in our model and experiments, also in [9,18], is characterized by a three-to fivefold increase in density, but also strong fibre alignment toward the direction of maximum stretch, as compression squeezes an originally uniform angular distribution of fibres away from the compressive direction and toward the orthogonal tensile direction (see electronic supplementary material, Fiber alignment in the densified phase). Hence, within a tether formed between contracting particles (figure 6e) the fibre angle distribution we predict, electronic supplementary material, equation (S2), is as in figure 3e, with a peak along the tether axis, normal to the contracting spheroid boundary. This compares well with figure 3f (reproduced from [15]) which shows a typical fibre distribution associated with TACS-3 [12,13,15,16] or tumourassociated collagen signature 3 (aligned fibres normal to contractile tumour spheroid boundary) a biomarker for invasion of tumour cells along tracts of aligned fibres. 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 figures 5f and 3c,d of expanding particles (figure 5c) predicts a fibre distribution figure 3g peaked in the direction tangent to the (expanding) boundary. This agrees qualitatively with collagen signature TACS-2 [12,13,15,16] (fibres compacted and aligned parallel to a growing tumour boundary) with measured distribution figure 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 [13,15]) or tumour cell contraction (TACS-3) [16].
Other than serving as a track for cell migration, the aligned and densified collagen fibres 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,[59][60][61][62][63][64]. Remarkably, we find that such non-convex 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, patternforming 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.

Experiments with active particles
5.1.1. PNIPAAm particle generation PNIPAAm particles were generated using two different recipes, one giving larger particles (diameter 150 μm) and the other giving smaller particles (diameter 75 μ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 localized 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 h. It was then maintained under nitrogen for 10 min before stirring at 350 r.p.m. 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 × tris-buffered saline was prepared. Ten microlitres of 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 r.p.m. instead of 350 r.p.m. A solution of 1 g Nisopolyacrylamide (Sigma), 3.75 ml of 2% bis-acrylamide solution (Bio-Rad), 0.05 g ammonium persulfate (Bio-Rad), 1.5 ml of 1 × tris-buffered saline and 3.75 ml deionized water was prepared. Again, 10 μ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 deionized water. The particles were then filtered with a cell strainer to keep particles of diameter > 40 μm. The particles were resuspended in 1 × PBS.
The PNIPAAm particles were embedded in networks of rat tail collagen I (Corning), which was fluorescently labelled 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 neutralized using 100 mM HEPES buffer to a final concentration of 3 mg ml −1 and mixed with the PNIPAAm particles. The collagen then polymerized at 27°C for 1 h. After polymerization, the collagen formed a network of fibres that branch into multiple fibres at discrete nodes. The average distance along a fibre between nodes, which we refer to as the fibre length, was measured as in our previous work [39]. Briefly, images of the fibres were segmented [69], giving segmented areas of the void space between fibres. The average of the segmented areas was 3.9 μm 2 . Hence, the mesh size, which is proportional the square root of the segmented areas, was 1.97 μm. Next, we accounted for the fact that some fibres in the image cross without connecting [70], by selecting representative fibres and manually counting the number of unconnected crossings. The product of the mesh size and number of unconnected crossings gave the fibre length, of 19.7 ± 2 μm (mean ± s.d.).

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 × 0.75 NA air objective (Nikon) using a Zyla sCMOS camera (Andor). The particles diameters and the distances between particle centres 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).

Energy density
For a single fibre, we introduce the effective stretch λ, which equals the distance between its endpoints divided by its undeformed, or relaxed, length. The energy of a single fibre can be written as w(λ) as a function of effective stretch λ. When the fibre is in tension, it is straight and λ equals the actual stretch (strain +1), while w(λ) equals the elastic energy due to stretching of the fibre. When it is in compression, it may be buckled, in which case the elastic bending energy of the fibre can still be expressed as a function w(λ) of the distance between its endpoints, hence of the effective stretch λ. See electronic supplementary material, figure S1a. In order to model a 1D two well energy w(λ) for a single fibre like the red curve in figure 1d, we start with the derivative S(λ) = d w(λ)/dλ, which represents force as a function of stretch. We choose a polynomial that vanishes at 0, 1 and an intermediate value, to qualitatively represent the curve in figure 2e. Here a m is a parameter that controls the relative height of the two wells (minima) of w(λ). Integrating this with respect to λ gives an energy We model the ECM as a 2D nonlinear elastic continuum undergoing possibly large deformations y(x) where a particle with position vector x in the undeformed state is mapped to deformed position y = y(x). The elastic strain energy density of the material can be written as a function W(F) of the deformation gradient F ¼ ry, which is a 2 × 2 matrix. We model the random fibre network as an isotropic material, which means that W depends on F only through the principal stretches λ 1 , λ 2 , whose squares are the eigenvalues of the Cauchy-Green deformation matrix F T F. Equivalently W is a function of the two deformation invariants Here, J is the Jacobian determinant of the deformation and the ratio of deformed density ρ to undeformed density ρ 0 satisfies ρ/ρ 0 = 1/J.
To connect the single fibre energy with the 2D strain energy density W, we follow [34,35]. We assume the network contains a uniform distribution of fibres. A homogeneous deformation (with constant strain) is equivalent to a biaxial stretch in two orthogonal directions with stretches λ 1 , λ 2 . A fibre of undeformed length l that makes an angle θ with the principal stretch axes in the undeformed state, will have endpoints at (0, 0) and (l cos θ, l sin θ). After deformation the latter will become (λ 1 l cos θ, λ 2 l sin θ). As a result, the stretch ratio of the fibre will be l(u) ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi (l 1 cos u) 2 þ (l 2 sin u) 2 p , and its energy will be w(λ(θ)). Summing over all fibre orientation angles θ, we find the elastic energy density of the network which gives equation (2.1). It turns out that for w(λ) given by equation (5.1), the integral in equation (2.1) can be evaluated explicitlỹ Here, I 1 ¼ trFTF ¼ l 2 1 þ l 2 2 and J ¼ det F ¼ l 1 l 2 are the 2D deformation invariants. We then add a fibre volume penalty term to the energy to account for resistance of densified fibres to complete crushing by virtue of their non-zero volume. This term increases the energy abruptly when the Jacobian (volume ratio) J ¼ det 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

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 equation (2.5); see also electronic supplementary material, figure 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.
C{y; z 1 , . . . , z n } ¼ where Γ i is the boundary of the ith cavity (i = 1, …, n), k is the spring stiffness, z i the variable (deformed) centre position of the ith cluster, 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 electronic supplementary material, figure S1d.

Numerical method
We employ the finite-element method, based on a triangulation of the domain Ω. 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 minimization problem on the finiteelement space. In order to capture areas of high densification royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 18: 20200823 accurately a local mesh refinement strategy close to the areas of phase transition is adopted. Challenges include the subtle nonlinear 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, ε > 0 in equation (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 regularization by higher gradients. We thus use the same discrete spaces in all models considered in our study (ε = 0 and ε > 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 y h denote a function of the discrete finiteelement space of piecewise polynomials (C 0 -conforming), then the discretized functional for has the form: 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 ⊗ denotes the tensor product.
The discretization has been implemented in FEniCs [71]. For the minimization of the discrete energy functional a parallelized 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 minimized using a parallelized version of the nonlinear conjugate gradient method [72], a successful iterative method for large-scale nonlinear optimization problems. When ellipticity of the corresponding Euler-Lagrange equation for the unregularized problem holds, Newton's method is an efficient nonlinear minimization 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 regularized energy functionals.