# Non-linear visco-elasto-plastic rheology of a viscous vertex model

Shalabh Anand, Matthias Merkel

► To cite this version:

Shalabh Anand, Matthias Merkel. Non-linear visco-elasto-plastic rheology of a viscous vertex model. 2026. hal-05530297

HAL Id: hal-05530297

https://hal.science/hal-05530297v1

Preprint submitted on 27 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Non-linear visco-elasto-plastic rheology of a viscous vertex model*

Shalabh K. Anand¹,²,† and Matthias Merkel¹,‡

¹Aix Marseille Univ, Université de Toulon, CNRS, CPT (UMR 7332),
Turing Centre for Living systems, Marseille, France.

²Institut für Theoretische Physik, Universität Innsbruck,
Technikerstrasse 21A, A-6020, Innsbruck, Austria

Morphogenesis involves complex shape changes of biological tissues. Yet, tissue shape changes depend on tissue rheology, which in turn arises from the interplay of large numbers of cells. Here, we link cell- and tissue-scale mechanics by constructing mean-field rheological relations for the vertex model. In contrast to past work in the field, we study a vertex model with an explicit viscous friction. We also include two different cellular mechanisms creating active, anisotropic stresses. Our mean-field model accounts for cell shape and the non-linear elastic and visco-plastic regimes. We validate our results by predicting the response to large-amplitude oscillatory shear. There are several vertex model variants, and comparing to results from the literature, we show that their rheology depends on a number of model details. Our approach should be sufficiently general to construct non-linear mean-field constitutive relations for any cell-based tissue model.

# I. INTRODUCTION

Understanding the mechanical behavior of biological tissues is crucial for deciphering complex processes in developmental biology, tissue engineering, and cancer metastasis. Yet, how tissue-scale mechanics arises from the collective behavior of many cells is subject to active research, including experimental work [1–7], theoretical work [8–11], and both combined [12–15].

There are several tissue models that could be used to examine how tissue dynamics depends on cellular mechanics, i.e. where cells are explicitly represented [16]. These include particle-based models [17, 18], the Cellular Potts Model [19, 20], the vertex model [21–23], and phase-field models [24, 25]. In order to achieve a good compromise between the ability to represent features such as cell shape, not too unrealistic time-dependent dynamics, and reasonable computation time, we use here the vertex model.

Previous work has already studied some aspects of vertex model rheology, including its linear rheology [26–29]. Of note, Tong et al. [29] analytically connect tissue-scale linear rheology to vertex model details using the elegant normal mode formalism, which was also used on particle-based tissue models before [30]. Moreover, much work has focused on tissue-scale elastic properties. For instance, prominently, several vertex model variants are known to exhibit a rigidity transition, which was studied for both ordered [22, 31] and disordered [8] packings. This transition, and the related linear and non-linear elastic properties, have subsequently been studied

in a series of papers [10, 32–41]. However, stress dissipation and the non-linear visco-plastic regime of the vertex model is less well explored.

In the vertex model, there are several ways to implement mechanical stress dissipation on the cell scale. Most often, a substrate friction is implemented, where the frictional force on a vertex is proportional to the vertex velocity. In this case, to probe rheology, care needs to be taken to use the vertex velocity *relative* to a deforming substrate [29]. However, many tissues do not have any substrate – take for instance 3D stem cell aggregates [42], where friction is dominated by internal dissipative processes such as viscous friction. To properly describe such cases in a vertex model, a friction is needed that is Galilean invariant: whether a tissue is moving at some constant velocity or not, nothing should change about the forces the cells experience. Therefore, several models include friction terms that only depend on velocity vector *differences* [29, 43–46]. Yet, without a substrate and in the over-damped limit, one can also globally rotate a tissue at constant rate without any relevant additional force. In other words, appropriate models should also display a “Galilean invariance” with respect to rotation. Both invariances should be required to describe any tissue-internal dissipation processes. But, as far as we know, very little work so far satisfies both invariances for their internal frictions [47–51]. Note that, as an alternative to explicitly taking dissipative processes into account, also quasi-static relaxation obeys these two invariances [26, 52, 53]. In this article, we will explicitly account for the internal dissipation generated in cell-cell interfaces, cell perimeters, and cell areas, broadly following ideas from Ref. [48].

On the tissue scale and at long times, mechanical stresses relax mostly through T1 transitions. This was studied in the vertex model when triggered by externally applied stresses [54–57], active fluctuations [53, 58], or thermal noise [59–62]. Here, to systematically characterize the tissue-scale stress relaxation due to T1 transitions, we will employ a triangle-based method to de-

<figure><img src="image_2.png" /></figure>

<footnote>*</footnote>

<footnote>This work is licensed under a Creative Commons “Attribution 4.0 International” license, https://shalabh.anand@uibk.ac.at</footnote>

<footnote>‡ matthias.merkel@univ-amu.fr</footnote>
compose shear of a cellular network into contributions by cell shape change and T1 transitions [63]. This approach was previously used by one of us to characterize tissue-scale rheology [52], but for a vertex model without explicit internal dissipation. One of us also applied a related approach to emulsion experiments, where we observed quasi-static, plastic stress relaxation. In these emulsions, much of the plastic behavior could be understood in terms of a reversible fraction function [57]. Here, we will follow a similar approach, but for a vertex model with explicit internal dissipation.

Many vertex models also include active terms [27, 64–68]. Notably, we recently studied two different cellular mechanisms of creating active, anisotropic stresses in biological tissues, motivated by active deformations observed in frog embryo explants [69]. Intriguingly, we found that cell shape $Q$ behaved differently in both cases, calling for a mean-field model of non-linear vertex model rheology that explicitly accounts for cell shape.

Here, we implement a vertex model with explicit frictional dynamics that fulfills translational and rotational Galilean invariance. We show how, using the strain decomposition from Ref. [34], tissue shear rheology, i.e. the complex relation between shear stress $\tilde{\sigma}$ and shear rate $\tilde{v}$ can be separated into two relations: (i) shear stress depending on cell shape and shear rate, $\tilde{\sigma}(Q, \tilde{v})$, and (ii) the shear rate contribution by T1 transitions, $R(Q, \tilde{v})$. To quantify these relations, we first use small-amplitude oscillatory pure shear simulations to study the linear rheology. Subsequently, we use constant pure shear rate simulations, to include effects from non-linear elasticity, visco-plastic stress relaxation by T1 transitions, and activity. Our approach can be used to develop a mean-field rheology for any vertex model variant.

## II. METHODS

### A. Vertex model

The 2D vertex model describes an epithelial tissue as a tiling of polygons, each of which represents a cell. The state of the model depends on the positions $\mathbf{r}_i$ of the polygon corners, called vertices, which will be denoted by Latin indices $i, j, \dots$ throughout this article. Elastic forces on the vertices, $\mathbf{F}_i^{\text{el}} = -\partial E/\partial \mathbf{r}_i$, are defined by an energy functional $E(\{\mathbf{r}_i\})$. We use a classical form for the energy functional [22, 31], which reads in dimensionless form:

$$E = \frac{1}{2} \sum_{c=1}^{N} \left[ k_A (a_c - 1)^2 + (p_c - p_0)^2 \right]. \quad (1)$$

Here, the sum is over all $N$ cells $c$. The first term in the sum represents area elasticity of cell $c$, where $k_A$ is an elastic modulus and $a_c$ the cell area. The second term represents cell perimeter elasticity, where $p_c$ denotes the perimeter of cell $c$, and $p_0$ is a preferred perimeter.

We use periodic boundary conditions with periodic box dimensions $L_x \times L_y$. We impose time-dependent pure shear, $\gamma(t)$, such that $L_x(t) = L_{x0}e^{\gamma(t)}$ and $L_y(t) = L_{y0}e^{-\gamma(t)}$ with $L_{x0}L_{y0} = N$ in dimensionless units. Note that this dynamics keeps the total box area constant with $L_x L_y = N$.

For the time-dependent dynamics of the vertex model, we focus on the over-damped limit, which implies force balance on each vertex $i$:

$$0 = \mathbf{F}_i^{\text{el}} + \mathbf{F}_i^{\text{fr}} + \mathbf{F}_i^{\text{act}}. \quad (2)$$

Here, $\mathbf{F}_i^{\text{el}}$ is the elastic force defined right above Eq. (1), $\mathbf{F}_i^{\text{fr}}$ are friction forces, and $\mathbf{F}_i^{\text{act}}$ are active forces. Since many tissues or parts of tissues are not directly in contact with any solid substrate, we want to study a model that fulfills translational and rotational Galilean invariances. This corresponds to the local conservation of linear and angular momentum. We note that the elastic forces $\mathbf{F}_i^{\text{el}}$ already fulfill these invariances.

To introduce friction forces $\mathbf{F}_i^{\text{fr}}$ that fulfill translational and rotational Galilean invariances, we introduce scalar frictional elements. In this article, we introduce frictions on all (i) cell areas, (ii) cell perimeters, and (iii) cell-cell interface lengths. Specifically, we introduce (i) an area stress $\sigma_c^a = (\bar{\eta}^a/a_c)\dot{a}_c$ for any cell $c$, where the dot indicates the time derivative, (ii) a perimeter tension $T_c^p = (\bar{\eta}^p/p_c)\dot{p}_c$ for any cell $c$, and (iii) an interface tension $T_{ij}^\ell = \bar{\eta}^\ell \dot{\ell}_{ij}$ for any interface of length $\ell_{ij}$ between any two connected vertices $i$ and $j$. Following Ref. [48], we obtain for the friction forces:

$$\mathbf{F}_i^{\text{fr}} = - \sum_{m=1}^{N_{\text{fr}}} \frac{\partial g_m}{\partial \mathbf{r}_i} \eta_m \dot{g}_m, \quad (3)$$

where $m$ goes over all $N_{\text{fr}} = 5N$ frictional elements $g_m$. These include: (i) all $N$ cell areas $a_c$, where, for $m = c = 1, \dots, N$, we set $g_m = a_c$, $\eta_m = \bar{\eta}^a/a_c$, and thus $\eta_m \dot{g}_m = \sigma_c^a$; (ii) all $N$ cell perimeters $p_c$, where, for $m = N+c = N+1, \dots, 2N$, we set $g_m = p_c$, $\eta_m = \bar{\eta}^p/p_c$, and thus $\eta_m \dot{g}_m = T_c^p$; (iii) all $3N$ cell-cell interface lengths $\ell_{ij}$, where, for $m = 2N+1, \dots, 5N$, we set $g_m = \ell_{ij}$, $\eta_m = \bar{\eta}^\ell$, and thus $\eta_m \dot{g}_m = T_{ij}^\ell$. Furthermore, we use that $\dot{g}_m$ can be expressed in terms of vertex position movements and changes of the boundary dimensions, i.e. of the pure shear variable $\gamma$: $\dot{g}_m = (\partial g_m/\partial \gamma)\dot{\gamma} + \sum_j (\partial g_m/\partial \mathbf{r}_j) \cdot \dot{\mathbf{r}}_j$. We thus have:

$$F_{i\alpha}^{\text{fr}} = - \sum_{m,j,\beta} \frac{\partial g_m}{\partial r_{i\alpha}} \eta_m \frac{\partial g_m}{\partial r_{j\beta}} \dot{r}_{j\beta} - \sum_m \frac{\partial g_m}{\partial r_{i\alpha}} \eta_m \frac{\partial g_m}{\partial \gamma} \dot{\gamma}, \quad (4)$$

where here and throughout this article, Greek letters denote dimension indices, $\alpha, \beta, \dots \in \{x, y\}$. The first term in Eq. (4) corresponds to what has been obtained previously [48, 49]. The additional second term that we obtain here is due to the shear of the periodic box.

To ensure that also the active forces $\mathbf{F}_i^{\text{act}}$ fulfill translational and rotational Galilean invariance, we use a similar
<figure><img src="image_3.png" /><figcaption>FIG. 1. A schematic showing the implementation of the two active cellular mechanisms: (a) anisotropic interface tension and (b) active crawling forces.</figcaption></figure>

form as for the friction force, Eq. (3):

$$F_i^{\text{act}} = - \sum_{n=1}^{N_{\text{act}}} \frac{\partial h_n}{\partial \mathbf{r}_i} T_n^{\text{act}}, \quad (5)$$

where the sum is over all $N_{\text{act}}$ active elements $h_n$, each of which experiences an active tension $T_n^{\text{act}}$. In this article, we include two different types of activity [69] (Figure 1): (a) We study anisotropic interface tensions, where the active elements are the $N_{\text{act}} = 3N$ cell-cell interface lengths $h_{ij} = \ell_{ij}$ with $n \equiv ij$, between vertices $i$ and $j$, respectively. For each interface $ij$, we introduce an active tension $T_{ij}^{\text{act}}$ that depends on the absolute angle $\phi_{ij}$ of the interface with respect to the $x$ axis:

$$T_{ij}^{\text{act}} = -T_{\text{ait}}^{\text{act}} \cos 2\phi_{ij} \quad (6)$$

Here, $T_{\text{ait}}^{\text{act}}$ is a parameter that scales the anisotropy in the interface tension. (b) We study an active “crawling” mechanism, similar to a known mechanism creating convergent extension in the Xenopus laevis embryo [70]. To this end, within a cell given cell c, we include two pulling forces, which act on the two vertices with minimal and maximal y position, respectively (Figure 1b). In any of those two vertices $i$, the pulling force has magnitude $T_{\text{crawl}}^{\text{act}}$ and acts along the interface pointing away from the cell c (Figure 1b). To ensure local conservation of linear and angular momentum, we include additional forces that act on the cell centers adjacent to vertex $i$ and cell $c$. In terms of Eq. (5), this is equivalent to $N_{\text{act}} = 2N$ active length elements $h_{ci}$ (2 per cell) with tensions $T_{ci}^{\text{act}} = T_{\text{crawl}}^{\text{act}}$.

To numerically solve the time-dependent dynamics, we transform the force balance equation, Eq. (2), into a sparse linear problem:

$$\sum_{j,\beta} \Gamma_{i\alpha j\beta} \dot{r}_{j\beta} = F_{i\alpha}^{\text{eff}}. \quad (7)$$

Here, the friction matrix $\Gamma_{i\alpha j\beta}$ and the effective force $F_{i\alpha}^{\text{eff}}$

depend only on the vertex positions. They are given by:

$$\Gamma_{i\alpha j\beta} = \sum_{m=1}^{N_{\text{fr}}} \eta_m \frac{\partial g_m}{\partial r_{i\alpha}} \frac{\partial g_m}{\partial r_{j\beta}}, \quad (8)$$

$$F_{i\alpha}^{\text{eff}} = -\frac{\partial E}{\partial r_{i\alpha}} - \sum_{m=1}^{N_{\text{fr}}} \eta_m \frac{\partial g_m}{\partial r_{i\alpha}} \frac{\partial g_m}{\partial \gamma} \dot{\gamma} - \sum_{n=1}^{N_{\text{act}}} t_n^a \frac{\partial h_n}{\partial r_{i\alpha}}. \quad (9)$$

To ensure that the friction matrix in Eq. (7) is regular, we also include a very small substrate friction $\eta^s$: $F_i^{\text{fr}} = \dots + \eta^s \dot{r}_i$. This allows to regularize the numerical inversion of the linear problem; yet, we choose $\eta^s$ small enough to leave the physics unaffected.

For any given simulation run, we initialize the system with a Voronoi tiling of $N = 200$ uniformly distributed cell centers, and then minimize the energy $E$. Starting from this initial configuration, we then solve the equations of motion, Eq (7), using an explicit Euler integration scheme with an adaptive time step. Throughout the manuscript, we choose the parameter values $k_A = 10$, $p_0 = 3.0$, $\bar{\eta}^a = \bar{\eta}^p = \bar{\eta}^l = 1$, and $\eta^s = 10^{-3}$. The T1 cutoff is chosen to be $10^{-2}$. The two activity parameters are varied in the ranges $T_{\text{ait}}^{\text{act}} = 0 \dots 2$ and $T_{\text{crawl}}^{\text{act}} = 0 \dots 2$. For each parameter set, we run 8...10 simulations with different random initial configurations.

Generalizing the Batchelor stress formula [71], we compute the mean-field shear stress $\tilde{\sigma}_{xx}$ as:

$$\tilde{\sigma}_{xx} = \frac{1}{2N} \left( \frac{\partial E}{\partial \gamma} + \sum_{m=1}^{N_{\text{fr}}} \eta_m \dot{g}_m \frac{\partial g_m}{\partial \gamma} + \sum_{n=1}^{N_{\text{act}}} t_n^a \frac{\partial h_n}{\partial \gamma} \right). \quad (10)$$

In the following, we will focus on the $xx$ component exclusively. We will use the simplified notation $\tilde{\sigma} \equiv \tilde{\sigma}_{xx}$.

## B. Strain decomposition

We use the framework from Ref. [63] to decompose the overall tissue shear rate tensor $\tilde{v}_{\alpha\beta}$, i.e. the anisotropic part of the strain rate tensor, into contributions by cell shape changes and cell rearrangements:

$$\tilde{v}_{\alpha\beta} = \frac{DQ_{\alpha\beta}}{Dt} + R_{\alpha\beta}. \quad (11)$$

Here, $Q_{\alpha\beta}$ is a symmetric, traceless tensor that describes the average cell shape anisotropy, $DQ_{\alpha\beta}/Dt$ is a co-rotational derivative of that tensor, and $R_{\alpha\beta}$ describes the shear rate contribution by T1 transitions and correlation effects. The latter arise from interactions between local rotation and area expansion with fluctuations in $Q_{\alpha\beta}$ [63].

Importantly, all quantities in Eq. (11) can be quantified in the simulations. Briefly, the overall shear rate $\tilde{v}_{\alpha\beta}$ corresponds to the changes imposed through the boundaries of the periodic box [52]. In this article, we have $\tilde{v}_{xx} = \dot{\gamma}$ and $\tilde{v}_{xy} = 0$. Following Ref. [63], to quantify average cell shape $Q_{\alpha\beta}$, we first construct a trian- gulation of the cellular network that corresponds to its
dual, where the triangle corners are the cell centers. The shape of each triangle is then characterized by a symmetric, traceless tensor, from which $Q_{\alpha\beta}$ is computed as their area-weighted average. Finally, while $R_{\alpha\beta}$ can be computed explicitly from the triangulation [63], here, we compute it using Eq. (11).

In this paper, we focus on pure shear deformation, and thus the $xx$ components of all symmetric, traceless tensors. To simplify notation, we set: $\tilde{v} \equiv \tilde{v}_{xx} \equiv \dot{\gamma}$, $Q \equiv Q_{xx}$, and $R \equiv R_{xx}$.

### III. RESULTS

In order to obtain a mean-field rheological model for our vertex model, we will focus on the symmetric-traceless part of the stress tensor, relating the anisotropic strain rate, i.e. the pure shear rate $\tilde{v}$, to the shear stress $\tilde{\sigma}$. Yet, the instantaneous shear stress $\tilde{\sigma}$ in the vertex model does not only depend on the shear rate, but also on cell shape [52, 67]. Thus, we aim to construct a shear stress function $\tilde{\sigma}(Q, \tilde{v})$.

For a complete description of the relation between $\tilde{v}$ and $\tilde{\sigma}$, the relation $\tilde{\sigma}(Q, \tilde{v})$ needs to be complemented by an additional relation that allows to predict $Q$. To this end, we use the shear rate decomposition, Eq. (11), in the form:

$$ \dot{Q} = \tilde{v} - R(Q, \tilde{v}). \qquad (12) $$

Note that we did not include the corotational term from Eq. (11) here, since in our pure shear simulations, there is no average vorticity. Eq. (12) can be used to predict the time-dependent $Q$ given the imposed shear rate $\tilde{v} \equiv \dot{\gamma}$, once the rearrangement contribution $R$ is known. We postulate that $R$ mostly depends on the average $Q$ and the strain rate $\tilde{v}$. A similar assumption allowed to reconstruct droplet shapes in emulsions relatively well [57].

Because cell shape changes are generally time reversible, while cell rearrangements are not, Eq. (12) decomposes tissue shear into what we call a reversible fraction $\dot{Q}/\tilde{v}$ and an irreversible fraction $R/\tilde{v}$, such that $\dot{Q}/\tilde{v} + R/\tilde{v} = 1$ [57].

In the following, we will first construct the relation $\tilde{\sigma}(Q, \tilde{v})$ in section III A. Afterwards, we fit phenomenological relations to $R(Q, \tilde{v})$, which describes yielding in the vertex model, in section III B. This approach allows to study the effect of vertex model details on the shear stress and the yielding behavior separately. Also, it allows us to include cell shape in the rheological equations, which can often be quantified both in simulations and experiments. Finally, in section III C, we apply our mean-field model to predict cell shape and shear stress in response to large-amplitude oscillatory shear.

<figure><img src="image_4.png" /><figcaption>FIG. 2. Linear shear rheology. (a) Shear protocol, $\gamma(t)$, and corresponding cell shape $Q(t)$ and shear stress $\sigma(t)$. Parameter values: $\gamma_0 = 10^{-4}$ and $\omega = 10^{-1}$. (b) Storage and loss moduli, $G'$ and $G''$, as a function of $\omega$ for strain amplitude $\gamma_0 = 10^{-7}$. Error bars, which are mostly smaller than symbols size, show the standard error of the mean.</figcaption></figure>

#### A. Constitutive relation for shear stress, $\tilde{\sigma}(Q, \tilde{v})$

##### 1. Linear visco-elastic regime

We first study the linear regime of $\tilde{\sigma}(Q, \tilde{v})$, i.e. the regime for small $Q$ and $\tilde{v}$, by imposing oscillatory shear, $\gamma(t) = \gamma_0 \sin(\omega t)$ (Figure 2a, blue curve). We choose a sufficiently small amplitude so that there are no rearrangements, $R = 0$. Thus, according to Eq. (12), $Q(t) = Q_0 + \gamma(t) = Q_0 + \gamma_0 \sin(\omega t)$, where $Q_0$ is a small average cell shape anisotropy that could be present in the initial, energy-minimized, configuration (shift in $Q$ curve in Figure 2a). The shear rate is given by $\tilde{v} \equiv \dot{\gamma} = \gamma_0 \omega \cos(\omega t)$. Finally, for small amplitude strain, the shear stress is also sinusoidal with some unknown phase shift, i.e. it is the superposition of a sine and a cosine. Hence, for sufficiently-small-amplitude oscillatory shear, the functional form of the time-dependent part of $\tilde{\sigma}(Q, \tilde{v})$ is known à priori; it has to be of the form:

$$ \tilde{\sigma}(Q, \tilde{v}) = 2G_0Q + 2\eta_0\tilde{v}, \qquad (13) $$

where the unknown prefactors $G_0$ and $\eta_0$ are linear shear modulus and shear viscosity, respectively. The results for $\gamma_0 = 10^{-7}$ are shown in Figure 2b: The storage modulus corresponds to the shear modulus, $G' = G_0$, and the loss modulus is proportional to the shear viscosity, $G'' = \eta_0 \omega$. We find that both shear modulus and viscosity are constant for small $\omega \lesssim 1$. For larger $\omega$, the shear modulus starts to increase, while the viscosity decreases. This is consistent with earlier work [29], where it has been shown to be due to the internal degrees of freedom not having
<figure><img src="image_5.png" /><figcaption>FIG. 3. Non-linear elastic regime. (a) Shear protocol $\gamma(t)$ and measured $Q$ and $\sigma$. (b) Average steady-state stress $\tilde{\sigma}$ as a function of the average steady-state cell shape $Q$. The dashed line shows a fit to Eq. (14). Fixing the value $G_0 = 0.82$ from the linear rheology (Figure 2b), we obtain $G_2 = 0.86$. Error bars show the standard error of the mean.</figcaption></figure>

sufficient time to fully relax during an oscillation period. Here, we focus on the slow regime, where shear modulus and viscosity are constant, with values $G_0 = 0.82$ and $\eta_0 = 1.16$.

## 2. Non-linear elastic regime

We next extend $\tilde{\sigma}(Q, \tilde{v})$ to the non-linear elastic regime. To this end, we carry out simulations with an initial phase with a constant shear rate, $\dot{\gamma} > 0$, followed by a relaxation phase with zero shear rate, $\dot{\gamma} = 0$ (Figure 3a). To allow for sufficient pure shear that can be applied to the system, we start from a vertically elongated state with $L_{0x} = 10\sqrt{2}/e$ and $L_{0y} = 10e\sqrt{2}$, where $e$ is the Euler constant. We probe the non-linear, elastic behavior of the shear stress, i.e. $\tilde{\sigma}(Q, \tilde{v} = 0)$, using the final states in these simulations. Varying the duration of the initial phase and its shear rate allows us to vary $Q$ up to $Q \approx 0.4$.

The results are plotted in Figure 3b. We find that the data can be fitted to

$$ \tilde{\sigma}(Q, \tilde{v} = 0) \equiv \tilde{\sigma}_{\text{el}}(Q) = 2G_0Q + 4G_2Q^3, \quad (14) $$

where, imposing the known value for $G_0$, we obtain $G_2 = 0.86$. In our phenomenological fit function, we only included odd terms, since $\tilde{\sigma}$ should flip sign whenever $Q$ flips sign.

## 3. Non-linear elasticity with linear viscosity

Next, we combine the non-linear elastic regime with a finite shear rate by carrying out constant-shear-rate simulations (Figure 4a). We again start from the vertically elongated state.

The resulting final values for shear stress and average cell shape are plotted in Figure 4b. Probing the dependence of $\tilde{\sigma}$ on both $Q$ and $\tilde{v}$ is not obvious, since $Q$ and $\tilde{v}$ are correlated in these simulations (Figure 4b inset). We thus probe in how far a simplified relation for the shear stress holds, i.e. in how far it can be described by non-linear elastic behavior combined with a linear viscosity acting in parallel:

$$ \tilde{\sigma}(Q, \tilde{v}) \equiv \tilde{\sigma}_{\text{passive}}(Q, \tilde{v}) = \tilde{\sigma}_{\text{el}}(Q) + 2\eta\tilde{v}. \quad (15) $$

We first tested whether this equation holds using $\tilde{\sigma}_{\text{el}}$ from the previous section (section III A 2) and $\eta = \eta_0$ from the linear shear rheology (section III A 1). In Figure 4c, we plot $\tilde{\sigma}(Q, \tilde{v}) - \tilde{\sigma}_{\text{el}}(Q)$ over the shear rate $\tilde{v}$, and find indeed a linear relation between both (black filled circles). However, the predicted relation using $\eta = \eta_0$ has a much smaller slope (orange dashed line). In other words, for simulations where the tissue undergoes steady pure shear, the measured shear stress is larger than expected.

There are – at least – two possibilities for this deviation from the expected result, either the viscous or the elastic stress was under-estimated. First, viscous stress may have been under-estimated. Indeed, the blue line in Figure 4c represents a linear fit with much larger viscosity of $\eta_s \approx 3.12$. This would suggest that the viscosity $\eta$ in Eq. (15), which acts in parallel to the elasticity, depends on whether there are cell rearrangements in the tissue.

Second, we may have under-estimated the elastic stress, $\tilde{\sigma}_{\text{el}}(Q)$. Indeed, in the previous section, we estimated $\tilde{\sigma}_{\text{el}}(Q)$ based on values of $Q \lesssim 0.4$ (Figure 3b), which we used to predict the shear stress for $Q \gtrsim 0.5$ in this section (Figure 4b). Thus, maybe higher-order terms in $\tilde{\sigma}_{\text{el}}(Q)$ become relevant for $Q \gtrsim 0.5$. To obtain the corresponding $\tilde{\sigma}_{\text{el}}(Q)$ in Figure 4d, we plot $\tilde{\sigma}(Q, \tilde{v}) - 2\eta_0\tilde{v}$ over $Q$, including both the data from the $\tilde{v} = 0$ protocol of the previous section (blue squares) and the $\tilde{v} > 0$ protocol of this section (black circles). Taken together, this data can be well described by a fifth-order polynomial for $\tilde{\sigma}_{\text{el}}(Q)$ (fit shown as black solid line).

Thus, both hypotheses describe the stationary-state data reasonably well. To distinguish between them, we predict the time evolution of the shear stress $\tilde{\sigma}$ given the known time evolutions of $Q$ and $\tilde{v}$, and compare to the $\tilde{\sigma}$ measured in the simulations. To predict $\tilde{\sigma}$ for the first hypothesis, i.e. the shear viscosity is larger when there are T1 transitions, we use Eq. (15) with $\tilde{\sigma}_{\text{el}}(Q)$ from the 3rd-order polynomial fit in Figure 3b in the previous section. Furthermore, to make the viscosity dependent on whether there are T1 transitions, we use the irreversible fraction $R/\tilde{v}$ to linearly interpolate between
<figure><img src="image_6.png" /><figcaption>FIG. 4. Non-linear elasticity with linear viscosity. (a) Shear protocol $\gamma(t)$ and measured $Q$ and $\sigma$. (b) Average steady-state stress $\sigma$ and $Q$ as a function of shear rate $\tilde{v}$. (c) The part of the shear stress $\tilde{\sigma}$ not captured by the elastic stress $\tilde{\sigma}_{el}$ from Figure 3b increases linearly with the shear rate (black circles). Yet, the associated viscosity $\eta_s$ (blue line) is substantially larger than the viscosity $\eta_0$ obtained from the linear rheology (orange dashed line, cf. Figure 2b). (d) Conversely, the part of the shear stress $\tilde{\sigma}$ not captured by the viscous stress expected from the linear rheology, i.e. $\tilde{\sigma} - 2\eta_0\tilde{v}$, is larger than the elastic stress function extrapolated from Figure 3b (blue dashed line). The blue data points are the same as in Figure 3b and the black data points are those from the constant shear protocol, panel b. We take the union of both data sets and fit the 5th-order polynomial $\tilde{\sigma}_{el} = 2G_0Q + 4G_2Q^3 + 6G_4Q^5$ (black solid line). We obtain $G_2 = -0.85$ and $G_4 = 5.4$ when imposing $G_0 = 0.82$ from the linear rheology. (e) Comparison of the time-dependent shear stress $\tilde{\sigma}$ measured in the simulations (blue and orange solid lines) with predictions for two different shear rates $\tilde{v}$ (black solid and dashed lines). Solid lines show predictions assuming the scenario in panel c: the elastic stress from Eq. (14) and Figure 3b with an increased viscosity when there are T1 transitions, Eq. (16). To this end, $R(t)$ is computed using a binned average. Dashed lines show predictions assuming the scenario from panel d: the elastic stress from the 5th-order polynomial with the linear-rheology viscosity $\eta_0$. Shaded regions indicate the standard error of the mean.</figcaption></figure>

small-amplitude and steady-shear viscosities:

$$ \eta = \eta_0 + \frac{R}{\tilde{v}}(\eta_s - \eta_0). \qquad (16) $$

To predict $\tilde{\sigma}$ for the second hypothesis, i.e. the elastic stress becomes much larger beyond $Q \approx 0.4$, we use Eq. (15) with $\tilde{\sigma}_{el}(Q)$ from the 5th-order polynomial fit in Figure 4d, and we use the viscosity from the linear rheology, $\eta = \eta_0$. The resulting predictions are shown in Figure 4e by solid (first hypothesis) and dashed (second hypothesis) lines, respectively. Comparing both predictions to the measured stress (shaded area in Figure 4e), we find that the first hypothesis fits the data better than the second one. This suggests that the viscosity $\eta$ is strongly increased when the tissue is undergoing T1 transitions.

Note that we find an initial transient deviation between

simulation and predictions in Figure 4e: Both predictions show a finite stress already initially – this is the expected viscous stress, due to the tissue being sheared at a constant rate. Yet, in the simulations, the stress is building up over time. We believe that this stress buildup is related to the small, but finite relaxation time of the internal degrees of freedom [29]. This already created the deviation for large frequencies in the linear shear rheology (Figure 2b), and it is not captured by the simple rheological relations that we aim for in this article.

## 4. Activity

We next study how the shear stress $\tilde{\sigma}(Q, \tilde{v})$ is changed when including two different active, anisotropic forces
<figure><img src="image_7.png" /><figcaption>FIG. 5. Shear stress created by activity. The part of the shear stress $\bar{\sigma}$ not captured by the passive stress from Eq. (15) as a function of the respective activity parameter for varying shear rate $\tilde{v}$. The dashed lines represent linear fits. (a) Active, anisotropic interface tensions and (b) active crawl forces. (insets) The respective total shear stress $\tilde{\sigma}$ as a function of $\tilde{v}$ for different activity strengths. The dashed lines are the consistency checks using Eqs. (17) and (18), respectively.</figcaption></figure>

in the vertex model: anisotropic interface tensions, and crawling forces. To this, end, we carried out constant-shear-rate simulations as in section III A 3, while including either a finite interface tension anisotropy $T_{ait}^{act}$ or a finite crawl force $T_{crawl}^{act}$.

We found that in both cases, the activity changed $\tilde{\sigma}(Q, \tilde{v})$ in a way that shifts the stress curve (see insets to Figure 5a,b). To see how this shift depends on the activity parameter, we plotted $\tilde{\sigma}(Q, \tilde{v}) - \tilde{\sigma}_{passive}(Q, \tilde{v})$ over the activity parameter, respectively, and we found a linear relation in both cases (Figure 5a and b). This suggests the following rheological relations, for active, anisotropic interface tensions:

$$ \tilde{\sigma}(Q, \tilde{v}) = \tilde{\sigma}_{passive}(Q, \tilde{v}) - \rho_{ait} T_{ait}^{act}, \quad (17) $$

and for an active crawling force

$$ \tilde{\sigma}(Q, \tilde{v}) = \tilde{\sigma}_{passive}(Q, \tilde{v}) - \rho_{crawl} T_{crawl}^{act}, \quad (18) $$

where $\tilde{\sigma}_{passive}(Q, \tilde{v})$ is given by Eq. (15) with elastic stress and viscosity given by Eqs. (14) and (16), respectively.

Linear fits yield the values $\rho_{ait} = 0.96$ and $\rho_{crawl} = 0.8$ for the activity coefficients (Figure 5a and b).

## B. Constitutive relation for rearrangements, $R(Q, \tilde{v})$

### 1. Without activity

To probe the dependence of $R$ on $Q$ and $\tilde{v}$, we use the constant-shear-rate simulations from section III A 3. For each simulation run and time point, we use the $Q$ and $R$ data, and bin them with respect to $Q$, while averaging $R$ within each bin. The resulting data is shown in Figure 6a and inset.

In past work, where sufficiently strong fluctuations were included, $R$ depended mostly only on $Q$ [52, 53]. However, this is different in our case, where the function $R(Q)$ strongly depends on the shear rate (Figure 6a inset). However, the irreversible fraction $R(Q, \tilde{v})/\tilde{v}$ depends much more weakly on $\tilde{v}$ (Figure 6a). This is expected for a regime where the shear rate is smaller than the inherent relaxation times. In other words, we are closer to the quasi-static limit, where $R(Q, \tilde{v})/\tilde{v}$ should become independent of the shear rate [57].

We find that for a given shear rate $\tilde{v}$, we can phenomenologically describe the dependence of $R$ on $Q$ by exponentials (Figure 6a, dashed lines):

$$ \frac{R(Q, \tilde{v})}{\tilde{v}} = \frac{1}{2} \exp \left( b(\tilde{v}) [Q - Q_0(\tilde{v})] \right), \quad (19) $$

where the two fit parameters $b$ and $Q_0$ depend on $\tilde{v}$ (Figure 7a and b, blue data points). Note that, $b$ and $Q_0$ are odd in $\tilde{v}$, i.e. $b(-\tilde{v}) = -b(\tilde{v})$ and $Q_0(-\tilde{v}) = -Q_0(\tilde{v})$. This is because the passive system is invariant with respect to a rotation by $\pi/2$, and thus $R/\tilde{v}$ has to remain the same when simultaneously flipping the signs of $\tilde{v}$ and $Q$.

### 2. With activity

When including activity, we find again that $R/\tilde{v}$ can be fitted by exponentials (Figure 6b, c), at least as long as the activity is not too large. Here, we thus focus on small activity parameters, $T_{ait}^{act} < 1$ and $T_{crawl}^{act} < 1$. We find that the values of $b$ are largely independent of activity (Figure 7a and c), while $Q_0$ does depend on activity (Figure 7b and d). To further simplify the dependence of $Q_0$ on $\tilde{v}$ and activity, we Taylor-expand it in terms of the activity parameters. For active, anisotropic interface tensions:

$$ Q_0(\tilde{v}, T_{ait}^{act}) = Q_0^{passive}(\tilde{v}) + T_{ait}^{act} \Delta q_0^{ait}(\tilde{v}). \quad (20) $$

Note that while $b(\tilde{v})$ and $Q_0^{passive}(\tilde{v})$ are odd functions, $\Delta q_0^{ait}(\tilde{v})$ has to be an even function. This is because a rotation by $\pi/2$, which should leave $R/\tilde{v}$ invariant, corresponds to a simultaneous sign flip of $\tilde{v}$, $Q$, and $T_{ait}^{act}$ (cf.
<figure><img src="image_8.png" /><figcaption>FIG. 6. Irreversible fraction $R(Q, \tilde{v})/\tilde{v}$ as a function of cell shape $Q$ and shear rate $\tilde{v}$, for (a) passive vertex model tissue, (b) with active, anisotropic interface tension, where $T_{ait}^{act} = 0.2$, and (c) with active crawl force, where $T_{crawl}^{act} = 0.2$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIG. 7. Parameters of the irreversible fraction function fits in Figure 6 and Eq. (19). (a,c) $b$ and (b,d) $Q_0$, and their dependence on shear rate $\tilde{v}$ and activity. (a,b) For active, anisotropic interface tension. (c,d) For active crawling forces. Dashed lines in (b, d) show $Q_0$ computed from Eqs. (20) and (21) using the $\Delta q_0(\tilde{v})$ functions for activity strength 0.5, respectively.</figcaption></figure>

Eq. (6)). Similarly, the relation for the crawl forces read:

$$Q_0(\tilde{v}, T_{\text{crawl}}^{\text{act}}) = Q_0^{\text{passive}}(\tilde{v}) + T_{\text{crawl}}^{\text{act}} \Delta q_0^{\text{crawl}}(\tilde{v}). \quad (21)$$

While the symmetry argument cannot be easily made for the active crawling force, in the following we will test the assumption that $\Delta q_0^{\text{crawl}}(\tilde{v})$ is even as well. We find that Eqs. (20) and (21) characterize $Q_0$ quite well, as demonstrated by the dashed lines in Figure 7b and d. These lines were computed using Eqs. (20) and (21) with $\Delta q_0(\tilde{v})$ taken from the $Q_0$ data at activity strengths 0.5, respectively.

<figure><img src="image_10.png" /><figcaption>FIG. 8. Measured steady-state cell shape $Q_s$ and its prediction as a function of shear rate $\tilde{v}$ for different values of (a) active, anisotropic line tension and (b) active crawling force. Closed symbols show the simulation data and dashed lines are the predictions based on Eq. (22).</figcaption></figure>

### C. Predictions for cell shape $Q$ and shear stress $\tilde{\sigma}$

We now use the relations we obtained for $\tilde{\sigma}$ and $R$ to make predictions for cell shape $Q$ and shear stress $\tilde{\sigma}$ for some given shear protocol. Specifically, for $\tilde{\sigma}$ we use Eqs. (14)–(18), while for $R$, we use Eqs. (19)–(21) together with the $b(\tilde{v})$, $Q_0^{\text{passive}}(\tilde{v})$, and $\Delta q_0(\tilde{v})$ data from plots Fig. 7, taking into account their respective symmetry properties.
<figure><img src="image_11.png" /><figcaption>FIG. 9. Predictions for (a,b,c) $Q$ and (d,e,f) $\tilde{\sigma}$ for large-amplitude oscillatory shear. Blue solid lines show simulation data, while dashed lines show predictions. Shaded regions represent the standard error of the mean. (a,d) Passive tissue, (b,e) active, anisotropic interface tension at $T_{ait}^{act} = 0.2$, and (c,f) active crawling at $T_{crawl}^{act} = 0.2$.</figcaption></figure>

1. Steady-state cell shape

First, as a consistency check for the obtained $R$ functions, we predict the steady-state value for cell shape, $Q_s$, for the constant-shear-rate simulations. At long times, under steady shear, all of the shear will be accommodated by cell rearrangements, i.e. the irreversible fraction is $R/\tilde{v} = 1$. Insertion into Eq. (19) and inverting yields:

$$Q_s = Q_0 + \frac{\log 2}{b}. \qquad (22)$$

Indeed, comparison to the values obtained in the simulation show an excellent fit for passive tissues, and tissues with active, anisotropic interface tension and active crawling forces (Figure 8a and b). The increase of $Q_0$ with shear rate (Figure 7b and d) also explains why we observed a larger steady-state cell shape, $Q_s$, in the constant-shear-rate simulations (section III A 3), where $Q_s \gtrsim 0.5$, than in the two-phase shear simulations with final shear rate zero (section III A 2), where $Q_s \lesssim 0.4$. Moreover, this also explains why the two active cellular mechanisms can create different steady-state cell shapes: We find that for the crawling mechanism, $\Delta q_0$ has a higher negative amplitude than for the anisotropic interface tensions, in particular for small $\tilde{v}$. This translates into smaller $Q_0$, and thus a smaller steady-state value $Q_s$. As a consequence, for larger activity magnitude as used in our earlier work [69], for the crawling mechanism cell shapes can even become perpendicular to the direction of shear.

2. Large-amplitude oscillatory shear

After this consistency check, we predict the result of a completely different type of simulation – large-amplitude oscillatory shear. For frequency $\omega = 0.1$ and the amplitude $\gamma_0 = 1$, we plot the trajectories for $Q(t)$ and $\tilde{\sigma}(t)$ in Fig. 9 for the purely passive case (panels a and d) as well as the two active cases (panels b, c, e, and f). We find that our predictions (black dashed lines) matched the respective simulated trajectories (blue solid lines) relatively well in all cases. There are some deviations before the reversal of the shear direction, when the tissue is yielding. This may be due to the fact that we apply our theory gauged from constant-shear-rate simulations to a simulation with a varying shear rate. The observed deviations may thus arise from some form of memory of the past deformation history that is not fully captured by the average cell shape $Q$ alone.

IV. DISCUSSION

We have developed a mean-field rheological model for the vertex model. While this paper is not the first one studying vertex model rheology [8, 10, 22, 26–29, 31–40, 52, 53, 72–74], we note that there are important differences between our work to what was done before.

First, we follow a systematic approach based on a decomposition of the tissue shear rate $\tilde{v}$ into a reversible contribution by change of cell shape $Q$ and an irreversible
contribution by cell rearrangements *R* [52, 57, 63]. This
allows us to split the question of obtaining a mean-field
rheology into two separate relations, one about the re-
arrangements, *R*(*Q*, *ṽ*), and one about the shear stress,
*ṽ*(*Q*, *ṽ*). Also, this automatically allows to include cell
shape *Q* in the rheological relations, which can often be
determined more easily in experiments than shear stress
[5, 12, 75, 76].

Second, we studied here a vertex model whose dis-
sipation dynamics fulfills translational and rotational
Galilean invariance, i.e. with exclusively tissue-internal
friction instead of the very common friction with an ex-
ternal substrate. This is important, since there are many
tissues without external substrate. Moreover, it is known
that the behavior of active materials can strongly de-
pend on the type of friction, e.g. an external friction can
suppress a well-known active-matter instability [77–81].
Only few vertex model implementations include a tissue-
internal friction without bulk forces or bulk torques, and,
as far as we know, this is the first time the rheology of
a vertex model with internal friction is systematically
probed.

Finally, while there are several papers studying the lin-
ear vertex model rheology, or that focus on the elastic be-
havior only, here we discussed the full active, non-linear
visco-elasto-plastic rheology of the vertex model.

In our mean-field rheological model, the apparent ver-
tex model viscosity on long times results from two ef-
fects: First, a larger shear rate leads to larger steady-
state cell shape $Q_s$, and thus a larger elastic stress $\tilde{\sigma}_{el}$ (cf.
Eq. (14)). This is the only mechanism that created an ef-
fective viscosity in some past work [26, 52, 53]. Yet, in our
paper, where we explicitly account for internal friction,
we also account for an explicit viscosity $\eta$ (cf. Eq. (15)).
A key result of our work is that this viscosity strongly
depends on whether the system shears visco-elastically
due to cell shape changes, or visco-plastically due to cell
rearrangements. In the latter case, the viscosity is three
times as large as in the former. We suspect that the ad-
ditional shear stress is created by the relaxation of cell

shapes right after the occurrence of T1 transitions.

Our work further illustrates that the vertex model
rheology strongly depends on the type of vertex model
dynamics used. For instance, studies with a dynamics
purely based by line tension fluctuations showed a linear
relation between shear stress and cell shape at least up
until $Q \approx 0.4$ [52, 53], while we find a non-linear elas-
tic stress already for these values. Moreover, in these
studies the fluctuations created a very different yielding
behavior, with linear $R(Q)$ for large fluctuation ampli-
tudes, while we observe exponential relations $R(Q)$ that
also depend on the shear rate $\tilde{v}$.

There are several ways in which our work can be ex-
tended. For instance, we needed to describe the irre-
versible fraction $R/\tilde{v}$ by some phenomenological fit func-
tion. In the future, it will be interesting to study how
these curves emerge, for instance using some mesoscopic
model. Such a mesoscopic model may also help refine
predictions for non-constant-shear-rate simulations (Fig-
ure 9). Moreover, we occasionally observed partial crys-
tallization, and future work could suppress this using
polydisperse cell parameters in the vertex model energy.
Finally, while the rheological model developed here and
other work focuses on the mean-field limit, it would be
interesting to include gradient terms in such a description
to describe the interactions of neighboring tissue regions.

ACKNOWLEDGMENTS

The project leading to this publication was supported
by the grant RobustTissue attributed to M.M. by the
French National Research Agency (ANR-22-CE30-0039).

Following the *Rights Retention Strategy of Plan S*, a
CC-BY 4.0 public copyright license has been applied by
the authors to the present document and will be applied
to all subsequent version, including the Author Accepted
Manuscript arising from this submission. This does not
apply to the Version of Record, for which this paragraph
can be removed.

[1] A. R. Harris, L. Peter, J. Bellis, B. Baum, A. J. Kabla, and G. T. Charras, Proceedings of the National Academy of Sciences **109**, 16449 (2012).

[2] F. Serwane, A. Mongera, P. Rowghanian, D. A. Kealhofer, A. A. Lucio, Z. M. Hockenbery, and O. Campàs, Nature Methods **14**, 181 (2017), ISSN 15487105.

[3] A. Mongera, P. Rowghanian, H. J. Gustafson, E. Shelton, D. A. Kealhofer, E. K. Carn, F. Serwane, A. A. Lucio, J. Giammona, and O. Campàs, Nature **561**, 401 (2018), ISSN 0028-0836.

[4] N. Khalilgharibi, J. Fouchard, N. Asadipour, R. Barrientos, M. Duda, A. Bonfanti, A. Yonis, A. Harris, P. Mosaffa, Y. Fujita, et al., Nature Physics **15**, 839 (2019), ISSN 1745-2473.

[5] M. Gómez-González, E. Latorre, M. Arroyo, and

X. Trepat, Nature Reviews Physics **2**, 300 (2020), ISSN 2522-5820.

[6] S. L. Tlili, F. Graner, and H. Delanoë-Ayari, Development **149**, dev200774 (2022), ISSN 0950-1991, 1477-9129.

[7] P. Guillamat, W. Mirza, P. K. Bal, M. Gómez-González, P. Roca-Cusachs, M. Arroyo, and X. Trepat, *Guidance of cellular nematics into shape-programmable living surfaces* (2025).

[8] D. Bi, J. Lopez, J. M. Schwarz, and M. L. Manning, Nature Physics **11**, 1074 (2015).

[9] M. Popović, A. Nandi, M. Merkel, R. Etournay, S. Eaton, F. Jülicher, and G. Salbreux, New Journal of Physics **19**, 033006 (2017), ISSN 1367-2630.

[10] M. Merkel, K. Baumgarten, B. P. Tighe, and M. L. Man-
ning, Proceedings of the National Academy of Sciences **116**, 6560 (2019), ISSN 0027-8424, 1809.01586.

[11] A. Bonfanti, J. Fouchard, N. Khalilgharibi, G. Charras, and A. Kabla, Royal Society Open Science **7** (2020), ISSN 20545703.

[12] R. Etournay, M. Popović, M. Merkel, A. Nandi, C. Blasse, B. Aigouy, H. Brandl, G. Myers, G. Salbreux, F. Jülicher, et al., eLife **4**, e07090 (2015), ISSN 2050084X.

[13] N. I. Petridou, B. Corominas-Murtra, C.-p. Heisenberg, and E. Hannezo, Cell pp. 1–15 (2021), ISSN 00928674.

[14] M. Serra, G. Serrano Nájera, M. Chuai, A. M. Plum, S. Santhosh, V. Spandan, C. J. Weijer, and L. Mahadevan, Science Advances **9**, eadh8152 (2023).

[15] C. A. Dessalles, N. Cuny, A. Boutillon, P. F. Salipante, A. Babataheri, A. I. Barakat, and G. Salbreux, Nature Physics pp. 1–10 (2025), ISSN 1745-2481.

[16] M. Merkel and M. L. Manning, Seminars in Cell & Developmental Biology **67**, 161 (2017), ISSN 10849521.

[17] D. A. Matoz-Fernandez, K. Martens, R. Sknepnek, J. L. Barrat, and S. Henkes, Soft Matter **13**, 3205 (2017), ISSN 1744-683X, 1610.09340.

[18] P. Germann, M. Marin-Riera, and J. Sharpe, Cell Systems **8**, 261 (2019), ISSN 2405-4712, 2405-4720.

[19] F. Graner and J. A. Glazier, Physical Review Letters **69**, 2013 (1992), ISSN 00319007.

[20] A. F. M. Marée, V. A. Grieneisen, and P. Hogeweg, in *Single-Cell-Based Models in Biology and Medicine*, edited by A. R. A. Anderson, M. A. J. Chaplain, and K. A. Rejniak (Birkhäuser, Basel, 2007), pp. 107–136, ISBN 978-3-7643-8123-3.

[21] H. Honda, R. Kodama, T. Takeuchi, H. Yamanaka, K. Watanabe, and G. Eguchi, Journal of embryology and experimental morphology **83 Suppl**, 313 (1984), ISSN 0022-0752.

[22] R. Farhadifar, J. C. Röper, B. Aigouy, S. Eaton, and F. Jülicher, Current Biology **17**, 2095 (2007), ISSN 09609822, 1112.5905v1.

[23] S. Alt, P. Ganguly, and G. Salbreux, Philosophical Transactions of the Royal Society B: Biological Sciences **372**, 20150520 (2017), ISSN 0962-8436.

[24] R. Mueller, J. M. Yeomans, and A. Doostmohamadi, Physical Review Letters **122**, 48004 (2019), ISSN 10797114, 1811.05040.

[25] B. Loewe, M. Chiang, D. Marenduzzo, and M. C. Marchetti, Physical review letters **125**, 038003 (2020), ISSN 10797114, 1912.10549.

[26] C. Duclut, J. Paijmans, M. M. Inamdar, C. D. Modes, and F. Jülicher, The European Physical Journal E **45**, 29 (2022), ISSN 1292-895X.

[27] O. K. Damavandi, S. Arzash, E. Lawson-Keister, and M. L. Manning, PRX Life **3**, 033001 (2025).

[28] S. Tong, N. K. Singh, R. Sknepnek, and A. Košmrlj, PLOS Computational Biology **18**, e1010135 (2022), ISSN 1553-7358.

[29] S. Tong, R. Sknepnek, and A. Košmrlj, Physical Review Research **5**, 013143 (2023).

[30] S. Henkes, K. Kostanjevec, J. M. Collinson, R. Sknepnek, and E. Bertin, Nature Communications **11**, 1 (2020), ISSN 20411723, 1901.04763.

[31] D. B. Staple, R. Farhadifar, J.-C. Röper, B. Aigouy, S. Eaton, and F. Jülicher, The European physical journal. E, Soft matter **33**, 117 (2010), ISSN 1292-895X, 21082210.

[32] M. Moshe, M. J. Bowick, and M. C. Marchetti, Physi-

cal Review Letters **120**, 268105 (2017), ISSN 0031-9007, 1708.07848.

[33] D. M. Sussman and M. Merkel, Soft Matter **14**, 3397 (2018), ISSN 1744-683X, 1708.03396.

[34] M. Merkel and M. L. Manning, New Journal of Physics **20**, 022002 (2018), ISSN 1367-2630, 1706.02656.

[35] X. Wang, M. Merkel, L. B. Sutter, G. Erdemci-Tandogan, M. L. Manning, and K. E. Kasza, Proceedings of the National Academy of Sciences **117**, 13541 (2020), ISSN 0027-8424, 2005.07283.

[36] J. Huang, J. O. Cochran, S. M. Fielding, M. C. Marchetti, and D. Bi, Physical Review Letters **128**, 178001 (2022), ISSN 10797114, 2109.10374.

[37] A. Hernandez, M. F. Staddon, M. J. Bowick, M. C. Marchetti, and M. Moshe, Physical Review E **105**, 064611 (2022).

[38] C.-T. Lee and M. Merkel, *Generic elasticity of thermal, under-constrained systems* (2023), 2304.07266.

[39] C.-T. Lee and M. Merkel, *Partition sum of thermal, under-constrained systems* (2023), 2304.07264.

[40] K. Kim, T. Zhang, and J. M. Schwarz, New Journal of Physics **26**, 043009 (2024), ISSN 1367-2630.

[41] S. Kim, R. Amini, S.-T. Yen, P. Pospíšil, A. Boutillon, I. A. Deniz, and O. Campàs, Nature Materials **23**, 1592 (2024), ISSN 1476-4660.

[42] S. Gsell, S. Tlili, M. Merkel, and P.-F. Lenne, Nature Physics **21**, 644 (2025), ISSN 1745-2481.

[43] S. Okuda, Y. Inoue, M. Eiraku, T. Adachi, and Y. Sasai, Biomechanics and Modeling in Mechanobiology **14**, 413 (2015), ISSN 16177940.

[44] S. Tlili, J. Yin, J.-F. Rupprecht, M. A. Mendieta-Serrano, G. Weissbart, N. Verma, X. Teng, Y. Toyama, J. Prost, and T. E. Saunders, Proceedings of the National Academy of Sciences of the United States of America **116**, 25430 (2019).

[45] J. Rozman, K. V. S. Chaithanya, J. M. Yeomans, and R. Sknepnek, Nature Communications **16**, 530 (2025), ISSN 2041-1723.

[46] H. Koyama and Y. Inoue, Biophysics and Physicobiology p. e230007 (2026), ISSN 2189-4779.

[47] G. W. Brodland, D. Viens, and J. H. Veldhuis, Computer Methods in Biomechanics and Biomedical Engineering **10**, 121 (2007), ISSN 1025-5842.

[48] D. B. Staple, Ph.D. thesis, Technische Universität Dresden (2012).

[49] A. Nestor-Bergmann, E. Johns, S. Woolner, and O. E. Jensen, Physical Review E **97**, 052409 (2018), ISSN 2470-0045, 1711.02909v2.

[50] C. Fu, F. Dilasser, S.-Z. Lin, M. Karnat, A. Arora, H. Rajendiran, H. T. Ong, N. Mui Hoon Brenda, S. W. Phow, T. Hirashima, et al., Proceedings of the National Academy of Sciences **121**, e2405560121 (2024).

[51] S.-Z. Lin, S. Tlili, and J.-F. Rupprecht, *Viscous vertex model for active epithelial tissues* (2026), 2602.13049.

[52] M. Merkel, Ph.D. thesis, Technische Universität Dresden (2014).

[53] C. Duclut, J. Paijmans, M. M. Inamdar, C. D. Modes, and F. Jülicher, Cells & Development **168**, 203746 (2021), ISSN 2667-2901.

[54] D. Bi, J. H. Lopez, J. M. Schwarz, and M. L. Manning, Soft Matter **10**, 1885 (2014), ISSN 1744-683X, 24652538.

[55] M. Popović, V. Druelle, N. A. Dye, F. Jülicher, and M. Wyart, New Journal of Physics **23**, 033004 (2021), ISSN 1367-2630.
[56] S. Kim, M. Pochitaloff, G. A. Stooke-Vaughan, and O. Campàs, Nature Physics 17, 859 (2021), ISSN 1745-2473, 1745-2481.

[57] Q. Guigue, M. Besse, R. Voituriez, A. M. Prevost, E. Wandersman, M. Merkel, and L.-L. Pontani, *Adhesion differentials control the rheology of biomimetic emulsions* (2025), 2503.18782.

[58] D. Bi, X. Yang, M. C. Marchetti, and M. L. Manning, Physical Review X 6, 021011 (2016), ISSN 2160-3308.

[59] D. M. Sussman, M. Paoluzzi, M. Cristina Marchetti, and M. Lisa Manning, EPL (Europhysics Letters) 121, 36001 (2018), ISSN 0295-5075, 1712.05758.

[60] Y.-W. Li, L. L. Y. Wei, M. Paoluzzi, and M. P. Ciamarra, Physical Review E 103, 022607 (2021), ISSN 2470-0045, 2470-0053.

[61] C. Li, M. Merkel, and D. M. Sussman, Physical Review Letters 134, 048203 (2025), ISSN 0031-9007, 1079-7114, 2409.18949.

[62] S. Pandey, S. Kolya, P. Devendran, S. Sadhukhan, T. Das, and S. K. Nandi, Soft Matter 21, 269 (2025), ISSN 1744-6848.

[63] M. Merkel, R. Etournay, M. Popović, G. Salbreux, S. Eaton, and F. Jülicher, Physical Review E 95, 032401 (2017), ISSN 24700053, 1607.00357.

[64] D. L. Barton, S. Henkes, C. J. Weijer, and R. Sknepnek, PLoS computational biology 13, e1005569 (2017).

[65] M. Krajnc, S. Dasgupta, P. Ziherl, and J. Prost, Physical Review E 98, 022409 (2018), ISSN 2470-0045, 1805.06500.

[66] R. Sknepnek, I. Djafer-Cherif, M. Chuai, C. Weijer, and S. Henkes, eLife 12, e79862 (2023), ISSN 2050-084X.

[67] S.-Z. Lin, M. Merkel, and J.-F. Rupprecht, Physical Review Letters 130, 058202 (2023).

[68] J. Rozman and J. M. Yeomans, Physical Review Letters 133, 248401 (2024).

[69] K. Barrett, S. Anand, V. Thome, L. Kodjabachian, M. Merkel, and P.-F. Lenne, Development pp. dev-204344 (2025).

[70] A. Shindo, Wiley Interdisciplinary Reviews: Developmental Biology 7, e293 (2018), ISSN 17597684.

[71] G. K. Batchelor, Journal of Fluid Mechanics 41, 545 (1970), ISSN 0022-1120.

[72] N. Noll, M. Mani, I. Heemskerk, S. J. Streichan, and B. I. Shraiman, Nature Physics 13, 1221 (2017), ISSN 1745-2473, 1508.00623.

[73] D. Grossman and J.-F. Joanny, Physical Review Letters 129, 048102 (2022).

[74] D. Grossman and J.-F. Joanny, Physical Review Research 7, 013039 (2025), ISSN 2643-1564.

[75] B. Guirao, S. U. Rigaud, F. Bosveld, A. Bailles, J. López-Gay, S. Ishihara, K. Sugimura, F. Graner, and Y. Bellaïche, eLife 4 (2015), ISSN 2050-084X.

[76] S. Tlili, M. Durande, C. Gay, B. Ladoux, F. Graner, and H. Delanoë-Ayari, Physical Review Letters 125, 88102 (2020), ISSN 10797114, 1811.05001.

[77] R. A. Simha and S. Ramaswamy, Physical review letters 89, 058101 (2002), ISSN 00319007, cond-mat/0108301v2.

[78] R. Voituriez, J. F. Joanny, and J. Prost, Europhysics Letters 70, 404 (2005), ISSN 02955075, q-bio/0503022.

[79] M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao, and R. A. Simha, Reviews of Modern Physics 85, 1143 (2013), ISSN 0034-6861, 1207.2929v1.

[80] M. Ibrahimi and M. Merkel, New Journal of Physics 25, 013022 (2023), ISSN 1367-2630, 2206.12850.

[81] M. Ibrahimi and M. Merkel, PRX Life 3, 043013 (2025), 2412.15774.