Modeling epithelial tissues as active-elastic sheets reproduce contraction pulses and predict rip resistance

Confluent epithelial tissues can be viewed as soft active solids, as their individual cells contract in response to local conditions. Little is known about the emergent properties of such materials. Empirical observations have shown contraction waves propagation in various epithelia, yet the governing mechanism, as well as its physiological function, is still unclear. Here we propose an experiment-inspired model for such dynamic epithelia. We show how the widespread cellular response of contraction-under-tension is sufficient to give rise to propagating contraction pulses, by mapping numerically and theoretically the consequences of such a cellular response. The model explains observed phenomena but also predicts enhanced rip-resistance as an emergent property of such cellular sheets. Unlike healing post-rupture, these sheets avoid it by actively re-distributing external stresses across their surface. The mechanism is relevant to a broad class of tissues, especially such under challenging mechanical conditions, and may inspire engineering of synthetic materials. Observations on confluent epithelial tissues show the emergence of dynamic contraction patterns that are suspected to be governed mechanically. Here, the authors propose a model for epithelial sheets and show that cells’ Extension-Induced-Contraction response explains experimentally-observed contraction pulses, that along with a cell softening response enhances epithelial resistance to rupture.

pithelial tissues are confluent cell sheets, made of individual active units (the cells) that are mechanically coupled. Such living materials far from equilibrium have been suggested to exhibit unexpected material properties, ranging from solid-fluid transitions [1][2][3] , active super-elasticity 4 , and exotic forms of viscoelasticity 5,6 . These active cellular sheets act as chemical barriers and apply mechanical protection and shape to the underlying organ. As such, they are required to keep their stability and integrity at all times. Evidence shows that when put under tension, these substrates can quickly and reversibly rigidify [7][8][9][10][11][12][13] or soften 4,10,12,14,15 (aside slower remodeling via plastic deformations like oriented divisions or cell intercalations 16,17 ). The interplay between these responses, orchestrated to keep integrity across millions of cells and various mechanical threats, is still unclear.
In recent years observations are accumulating showing dynamic contraction patterns in confluent epithelia that are suspected to be governed mechanically. These evidence include contraction fronts in drosophila embryo during development 18 and ongoing density waves seen in migrating cell monolayers in vitro, either in confinement 11,19 , during expansion 11,20 or in response to substrate-shear 21 . Theoretical works suggested models to explain the contractile patterns [22][23][24][25] , models that suggest an interplay between mechanical, chemical, and geometrical fields, diffusion or active transport, and cell polarity. A few works 11,18,25 suggested that a cellular response of contractionunder-tension may be the only requirement to create contraction waves that propagate mechanically. Both the mentioned observations and the theoretical models exhibit time scales in the range of minutes to hours.
Recently we reported the observation of ultrafast contraction dynamics in the thin, suspended epithelium of the marine animal T. adhaerens 26 . These traveling pulses initiate spontaneously, irregularly, propagate at high speed (1-3 cells or 10-30 μm in 1 s), and transmit a fast contraction event (50% cell area in 1 s) across the tissue (Fig. 1a-d, Supplementary Note 1). During these dynamics, and in time scales of up to 10 min, no cell rearrangements (neighbor exchange) have been seen in the epithelium 26 . The contraction pulses can propagate radially or uniaxially, split at the propagating fronts, and annihilate each other. Importantly, this early-divergent animal has no reported muscles, neurons or synapses, and its epithelium has no gap junctions that can support chemical transport between cells. In addition, the contraction propagation speeds exclude too slow biological processes, such as transcription/translation, from being involved (Supplementary Note 1). All these raise the speculation that mechanics governs the contraction propagation. The pulses are seen in the suspended epithelium, that has no known ECM, while the animal is freely moving. Due to its erratic, locally-driven ciliary locomotion, the tissue is found constantly under alternating tensile/ compression stresses. Despite these external stresses, and the intrinsic contractions, the tissue does not rupture 26 .
Extension-induced contraction (EIC) is a common, well known cellular response in epithelia. EIC has been demonstrated experimentally in various cells, and different molecular mechanisms have been shown to participate in it (e.g., mechanosensitive calcium channels, ERK protein activation, actin alignment, myosin recruitment, conformational changes in adherens junction, and more) 5,[7][8][9][10]25,[27][28][29][30] . The diversity of mechanisms for EIC generation hints towards its crucial evolutionary advantages. In addition, in many systems, it was shown that aside EIC, cells also respond to stretch by softening/yielding 4,10,12,14,15 . The two antagonist responses, of contraction/stiffening and yielding/softening, have been shown to coexist in cytoskeletal assays 7,14 , and were demonstrated, separately, in the same cell type (MDCK cells 4,11 ), but their exact interplay in live cells and tissues is still unknown.
Here we suggest a model for dynamic contraction patterns emerging from the cellular response of EIC. We frame our model regardless of the specifics of the intra-cellular process, and show that cellular EIC may be the only requirement for contraction propagation and the observed dynamics. This is in agreement with previously suggested models 11,18,25 . Built from direct, high spatial-temporal measurements in T. adhaerens, and using numerical and analytical considerations, we construct our model which is very general and purely mechanical. Its basic feature is longitudinal density waves that propagate utilizing local contractile activity, hence can be described as "active-acoustics". These nonlinear pulses propagate without dissipation, hence can be viewed also as mechanical "trigger waves". Primarily, we write scaling laws for these contraction pulses and quantitatively map the different pulsatile modes in parameter space. Our results may explain many phenomena seen in T. adhaerens and other epithelia, such as single or repeated contractile pulses, either oscillatory or in response to external stress, and pulse annihilation. Finally, the model predicts another emergent outcome of EIC in such cellular sheets: enhanced resistance to rupture under tension. The contraction dynamics increases tissue cohesion by distributing strains homogenously throughout the surface. By numerically introducing a cell softening response to the model, on top of EIC, we see further enhancement of cohesion in the cellular sheet, via elimination of high-stress values at cell-cell junctions. Together, the two cellular responses prevent tissue failure under tension in both cells and junctions simultaneously in a process we name "active cohesion".
As both EIC and yielding are common epithelial responses, and as tissue integrity is at the heart of any epithelium function, "active cohesion" may be relevant in a wide range of epithelia, during development, physiology, pathology, and regeneration, even if manifested in different time and length scales (Supplementary Note 1c). In addition, our work may contribute to the understanding of physical principles in active solids, and to the engineering of tissues and 'smart' materials.

Results
Single-cell dynamics. We begin by looking at the recentlymeasured contraction profile in T. adhaerens 26 (Fig. 1c). The profile shows the strain evolution of a single cell during a single contraction event in vivo, as averaged from multiple cellular events that did not propagate as pulses in the tissue. On average, during such a contraction event, the cell area increases gradually to a critical point, at~110% of its initial value, at which it abruptly decreases to~50%. Finally, the cell relaxes slowly towards its initial area. Despite the overdamped conditions, the contraction reaches significantly below the cell's steady-state size and beyond the~10% level of expansion ("overshoot"). This shows that the contraction is an active process that consumes internal cell energy, and explains why it involves observed force and time scales that are different from the viscoelastic ones. As the contraction seem to happen at a critical cell size, an additional length scale is required, to describe the extension needed for activation. These three scales-the active force, its duration, and the critical cell size-are the minimal requirements for the model.
We, therefore, model a single cell as an overdamped elastic entity (a spring with rest length l 0 and elastic modulus k, in a media with viscosity γ) that is connected in parallel to an active contractile unit (Fig. 1e). The active unit implements a simple EIC, using the three scales discussed above: when the cell reaches a critical length l c , it immediately applies a fixed compression force f c for a duration of t c (see "Discussion" for further reasoning and alternative options, see methods for implementation). Inspired by T. adhaerens data, we take the active forces to be larger than the elastic force at criticality (f c > kl c ) and the active time scale shorter than the passive one (t c < k/γ). Finally, we assume the spring to be linearly elastic. Effects of cell-volume conservation, though shown irrelevant to T. adhaerens 26 , may require additional nonlinearities.
The mechanical circuit proposed in Fig. 1e can be found in one of two modes: an excitable mode, i.e., contract only in response to external stretch (that is when l c > l 0 ), or in an oscillatory mode, i.e., oscillate spontaneously without any external stimulation, as a relaxation oscillator (that is when l c < l 0 ) (Fig. 1f). Evidence for such an oscillatory cell behavior have been reported in T. adhaerens (Fig. 1d), and in MDCK monolayers 11,[19][20][21] . These isolated-cell behaviors can be mathematically described in a piecewise manner (Supplementary Note 2).
We will now examine the emergent behavior of a finite 1D chain of such cells (Fig. 1g) or a 2D sheet of cells with similar properties (Fig. 1h).
Dynamics in one dimension. In the 1D case, we take each cell to be a spring connected in parallel to the active-contractile unit. All cells experience viscous drag from the media. Cell-cell junctions are the nodes, that feel tension due to forces from the nearest cells. Throughout our investigation, we choose a finite yet large number of cells, N, and free boundaries, in order to imitate T. adhaerens, (details in Supplementary Note 3).
First, we consider a one-dimensional chain of N identical excitable cells (i.e. l c > l 0 ). Initiating all cells at length l 0 , the chain will remain at rest. We now introduce an initial stretch, by initiating a rim-cell at l i = l c + ε, after which it is set free. The dynamics starts as the stretched cell actively contracts, stretching the next cell in line, due to the spatial coupling in the overdamped conditions. If the induced stretch in the neighbor cell is sufficient, it triggers active contraction, that in turn activates the next cell and so forth. A contraction pulse then propagates throughout the tissue at a constant speed, V (Supplementary Movie 1, Fig. 2a, b). The emerging longitudinal perturbation is not an acoustic (inertial) wave, but a slower, trigger-wave in an excitable media, that consumes cellular energy.
The symmetric shape of the strain signal (Fig. 2b, c) may be surprising, as we defined an asymmetric EIC profile for the individual cell: the expansion induces a much larger contraction. In fact, a cell participating in a contraction pulse shows a symmetric strain profile, with a delayed, reduced shrinkage. To understand this emerging effect, let us examine the pulse propagation at the bulk of this large, deeply overdamped tissue. We notice that essentially the only nodes moving are the ones at the interface between active and passive cells-all other nodes are at rest due to force balance. Therefore, at the interface, the sum of the active and passive cell lengths is fixed. Thus, the dynamics in time of a bulk-cell is as follows: (i) starts at l 0 (ii) gets stretched by  26 , where more related data can be found). a Sequential snapshots from contraction pulses in the tissue, imaged using fluorescent membrane stain. The fluorophore intensity (greyscale) is related to membrane density hence to cell area: white areas are contracted cells and dark areas are relaxed, larger ones. White arrows indicate pulse propagation direction. Scale bar: 35 μm. b Sequential snapshots from a cellular contraction event. The contraction induces the expansion of neighboring cells. Scale bar: 3 μm. c The average contraction profile over 746 individual cells within the tissue (such that are not participating in pulses). The gray sections mark the expansion phase, followed by a quick contraction and a slower relaxation. The horizontal line represents the average cell area over a longer time (40 s). d A single event where oscillatory contraction activity is seen in a single cell (a cell that is not participating in a pulse). e The mechanical circuit representing the single cell in our model. A spring with rest length l 0 and stiffness k is connected in parallel to the active unit, which is operating the active response (represented in red arrows): when reaching above a critical length l c , a cell contracts with constant force f c for a given period of time t c . In the simulations we use a boxcar function to model immediate activation and termination of the force. The cell is found in a media with viscosity γ. f The behavior of such isolated cells, in either an excitable mode (l 0 < l c ), or oscillatory mode (l c < l 0 ). Dashed red lines represent asymptotic tendencies. Expressions for l min , l ∞ and the piecewise exponential solutions are in Supplementary Note 2. g-h Sketches of the 1D and 2D settings in our multi-cell simulations. The degrees of freedom are the vertices between cells that are free to move under cellular forces and viscous drag. The systems are large but finite and the boundary is free.
its active neighbor to exactly l c , and gets activated (iii) contracts until reaching l 0 , at which point its passive neighbor reaches exactly l c and gets activated (iv) our cell keeps applying contraction forces without changing its length, as both its neighbors are also active (v) only when an active neighbor deactivates at t c and starts to relax-our cell finally shrinks, until reaching t c itself. As a result, the emergent pulse in the tissue is composed of an extension front, followed by an identical (but opposite-sign) contraction front, and in between, all cells are actively contracting yet remain at their rest-length ( Fig. 2a- Hence, in order to derive a simple scaling for the pulse behavior, as a first approximation, we consider two cells with fixed boundaries (this is practically the case at the interface between active and passive cells). By calculating the time it takes for one cell, contracting with f c , to excite its neighbor t * ¼ , we estimate the propagation speed of the pulse V % f c l 0 γðl c Àl 0 Þ . The width of the pulse is the distance the pulse propagates in time t c W = Vt c . The pulse amplitude, which is the difference between maximal and minimal cell size, is found to be amp ≈ 2(l c −l 0 ). Our numerical results confirm all these scaling laws (Fig. 2h-j). In addition, these scaling laws fit previously measured data in epithelia and may explain both fast and slow pulse propagations: 100 μm/s in T. adhaerens and 2-3 μm/min in MDCK monolayers (Supplementary Note 1). We use the above scaling to find the requirements for pulses to appear: The initial stretch should excite the first cell (l i > l c ), and the impulse of active contraction should be large enough to excite the next cell. Using a set of three non-dimensional parametersnormalized timet t c k γ , normalized forcef f c kl c , and normalized strainε l c Àl 0 l 0 -this requirement estimates tot>ε f (Supplementary Note 3). When the system satisfies these criteria, a pulse propagates indefinitely in the tissue with a fixed speed. In the absence of these requirements, the initial stretch decays, and bulk cells stay at their rest lengths indefinitely. We notice, that changing the contraction parameters f c ,t c changes the pulse velocity V and width W, but not the amplitude, amp (Fig. 2d-e). Only by changing l c the amplitude is altered (Fig. 2f). As a result, for a given set of activation parameters (l c , f c , t c ), if the external stimulation lasts longer than t c (in our case, if the initial excitation, l i is large enough), it generates a "spike train"-a series of sequential pulses, all carrying the same quanta of strain-amp (Fig. 2g, Supplementary Note 3). This effect resembles neuronal action potential: sequential identical firing events in response to above-threshold stimulation 31,32 .
Next, we consider a row of N identical oscillatory cells (l c < l 0 ). All cells are initiated at l 0 and the boundaries are set to be free. Each cell is a relaxation oscillator, that would beat spontaneously   in isolation. However, now the cells are additionally subjected to forces coming from their neighbors. We show, that for a wide range of parameters, the system converges to exhibit ordered pulsations regardless of the initial condition. Initially, all cells contract in a transient irregular phase. Eventually, the system reaches a dynamic steady state (at characteristic time T ss ) where contraction pulses are initiated repeatedly and regularly at the edges and annihilate at the center (Fig. 3a, b, Supplementary Movie 1). Annihilation can be seen as a result of the shrinkage front-a "recovering" regime at the back of the pulse that makes the tissue harder to excite (again resembling the refractory period of a neuronal action potential). The overall tissue length L(t) shrinks from L(t 0 )=N*l 0 to a steady-state size L ss and oscillates around it with amplitude Amp and highest frequency Ω (Fig. 3c).
In this mode, cells in the bulk (far from both rim and annihilation points) are effectively fixed at a compressed size: due to viscosity and a large number of cells, the time scale for bulk cells relaxation is much longer than the typical interval between pulses. Therefore, the pulses propagate through a still and uniform background of cells, that are all at l~l ss < l c . As a result, a pulse propagates at a fixed speed. Rim cells are the least constrained and hence relax the fastest after a contraction. When they reach criticality, they initiate a new pulse. The pulse profile features are the same as in the excitable mode. The way the system parameters (l c , f c , t c , k, γ, N) relate to the emergent overall tissue measurables (L ss , T ss , V, W, Amp, Ω) is plotted from our numerical results in Fig. 3d, e, and more elaborately in Supplementary Fig. 1. Most measurables are invariant to the system size (except T ss and frequencies lower than Ω), hence are effective "material properties".
Other solutions exist in the oscillatory mode, aside the ordered pulsations, as shown in the phase diagram (Fig. 3f): When the contraction force is weaker than the elastic retraction at l c (i.e., f c < kl c ) the system will be "stuck", i.e. continuously apply compression forces but exhibit fixed cell size, that is above l c . These states may look "flickering" with constant local activation/ deactivation (Supplementary Movie 1). In contrast, when the active force is too high, we reach a non-physical regime of the model where cells collapse to negative size (in reality, nonlinear elasticity at the limit of compressibility will prevent that). When the viscosity is high, the time it takes to reach steady state is very long, hence in realistic time scales the dynamics may look irregular. Finally, when viscosity is low, irregular pulsations emerge, as a result of various elastic modes that are not overdamped (Supplementary Movie 1). The sharp transition between the regular and irregular traveling pulsations is depicted in Fig. 3g, as a transition from perfect limit cycles to "smeared" chaotic activity.
Dynamics in two dimensions. We now apply the same principles of EIC on 2D cellular sheets. We choose a hexagonal grid to describe the cells, while each node is an independent degree of freedom, and each edge represents a two-cell-junction. We use a vertex model, that assumes separate stiffness of a cell perimeter, k p , and a cell area, k a , and is a common modeling approach for confluent tissues 2,33,34 (See "Methods"). In the vast parameter space of 2D vertex models we focus on a regime that is the most relevant to the biological context of T. adhaerens (low area stiffness k a , high f c -capable of reducing cell area by 50%, ϵ c~1 0% Fig. 3 Contraction pulses in a 1D oscillatory tissue. a Sequential snapshots from a dynamic simulation (Supplementary Movie 1) of a 1D chain of oscillatory cells (in which the critical cell length l c , is shorter than its rest length l 0 ). Cells are shown at their location along the chain, x, greyscale represents their length and red asters represent activated cells. At t 0 all cells are set to their elastic rest length l 0 . At t 1 the system is during its transient shrinking phase. At t 2-5 the system is at its dynamic steady-state: contraction pulses are constantly propagating from the edges to the bulk of the tissue, where they collide and annihilate. b A kymograph representation of the dynamics shows the transient shrinking phase, converging to the dynamic steady-state, where pulses propagate at constant speed from the rims to their annihilating point at the center. c A sample measurement of the entire system length L(t). During the transient shrinkage, a kink may be seen. Then, further exponential decay (with time scale we label T ss ) brings the system to its asymptotic average length L ss , around which it oscillates with frequency Ω and amplitude Amp. d-e Time series of L(t) as a function of the contraction force f c and the viscosity γ. We show early and late intervals, taken from long simulation runs. While we focus on ordered pulsatile behavior, different possible types of dynamics are seen at the extremities of parameter space. These include all-contractile (no-deactivation) mode; cell collapse to negative area (non-physical); increasingly long T ss (long transient); and irregular traveling pulses (chaos). More details in Supplementary Fig. 1. f A phase diagram shows the different types of dynamics as a function of the non-dimensional parameters e f(related to force) andt (related to time). g A plot of the system length L vs. its variation in time ∂L=∂t shows the sharp transition from an ordered pulsation mode (limit cycles) to chaos as a function oft. The trajectory is depicted in a dashed line in (f). Plotted data points are taken from the steady-state. and 5t c $ γ k p ). We consider a disk-like shape of a tissue with free boundaries (Fig. 1h). For the single-cell rest shape we choose a regular, compatible hexagon, with rest perimeter p 0 and rest area a 0 À a 0 ¼ 3 p 3ðp 0 =6Þ 2 =2 Á . The cell's perimeter is controlled by an active unit: when reaching criticality (either a critical area a c or critical perimeter p c , showing slight differences between them) a compression force f c is acting for a duration of t c to shorten the perimeter (see "Methods").
The results are qualitatively similar to the 1D case: In the excitable mode (a 0 < a c ) a single pulse is propagating from the perturbed point across the tissue (Supplementary Movie 2). It propagates faster closer to the rim as boundary cells are less constrained and hence more responsive. In the oscillatory mode (a c < a 0 ), after an initial shrinking phase, the system is self-compressed and contraction pulses are propagating in a uniaxial, azimuthal or a spiral fashion (Fig. 4a,  Supplementary Movie 3 left). As in 1D, an expansion front is a precursor to a shrinkage front, while all cells in between are actively contracting. Static contraction and chaotic modes exist as well, resembling similar states in the 1D case (not shown).
The emergent dynamics we observe is analogous in structure to generalized reaction-diffusion dynamics. To observe that, we write a simple continuum model, written in the spirit of the discrete model we presented here (see Supplementary Note 4). The dynamics may be written as: u(r,t) is the displacement vector of each point from its original location r at time t. The diffusion coefficient D depends on viscoelasticity D j j / k=γ and the reaction term, R, depends on the cellular activity f c ; t c ; l c .
Despite the similarities to the 1D setting, a unique feature of the 2D case is the fact that the system is prone to geometrical frustration. An intuitive way to see it, is that rim-cells can only release stresses in the radial direction, but not in the azimuthal one. As a result, rim-cells are not beating like isolated cells, as in 1D. Although rim-cells relax some of their stress faster than bulk cells, complete relaxation pends on "waiting" for the entire system. This can be seen in the radial gradient of strain (Supplementary Movies 3, 4, Supplementary Fig. 2). In addition, bulk cells are relaxing slowly due to viscosity and due to the energy wells, they reach at concave shapes (their perimeter needs to temporarily decrease in order to go back to convexity). The result of these two types of frustration is long intervals of quiescence, with no active contractions. We show that asf increases, quiescence periods increase, while short bursts of activity occur between them (Supplementary Fig. 2).
Rip resistance mechanism. An intriguing feature of the 2D active tissue is its effective mechanical properties under tensile stress. To demonstrate them, we design a numerical experiment where we pull a 2D sheet at two opposite rim points with a constant force  3 left). Shown, is a case of area-trigger (where a cell's critical area a c is lower than its rest-area a 0 ). Similar patterns appear with a perimeter-trigger (where a critical perimeter, p c , is lower than the rest-perimeter p 0 ). Greyscale represents cell area; red dots represent actively contracting cells. The white arrows highlight the propagation of circular and spiral pulses that have features similar to the 1D case: an extension front, followed by a shrinking front, while all cells in between are actively contracting. b Adding a yielding response to the model: when the tension on a cell-cell junction is higher than a critical value, σ s , an immediate softening occurs in both neighboring cells: a reduction in stiffness from k to k s (In the 2D case, the softening can be either in the perimeter stiffness k p or area stiffness k a ). In the sketch, arrows mark contractions, dashed blue lines mark softening. c Snapshots from a simulation of a tissue with both contraction and softening thresholds (Supplementary Movie 3 right). Red dots mark contracting cells, blue hexagons mark softened cells. The white arrows highlight the propagation of pulses with softening fronts accompanying the extension and shrinkage fronts. d-g Testing the response of an excitable tissue (p c > p 0 ) to external stretch in the x-axis: d The pulling configuration: A constant force is pulling uniaxially on the sheet, acting directly on a fixed set of cells on each side (marked with red dashed circles). e The distance between the pulling points L, normalized to the distance at rest L 0 , as a function of time. f We present peak values of mechanical fields per pixel across the tissue. We show the maximal cell perimeter strain (p/p 0 ), area strain (a/a 0 ) and junction stress (σ) throughout the simulation. We compare a passive material, a material with contraction response (EIC) and a material with both EIC and yielding responses. The contractile material dynamically distributes the strains in the tissue but presents high levels of junction stress. A material with the added yield-response increases slightly the levels of cell strain but cuts-off high junction stress values. See also Supplementary Movie 4 and Supplementary Fig. 3. g Histogram-view of all data points in the pulling simulation. (Fig. 4d, Supplementary Movie 4). First, we look at a passive elastic sheet. At equilibrium, high values of stress and strain are focused near the pinching points or along the pulling axis (depending on specifics of the elastic model). Adding viscous behavior either postpones this elastic equilibrium or creates indefinite creep. All these cases impose a threat for tissue cohesion. However, constant pulling on the active tissue results in a dynamic steady state, where high strains vanish, and instead strains oscillate in time around a lower average, exhibiting both tension and compression (Fig. 4e). In addition, maximal strains are evenly distributed in space, as all cells are participating in the dynamics (Fig. 4f). This can be seen also in the histogram of all strain values at all times (Fig. 4g). When the contraction threshold value (a c or p c ) is lower than the actual cell rupture value, the tissue will successfully avoid failure due to high cell strains (see Supplementary Movie 4 and Supplementary Fig. 3).
Adding a second cellular response-softening/yielding due to junction stress-sets an upper bound for junctional stresses as well. We add to the simulation such a cell softening thresholddecrease in the elastic module k p by the arbitrary factor two upon reaching critical junction stress σ s (Fig. 4b, methods). For simplicity, we assume that once the junction stress is restored below σ s , the stiffness returns to its original value. Our results show that the added response does not change the spatiotemporal patterns significantly, except introducing a tailing softening front (Fig. 4c, Supplementary Movie 3). Under external tension, the softening threshold did not change the strain and stress distributions dramatically (Supplementary Movie 4, Supplementary Fig. 3), but it did cut-off high-stress values, trading them for localized high strains, as seen in the pictures of peak values (Fig. 4f) and at the distribution tails (Fig. 4g). When the threshold values p c ; σ s À Á are set below cell and junction rupture values respectively, the tissue avoids failure by suppressing both cell strain and junction stress simultaneously.

Discussion
In this work, we suggest a model for the propagation of contraction pulses in epithelia. The only requirement for pulse propagation is a single cell EIC-contraction due to expansion. Unlike passive elastic retraction after expansion, the model requires the contraction to be active, and include a memory time scale, in order to bring a contracting cell significantly below its rest length despite the overdamped conditions. We show that in order for a contraction to propagate in the excitable-cell mode, the contraction impulse should be strong enough relative to the excitation thresholdtf >ε (a qualitatively similar criterion will exist in the oscillatory-cell mode, with l ss in the role of l 0 ). The resulting nonlinear pulses travel in the excitable media via local energy injections. Despite the overdamped conditions they travel long distances as solitons that do not decay nor change their shape, and evolve dynamically following reaction-diffusion equations for the strain. In biological tissues, and specifically in T. adhaerens, we speculate that such long-range contraction pulses can be used as means for inter-cellular information transfer and tissue communication.
In parallel, we suggest in this work another physiological role for such contraction dynamics in tissues: increasing their resistance to rupture. We refer to the immediate need of a tissue to keep its integrity under acute tension and avoid damage. Under continuous external stretch, passive sheets reach a steady state of strain focusing in specific areas. In contrast, sheets with cellular EIC reach a dynamic state of oscillations between stretch and compression across the entire tissue, diminishing high strain values. This strain homogenization consumes cellular energy and cannot be achieved in a passive material. However, when a cellular sheet is under continuous tension, there is a trade-off between two failure modes. Cellular contractions can prevent rupture by high cell strain, yet it increases cell-cell junction stress. Cell softening prevents cell-cell detachment by reducing junction stress, but it increases cell strain. Using the two responses of contraction and softening, the tissue is spatially distributing, while setting an upper bound, for both cell strain and junction stress, preventing rupture by the two failure modes simultaneously. We show here this simultaneous protection is possible with no parameter fitting.
In biological tissues, it is unclear whether an EIC is triggered by strain, stress, or strain rate. The exact parameter is hard to distinguish experimentally, and currently remains unknown. For our model, we choose cell-strain and junction-stress as triggers for cellular response, as they threat tissue cohesion: cell cortex networks rearrange due to stress, but are probable to fail entirely in high strain. Cell-cell junctions, on the other hand, hardly expand, but may break under stress. In addition, in the presence of drag, a cellular contraction reduces cell strain but increases junction stress. Therefore, we choose high cell strain as the trigger for contraction and high junction stress as the trigger for softening, as two antagonistic actions that directly relate to the two failure modes. We demonstrated here that these feedback loops can indeed promote tissue stability. We choose the active contraction to be governed by a fixed compression force that lasts a fixed duration of time. This is a minimal assumption inspired by our experimental results from T. adhaerens, where we showed the average contraction profile, as well as a narrow distribution of contraction times relative to contraction amplitudes and speeds 26 . We chose the yielding effect to mean an immediate softening by a constant factor and immediate recovery when the junction-stress is restored below the threshold. This is again a simplifying assumption in the lack of accurate data.
Our work opens a family of models, that vary in the specifics of the cellular responses (exact mechanical trigger, level of response, the activation and recovery functions etc). How do these details change the model's results is a subject for future works, both numerical/analytical and experimental. However, as long as the activity thresholds are lower than the rupture values ðε < ϵ rupture ; σ s < σ rupture Þ and as long as the contraction is propagating to neighboring cells ðtf >εÞ the dynamic patterns and the spatial homogenization of the loads would be qualitatively the same.
Our model and the rip resistance mechanism are yet to be tested experimentally in epithelia. Specifically, accurate measurements of response to tension in cells will shed light on the specific parameters of EIC. Measurements at the tissue level will then test our model's predictions regarding the properties of the pulse dynamics. The minimal nonlinear model we presented fits the fast soliton nature of the pulses seen in T. adhaerens but may be applicable to other observed contraction waves, as in drosophila embryo and MDCK monolayers. Specifically, it would be interesting to test our predictions in various embryonic tissues, and in epithelia that is either contractile or prone to high, repeatable mechanical stresses (e.g., heart, lung, gut, bladder, vasculature), We suggest to further investigate dynamics in epithelia of earlydivergent animals, that bring tissue mechanics to the context of the evolution of multicellularity. Alongside the fundamental ability of cells to stay cohesive as tissues, it may shed light on the origin of the excitation-contraction coupling 35 , on mechanical information propagation in living materials, different architectures of information processing, and the evolution of neuro-muscular systems [36][37][38][39][40] . Such early epithelium dynamics may also be seen as embodied calculations yielding "behavior" and supporting physiological needs (e.g., locomotion and navigation, wound healing, and size control). Finally, to the best of our knowledge, our work highlights a new type of active materials, and may inspire the engineering of synthetic analogues that actively resist rupture.
We adopt a dynamical modeling paradigm built around gradient-descent on overdamped equations of motion. Unlike the conventional use of gradient descent, the energy functional in the algorithm is not constant; rather, it changes every time a cell is activated or deactivated. Inspired by T. adhaerens, we take the activation time t c to be shorter than the viscoelastic time K/γ, which is the time scale to approach cellular steady state. Therefore, our gradient descent algorithm will not necessarily reach steady state, and indeed in most interesting cases it does not.
The equation of motion is therefore where l i is the length of cell #i, and the overall free energy E is a sum of the energies of all cells, E t ð Þ ¼ ∑ i ε i t ð Þ. (Note that the cellular energy ε i is different than the cell strain ϵ i mentioned earlier). In 1D, the energy of each cell is given by i takes the value f c when l i > l c and maintains it for a duration of t c , after which it goes back to 0 (Fig. 1d).
To model 2D tissues, we generalize the above framework using a 2D vertexmodel, evolving under the cellular energy function The function F A c i now takes the value f c when a i > a c (or when P > P c ) and maintains it for a duration of t c . In this formulation, the force may be triggered by an area or perimeter threshold but acts on the perimeter.
We define the tensile stress on cell-cell junction for each edge in the 2D hexagonal grid (an edge is attaching two cells, and is composed of two vertices). We evaluate this stress using a three-step procedure. First, the forces on the two relevant vertices arising from the first cell are averaged, and their mean is projected to the axis normal to the edge. The same average and projection is calculated for forces arising from the other cell. Second, for each pair of cells with a shared edge, the two normal forces are subtracted to find the force required to maintain the constraint of confluency. This calculated force is then divided by the edge length to find the force per unit length associated with each junction.
To model tension-induced yielding we add a second cellular response-softening of both cells adjacent to an overstressed junction (we treat cell edges as junctions). The energy functional then becomes: where k P,i = k P if the edges of the ith cell are all under the critical tension σs, and k P;i ¼ k s < k P if one of its junctions is overtensed. In the numerical experiments presented here we choose k s ¼ k P =2. Once the normal stress is reduced below σ s , perimeter stiffness is set back to k p .

Data availability
All relevant data are available from the authors upon request.