Excitable dynamics driven by mechanical feedback in biological tissues

Pulsatory activity patterns, driven by mechanochemical feedback, are prevalent in many biological systems. Here we present a theoretical framework to elucidate the mechanical origin and regulation of pulsatile activity patterns within multicellular tissues. We show that a simple mechanical feedback at the level of individual cells - activation of contractility upon stretch and subsequent inactivation upon turnover of active elements - is sufficient to explain the emergence of quiescent states, long-range wave propagation, and traveling activity pulse at the tissue-level. We find that the transition between a propagating pulse and a wave is driven by the competition between timescales associated with cellular mechanical response and geometrical disorder in the tissue. This sheds light on the fundamental role of cell packing geometry on tissue excitability and spatial propagation of activity patterns.


I. INTRODUCTION
Multicellular systems exhibit a wide range of pulsatile and wave-like patterns during collective migration, development, and morphogenesis [1][2][3].The appearance of these patterns can be attributed to various biochemical factors, depending on the specific phenomenon.These include waves of extracellular signal-related kinase (ERK) [4,5], calcium waves [6], periodic assembly and disassembly of myosin motors [7,8], and the periodic release of chemoattractants [9].Reactiondiffusion models [10][11][12][13] and cellular automaton models [14][15][16] have been widely used to study the mechanisms underlying biochemical pattern formation in multicellular systems.Mechanochemical patterns, on the other hand, have necessitated the development of new classes of models that integrate mechanical forces with chemical reactions [17][18][19][20][21].For instance, the coupling of mechanical and chemical processes is particularly relevant in understanding the spatial propagation of contraction patterns in T. Adhaerens [22], oscillatory morphodynamics in Drosophila amnioserosa tissue [23], collective migration patterns [20] and mechanical waves in expanding MDCK cell monolayers [19,24,25].However, the role of cellular mechanics and geometry in the propagation of mechanochemical signals remains poorly understood.
One commonly observed mechanical feedback motif in cells is stretch-induced contraction, wherein a local stretching deformation triggers the recruitment of active components that induce contraction [4,[26][27][28][29].Recent studies have utilized the concept of stretch-induced contraction to elucidate phenomena such as wave propagation in active elastic media [19,25], contraction pulses in epithelial tissues [30], cell migration patterns in vitro [20,31], cell and tissue morphogenesis [23,32].Specifically, all these studies focused on dynamics in active elastic media, without considering the effects of geometric disorder and viscous dissipation on mechanochemical signal propagation.
In this study, we ask how cellular viscoelasticity and packing geometry regulate the propagation of active stresses at multicellular scales.To this end, we extended the framework of the cellular vertex models [33][34][35][36] to incorporate feedback between cell junction strain and contractility.In addition, viscous dissipation is implemented by continuous strain relaxation in cell junctions.We implement a simple feedback rule in which contractility in cell junctions is activated above a threshold junctional stretch.The junction remains active for a duration commensurate with the turnover rate of active elements.This is followed by a refractory period during which junction contractility remains inactive due to the presence of inhibitors of contractility.As a result of these rules, each cell junction behaves as an excitable unit that can exist in one of three states: active, inactive, and refractory.
Our proposed model elucidates the emergence of longrange propagation of contractile pulses and different patterns of self-sustained traveling waves, such as circular, elliptic, and spiral waves.We show that these tissue-level propagation patterns are controlled by the competition between the timescales associated with active and refractory states of the junction, and the characteristic timescale of junction strain relaxation.To explain these observations analytically, we develop an effective theory of coupled excitable junctions, capable of explaining the emergence of the quiescent, wave-like and pulselike patterns observed in vertex model simulations.Our theoretical framework predicts that shorter junctions promote reactivation of contractility, while larger junctions facilitate the propagation of activity over a broader region of the parameter space.We validate these predictions through simulations of disordered tissues in two dimensions.We find that geometrical disorder promotes sustained wave propagation at the tissue-level, and that the ability of junctions to locally propagate activity increases with its length.

A. Equations of motion
To elucidate the emergent dynamic patterns in an excitable tissue, we use the framework of the vertex model [33][34][35][36], where a monolayer tissue is modeled as a two-dimensional polygonal tiling.The polygons represent the cells, and the arXiv:2310.04950v1[physics.bio-ph]8 Oct 2023 edges represent the cell-cell junctions.Each vertex i, with position r i , is subject to friction with coefficient µ, and elastic forces and inter-cellular tensions arising from a Hamiltonian H.The Hamiltonian governing tissue mechanical energy is given by where the first energy term is a sum over all cells α, and the second term is a sum over the cell-cell junctions defined by the adjacent vertices i and j.The first term in Eq. ( 1) is the elastic energy that penalizes changes in cell area, where K is the bulk elastic modulus, A α and A 0 are the actual and preferred cell areas, respectively.The second term represents an interfacial energy, with tension Λ along each cell junction of length l i j .
Active contractile forces arise at each junction from the actomyosin cortex, generating an active force per unit length, Γ i j (t).Consequently, the active force at each vertex can be written as: As opposed to existing vertex models, here we consider a time-dependent contractility Γ i j (t), whose dynamics depend on junctional strain and memory of mechanical state.To compute the timeevolution of each vertex, we assumed an overdamped limit, such that the equations of motion are given by: The above equation of motion is coupled to the dynamics of junctional strain and contractility, as described below.

B. Cell junction dynamics under stretch-induced contractility
Stretch-induced contractility is a commonly observed regulatory mechanism for controlling the level of active contractile stress in cells [4, 26-28, 37, 38].A local stretch in cell junctions could trigger actin fiber alignment [29,39,40], myosin recruitment [28] and also the activation of the ERK signaling [4] that would promote contractility.We therefore implement a simple model of cellular junctions as viscoelastic materials subject to a strain-tension feedback (Fig. 1A).Here, a local stretch triggers the activation of contractility, which in turn reduces stretch via contractile forces.Additionally, junction strain continuously relaxes over time due to viscous dissipation and contractility undergoes turnover as part of a self-regulatory mechanism.
The mechanical strain in cell junctions is defined as ε i j = l i j − l 0 i j /l 0 i j , where l 0 i j is the rest length of the junction shared by the vertices i and j.The viscoelastic nature of the junctions is modelled through rest length remodeling at a rate k L , leading to continuous strain relaxation [41], Rest length remodeling is a natural consequence of actomyosin networks with turnover, where strained elements are replaced by unstrained ones [26].The feedback between junction strain and contractility is implemented as follows.Each cell junction can exist in one of three states: inactive (Γ i j = 0), active (Γ i j = Γ 0 ), and refractory (Γ i j = 0).While both inactive and refractory states lack contractility, refractory junctions are those that cannot be active for a duration τ ref , representing the timescale associated with the presence of inhibitors of contractility.The rules describing junction state changes are given below (Fig. 1C): • Inactive junctions become active if their strain ε i j exceeds a threshold value ε on .
• Active junctions become refractory after being active for a time period τ act .• Refractory junctions become inactive after a duration τ ref .
As an example, we'll examine the scenario depicted in Fig. 1B-D, where gray, red, and blue junctions represent the inactive, active, and refractory states, respectively.Initially, the junction marked by the black arrow in Fig. 1B is set to an inactive state.Contraction of the neighboring junction induce strain levels exceeding the threshold ε on , thereby triggering ERK signaling activation [4].This activation, in turn, leads to the production of ERK inhibitors [13].The threshold value ε on ensures immunity to small perturbations, allowing junctions to be excitable units.ERK activation induces contractility, changing the junction state from inactive to active, increasing contractility to Γ i j = Γ 0 .The active state persists for a time period τ act (red phase in Fig. 1C), which represents an effective timescale arising from the turnover time of actomyosin, and the inactivation of ERK.The remaining levels of ERK inhibitors keep the junction in a refractory state, in which it can not be re-activated (blue phase in Fig. 1C).Finally, after a time period τ ref , the inhibitors reach low enough levels to take the junction back to the inactive state (final state in Fig. 1C).
Although the activation time period (τ act ) and the refractory time period (τ ref ) are distinct parameters, a recent study [13] has demonstrated that, to adequately describe the observed ERK activity waves, the characteristic timescales for ERK activation and inactivation by inhibitors tend to be similar, typ-ically of the order of a few minutes.For simplicity, we will first focus on the case where τ act = τ ref = τ, and we will refer to this timescale simply as activation period unless otherwise specified.The strain relaxation rate (k L ) and the activation period (τ) will jointly impact the dynamics at both the junction and tissue scales, as elaborated in Sections IIIA-B.Later in Section IIID, we will delve into the effects of varying duration for refractory and active states.

A. Traveling pulse and waves in ordered tissues
To characterize the emergent dynamic states arising from junction-level mechanical feedback, we first simulated an ordered tissue, composed of 260 hexagonal cells, in a box of sides L x ∼ 14 √ A 0 and L y ∼ 18.6 √ A 0 , under periodic boundary conditions.In simulations, we non-dimensionalized force scales by K(A 0 ) 3/2 , length scales √ A 0 , and timescales by µ/(KA 0 ), setting A 0 = 1, K = 1, and µ = 0.636 (min).
We initiate our simulations with a mechanically equilibrated tissue, where all junctions are initially in the inactive state.We then perturb the equilibrium state by manually activating a single junction positioned near the center of the simulation window (Fig. 2A).When the rate of strain relaxation is sufficiently slow, corresponding to a small value of k L , we observe the emergence of two distinct activity patterns depending on the activation period τ.For small τ, we find waves of activity traveling radially outwards, as shown in Fig. 2B (Movie 1).These self-sustaining waves are characterized by alternating rings and regions of red (indicating activity) and blue (indicating refractory) junctions.The tissue activity reaches a steady-state when the wavefront traverses the entire tissue (around t ∼ 5τ in Fig. 2D).Conversely, for larger values of τ, we do not observe self-sustaining waves due to the lack of junction reactivation events.Instead, a single transient activity pulse travels across the tissue (Figs. 2B  and D, Movie 2).Over an extended time period, the tissue eventually becomes entirely inactive.
To quantify the mechanical deformation due to these traveling activity patterns, we calculated the total tissue strain as ε total = ∑ ⟨i, j⟩ ε i j .Fig. 2F shows the dynamics of the total strain for both wave (Fig. 2A) and pulse-like (Fig. 2B) patterns.The pulse causes a positive peak in strain, followed by a negative peak, ultimately returning to zero strain due to mechanical relaxation.Conversely, in the traveling wave pattern, while there is a peak in strain, it eventually stabilizes as a result of activity-induced mechanical fluctuations and the relaxation of strain at the junction level.
These propagating activity states are only observed when the value of k L is sufficiently small.A large k L causes the strain in the neighboring junctions of an active junction to relax before activation can occur, resulting in a quiescent state without any propagation.To quantify the extent of tissue-scale activity, we calculated the maximum fraction of active junctions throughout the simulation.This measurement enables us to identify the phase boundary, determined by the critical value of k L , that separates the regimes with activity propagation (either wave or pulse) from those without propagation (cyan-dashed boundary in Fig. 2G).Moreover, by quantifying the active junction fraction at the final steady state, we can differentiate between the propagating modes, leading to the delineation of the wave-to-pulse phase boundary (white-dashed boundary in Fig. 2G).

B. Effective theory predicts emergent dynamic states
To analytically predict the emergence of excitable pulses, quiescent states, and oscillatory patterns as functions of the strain relaxation rate k L and activation period τ, we developed an effective one-dimensional theory of coupled excitable junctions.Our minimal model consists of three interconnected junctions with fixed boundaries, as shown in Fig. 3A.Each unit comprises an elastic component with a spring constant k and natural length L (representing the one-dimensional version of cell elasticity), connected in parallel with a dashpot of friction coefficient µ, and an active element with contractility Γ 1,2 .If the junction is inactive or refractory then Γ 1,2 = 0, and Γ 1,2 = Γ 0 if the junction is active.These active and elastic elements are connected in parallel with a tensile element with line tension Λ 1,2 .The central junction has a length l 1 (t), while the outer junctions have lengths l 2 (t).The fixed boundary conditions ensure that l 1 (t) + 2l 2 (t) = 3L.Max.Active Junction Fraction The system is initialized in a mechanical equilibrium state, and we perturb it by activating the central junction (Γ 1 = Γ 0 , Γ 2 = 0).We then let the system to evolve following the equations of motion: µdl i /dt = −∂ H eff /∂ l i , Eq. ( 3), and the rules governing the junction states.The effective Hamiltonian governing the system is defined as: We initially considered the scenario of symmetric junctions, wherein Λ 1 = Λ 2 .This corresponds to an ordered tissue where junction tensions and lengths are uniform.To explore the behavior of the system, we numerically solve the system of equations for different values of τ and k L , from t = 0 to t = 2τ.The simulation outcomes can be categorized as follows: i) If the outer junctions remain inactive throughout the simulation, it is classified as a case of No propagation; ii) If the outer junctions become active but the central junction does not reactivate, we observe a single Pulse; and iii) finally, if the outer junctions become active and the central junction re-activates, it falls into the category of Re-activation.Fig. 3B shows the phase diagram of the model in τ-k L phase space showing the emergence of the three outcomes described above.A comparison with the phase diagram for the ordered tissue (Fig. 2E) reveals that the effective model successfully captures both key features of the vertex model: a critical value of k L for propagation of activity, which diminishes for small τ, and a small region of reactivation corresponding to wave-like states.
We then used the one-dimensional effective model (Fig. 3A) to investigate the role of disorder in the propagation of activity.Disorder was introduced by removing the condition of homogeneous line tension, letting Λ 1 ̸ = Λ 2 .First, we analyzed the case Λ 1 < Λ 2 .Due to identical mechanical properties of each junction before activation (other than tension values), the initial equilibrium state featured a central long junction (l 1 > L) flanked by two shorter junctions (l 2 < L) (Fig. 3C).By solving the system of equations numerically, we found that the larger junction (l 1 > L) could propagate activity over a broader region in the (τ, k L ) parameter space, while the re-activation region is substantially diminished.This is because larger junctions produce greater active contractile forces, while shorter neighboring junctions require a lower extension to achieve the strain threshold for activation ε on .Conversely, when Λ 1 > Λ 2 , the opposite behavior was observed.Our effective model thus reveals two main effects of the geometrical heterogeneity (or disorder) on cellular response to active contractility.Large junctions promote propagation of activity, while shorter junctions facilitated reactivation, leading to oscillatory patterns.

C. Tissue disorder promotes self-sustained wave propagation
Motivated by the predictions of the effective model on the impact of geometric heterogeneity, we now investigate the effect of disorder in cell packing geometry on activity propagation in two-dimensional tissue simulations.To this end, we constructed a tissue comprising 208 cells, within a rectangular box with dimensions approximately equal to L x ∼ 14 √ A 0 and L y ∼ 15 √ A 0 , subject to periodic boundary conditions.In these simulations, all mechanical properties at cell and junction levels are the same, with disorder restricted to geometric heterogeneity only.The initial state of the tissue corresponded to a state of mechanical equilibrium, characterized by varying junction lengths and polygon sidedness, as depicted in Fig. 4A.
As previously, we activated a randomly chosen cell junction (see Fig. 4A, t = 0.0τ), and let the tissue evolve from t = 0 to t = 20, for different values of strain relaxation rate k L ∈ (0, 1.5) and activation period τ ∈ (0.2, 5.0).By measuring the maximum active junction fraction (Fig. 4B), we again observe that propagation occurs below a critical k L , for sufficiently large τ.Unlike in ordered tissues (Fig. 2G), wave states are now possible for a wide range of τ values, and propagating solitary pulse only occurs in particular cases with exceedingly large activation periods.Consistent with the predictions of the one-dimensional effective model, we find that the presence of short junctions in disordered tissues promotes junction reactivation, thereby facilitating the emergence of self-sustaining wave-like states.As an illustrative example, Fig. 4A (corresponding to the white dot in Fig. 4B) represents a wave-like state arising in a tissue with parameters (τ = 1.6, k L = 0.7) (Movie 3), which led to pulse propagation in the ordered tissue (Fig. 2E).Interestingly, the junction that is reactivated by the end of an oscillatory cycle need not necessarily be the same one initially chosen for activation.This introduces a non-local effect of disorder in promoting sustained wave-like patterns.We find that the critical k L required for wave propagation increases with the length of the initially activated junction (Fig. 4C), as predicted by the one-dimensional effective model.

D. Controlling the geometry of wavefronts
Our theory and simulations have elucidated that the propagation of activity at the tissue scale is governed by two distinct characteristic timescales of the system: the activation period (τ) (taken to be equal to the refractory time) and the rest length remodeling timescale (k −1 L ).We now investigate the impact of varying the activation (τ act ) and refractory periods (τ ref ) on the resulting dynamic patterns that emerge within the tissue.In particular, we show that the ratio of activation to refractory period controls the geometry of wave patterns.
As in previous simulations, we initiated the simulation by activating a single junction within the tissue, while the remaining junctions remained in the inactive state (see Fig. 5A).We find that the ratio of the activation period to the refractory period, ∆ = τ act /τ ref , controls the wavelength of the propagat- ing waves (see Figs. 5C-D).Specifically, at higher values of ∆, propagating waves fail to materialize, and instead, we observe the presence of a solitary traveling pulse of activity (Fig. 5E).
Inspired by self-sustained spiral patterns observed in excitable systems [16,[42][43][44], we inquired whether we could design an initial state that would break the circular symmetry of the emergent wavefronts.Previous theoretical work using a three-state (inactive, active, refractory) cellular automaton model has shown that spiral waves can emerge from an initial state consisting of a layer of excited cells and an adjacent layer of refractory cells [45][46][47].We therefore initialized our simulations by activating a partial row of junctions, with the neighbors underneath active junctions initialized in a refractory state (Fig. 5B).This initial condition leads to elliptical wavefronts for the case of ∆ = 1 (Fig. 5G, Movie 4).Similar to the case of single junction activation, smaller values for ∆ decrease the wavelength of the traveling wavefront (Fig. 5F).For higher values of ∆ we observe the emergence of a pair of self-sustaining spirals (Fig. 5H, Movie 5).This can be explained as follows.The initial condition of a partial row of active junctions followed by a layer of refractory junctions instigates two distinct pattern formations.Firstly, it leads to the emergence of a propagating wavefront.Secondly, it initiates the formation of two open ends within the wavefront, resulting initially in the addition of an excited element and subsequently in the development of a curved wave segment.This curved segment then propagates outward, adopting a spiral shape.Thus, each open end leads to the emergence of an spiral.Note that due to the periodic boundary conditions in our model, the least number of open ends that can be created is two.These results show that by designing appropriate initial states and disparate timescales for junction activation and refractory periods, the geometry and wavelength of the emergent wavefronts can be precisely controlled in our model.

IV. CONCLUSIONS
In this manuscript, we have introduced a minimal model designed to elucidate the behavior of mechanically-excitable tissues and investigated the role of viscous dissipation and geometrical disorder on tissue-level pattern formation.Traditional approaches for studying dynamic patterns, such as reaction-diffusion equations and cellular automata models, are constrained in their ability to account for spatial deformation in their standard formulations.To address this limitation, we have employed a vertex-based model incorporating agentbased rules at the junction level, representing coarse-grained biochemical reactions that connect junction deformations with the activation of contractility.
Prior work has examined pattern formation in excitable tissues, considering various triggering factors, including celllevel tension [30], cell size [31], and active ERK concentration [4].However, none of these models have considered the effects of viscous dissipation or explored the potential roles of geometrical disorder in pattern formation dynamics.Within our model, which encompasses three distinct junction-level states (active, inactive, and refractory) and considers mechanical strain as the triggering quantity, we observed the emergence of three tissue-level states: quiescent (no propagation), traveling waves, and traveling pulses.These states arise from the interplay between the characteristic timescales associated with junction activation, inactivation, and refractory states.
To explain these emergent dynamics, we have developed an effective junction-scale theory that qualitatively captures the observed behaviors in the vertex model.Our model also provides insights into the impact of geometrical tissue disorder on tissue-level activity states, demonstrating that large junctions promote propagation, while small junctions facilitate re-activation.These predictions have been corroborated through two-dimensional vertex-like simulations, although experimental validation in epithelial tissues remains an avenue for future exploration.
Furthermore, we have demonstrated that the geometry of the emerging traveling wavefronts is influenced by the initial state of junctions and the ratio between the durations of junction active states and refractory states.This intricate interplay results in variations in wavelengths, transitions from waves to pulses, formation of elliptic wavefronts, and pairs of self-sustained spiral wavefronts.The predicted patterns arising from specific initial junction states could potentially be experimentally tested using optogenetic tools to spatially activate myosin contractility [48], ERK [49], and FRET-imaging to visualize the resulting patterns [4].

METHODS
The custom simulation code for the vertex model was implemented using Python 3. Specifically, the simulation code for an ordered tissue featuring traveling wave states can be accessed on GitHub (https://github.com/BanerjeeLab).In implementing T1 transitions, a similar approach to that described in Ref. [50] was adopted.This involves enforcing the creation and instantaneous resolution of a 4-fold vertex whenever a junction's length becomes smaller than l T1 .A newly created junction is set to have l = l 0 = 1.5lT1 .The simulations encompassed tissues of varying sizes, as specified in the respective figure captions.Default model parameters used in the simulations are listed in Table IV.The numerical analysis of the one-dimensional effective model was done in Mathematica 12.The code for this analysis is also available on GitHub (https://github.com/BanerjeeLab).

FIG. 1 .
FIG. 1. Model for stretch-induced contraction and activity dynamics in cell junctions.(A) Junction stretch induces contractility via ERK activation, while contraction reduces stretch.The mechanical strain relaxes at a rate k L , and the active contraction pulse has a lifetime τ. (B) Representative section of a simulated tissue with hexagonal cells, using k L = 0.5 and τ = 2. Red junction color denotes active state, which spreads to neighboring junctions as they are stretched.(C) Junction-level rules to change its state.Gray junction color denotes inactive state, while blue denotes a refractory state.(D) Representative dynamics of junction length, strain and active contractile for an initially inactive junction (gray), indicated by the arrow in panel (B), with ε on = 0.1.

•FIG. 4 .
FIG. 4. Sustained propagation of activity waves in disordered tissues.(A) Snapshots of activity waves traveling across a tissue, using k ) = 0.7) (white dot in (B)).Colored segments represent inactive (gray), active (red), refractory junctions.(B) Phase diagram for the emergence of waves, pulses and quiescent states, varying strain relaxation rate k L and activation period τ. (C) Scatter plot of the critical rest length remodeling rate k L for a fixed τ = 3, as a function of the initially activated junction length, for 100 different chosen junctions.The critical k L is defined as the maximum k L that allow propagating states.Red-vertical line represents the junction length in the ordered hexagonal tissue.Blue-horizontal line represents the critical k L in the ordered tissue.

FIG. 5 .
FIG. 5. Emergent dynamical patterns arising in an ordered tissue composed of 2759 cells, varying τ ref /τ act .Colored segments represent inactive (gray), active (red), and refractory junctions.(A-B) Zoomed in configuration of two different initial conditions for junction states: (A) single junction activation, and (B) partial row activation with neighbours underneath in refractory states.(C-D-E) Collective dynamics arising from the initial condition in (A), for different values of τ ref , with τ act = 0.6.(F-G-H) Collective dynamical states arising from the initial condition (B), for different values of τ ref , with τ act = 1.(C-D-F-G-H) represent self-sustained waves, while (E) represents a solitary traveling pulse.All these simulations consider k L = 0.5.