# Bridging microscopic cell dynamics to nematohydrodynamics of cell monolayers

Aleksandra Ardaševa  The Niels Bohr Institute, University of Copenhagen,
Copenhagen, Denmark  Romain Mueller  Rudolf Peierls Centre for Theoretical
Physics, University of Oxford, UK  Amin Doostmohammadi
<EMAIL_ADDRESS>The Niels Bohr Institute, University of Copenhagen,
Copenhagen, Denmark

###### Abstract

It is increasingly being realized that liquid-crystalline features can play an
important role in the properties and dynamics of cell monolayers. Here, we
present a cell-based model of cell layers, based on the phase-field
formulation, that connects mechanical properties at the single cell level to
large-scale nematic and hydrodynamic properties of the tissue. In particular,
we present a minimal formulation that reproduces the well-known bend and splay
hydrodynamic instabilities of the continuum nemato-hydrodynamic formulation of
active matter, together with an analytical description of the instability
threshold in terms of activity and elasticity of the cells. Furthermore, we
provide a quantitative characterisation and comparison of flows and
topological defects for extensile and contractile stress generation
mechanisms, and demonstrate activity-induced heterogeneity and spontaneous
formation of gaps within a confluent monolayer. Together, these results
contribute to bridging the gap between micro-scale cell dynamics and tissue-
scale collective cellular organisation.

##  I Introduction

Collective cell migration is crucial for a wide variety of biological
scenarios, ranging from development Heisenberg and Bellaïche (2013) to wound
healing Brugués  _et al._ (2014) and pathological conditions, such as cancer
metastasis Friedl  _et al._ (2012). The profound role of mechanical forces in
regulating such behaviour is being increasingly recognized in a variety of
cell types Ladoux and Mège (2017). On the other hand, these collective
behaviours are to a large degree generic and can be described by different
classes of theoretical models that target specific levels of description. In
particular, there is a well-documented connection between liquid crystal
theories and the description of epithelial cells at large scales
Doostmohammadi and Ladoux (2021); Duclos  _et al._ (2018); Kawaguchi  _et al._
(2017); Prost  _et al._ (2015); Doostmohammadi  _et al._ (2018).

Nematic liquid crystals consist of rod-like particles, which under certain
conditions align in a given direction. This gives rise to quasi long-range
nematic order while singularities in orientational alignment lead to
topological defects. Theories of nematic liquid crystals have recently found
application in different biological systems. For example, the local
orientation of the actin fibres in the tissue can be effectively described as
a nematic liquid crystal Maroudas-Sacks  _et al._ (2021), and elongated
fibroblasts at high densities show nematic alignment and topological defects
due to steric interactions Duclos  _et al._ (2018). Even though epithelial
cells have no obvious microscopic nematic degree of freedom, several flow
features of cell layers have been successfully reproduced in continuum
theories of active liquid crystals. This includes flow fields around cell
division events Doostmohammadi  _et al._ (2015), long-range vortex lattices in
dividing cell layers Rossen  _et al._ (2014); Doostmohammadi  _et al._ (2016),
chaotic flow features within the active turbulent state Blanch-Mercader  _et
al._ (2018); Lin  _et al._ (2021), and formation of topological defects, that
have been linked with death and extrusion events in monolayers of MDCK cells
Saw  _et al._ (2017). However, there is still little understanding on how
emergent flows relate to microscopic dynamics of cells.

The mechanisms by which hydrodynamic modes emerge from overdamped microscopic
dynamics is well understood for self-propelled point particles with aligning
interactions, and there are now powerful methods that allow consideration of a
large class of such models Marchetti  _et al._ (2013); Peshkov  _et al._
(2014). There have been efforts to apply similar ideas to cellular layers Lee
and Wolgemuth (2011a, b); Gov (2009) but how these models relate to a
hydrodynamic description is less clear because dense cellular systems are far
from the dilute limit and steric interactions are important. More complex
physical descriptions of dense layers have also been explored Camley and
Rappel (2017) but most of the investigations concentrated around rigidity
transition and jamming Bi  _et al._ (2016); Henkes  _et al._ (2011) or tissue
reshaping by the extracellular matrix Etournay  _et al._ (2015).

Lately, the phase-field formalism has proven successful in reproducing various
experimental phenomenologies, including mechanical stress patterns Saw  _et
al._ (2017), as well as the cell-substrate and cell-cell interaction forces in
epithelial cell layers Balasubramaniam  _et al._ (2021). In these methods, the
cells are modelled as active deformable droplets in two-dimensions whose
boundaries are defined using a phase-field potential which allows for the
precise description of the intercellular interfaces and forces. Compared to
models where the individual cells are described using a tesselation of the
plane Bi  _et al._ (2016); Barton  _et al._ (2017), it has the advantages that
it can accommodate large shape deformations of the cells independently of the
location of their center-of-masses and that it naturally allows for the
description of tissue boundaries, including gaps within the tissue.

The phase-field approach has been used to study the behaviour of cellular
systems ranging from single cells Aranson (2016) to multicellular structures
Najem and Grant (2016) and cell monolayers Basan  _et al._ (2013); Sepúlveda
_et al._ (2013); Monfared  _et al._ (2021); Mueller  _et al._ (2019); Zhang
_et al._ (2020). However, most of these models introduce polar driving of the
individual cells as a source of activity and/or polar interactions between
cells, and it is unclear if such microscopic dynamics can be properly modelled
by a hydrodynamic continuum model of nematic liquid crystals at macroscopic
scales. A formulation based on coupling active stresses to cell shape
deformations showed the possibility of capturing some aspects of topological
defect formation in isotropic cell layers Mueller  _et al._ (2019). Recently,
this approach has been extended to account for both extensile and contractile
stresses, noting how these stresses contribute to cell shape deformations
Zhang and Yeomans (2021). However, the mechanisms of topological defect
formation and the link with the well-established instabilities in continuum
formulations of active nematics remain unresolved.

In this paper, we present a simple microscopic model of a dense 2D cellular
layer that reproduces the phenomenology of hydrodynamic theories of active
nematic liquid crystals and introduce a generic model of force transmission at
the cellular interfaces, which allows us to connect tissue-level quantities to
the microscopic description. We then use insights from the continuum theory of
active nematic liquid crystals to introduce a minimal form of local dipolar
interaction between cells and show that our model develops the correct
hydrodynamic behaviour at macroscopic scales. In particular, we show that this
formulation reproduces the typical bend and splay instabilities of active
nematic liquid crystals theories and that it produces topological defects and
associated flow features that are similar to the ones observed experimentally.
Furthermore, within a certain parameter regime, we show that the competition
between active stress generation and the elasticity of individual cells
results in the spontaneous formation of gaps within the confluent layers.

The paper is organized as follows: In Section II we describe the details of
the model and the minimal description of the cell-generated forces and
intercellular interactions. In Section III.1 we present numerical simulations
of the system that reproduce fundamental hydrodynamic instabilities of
continuum active nematics, followed by a linear stability analysis in Section
III.2 that sheds light on the mechanism of the generation of such hydrodynamic
modes within our cell-based description. In Section III.3 we employ numerical
simulations to move beyond these instabilities and analyse more closely the
dynamics of topological defects, followed by characterization of the flow
features of active turbulence in the cell-based model in Section III.4.
Finally, we present quantitative results of gap formation in extensile and
contractile systems in Section III.5, followed by concluding remarks and
discussions in Section IV.

##  II The phase-field model

Following Mueller et al. Mueller  _et al._ (2019), we represent each cell in
an epithelium monolayer via a separate phase field, $\phi_{i}$. The dynamics
of each phase field is governed by

$\partial_{t}\phi_{i}+\vec{v}_{i}\cdot\vec{\nabla}\phi_{i}=-\frac{\delta\mathcal{F}}{\delta\phi_{i}},\qquad
i=1,\ldots,N,$ (1)

where $N$ is the total number of cells and $\vec{v}_{i}$ is the total velocity
of the cell $i$. The free energy, $\mathcal{F}$, governs the dynamics of the
cell interfaces and has three contributions:

$\mathcal{F}=\mathcal{F}_{\text{CH}}+\mathcal{F}_{\text{area}}+\mathcal{F}_{\text{rep.}}.$
(2)

The first term in eqn (2) is the Cahn-Hilliard free energy, which stabilizes
the interface between exterior and interior of each cell. It is defined as

$\mathcal{F}_{\text{CH}}=\sum_{i}\frac{\gamma}{\lambda}\int\textup{d}\vec{x}\left\\{4\phi_{i}^{2}(1-\phi_{i})^{2}+\lambda^{2}(\vec{\nabla}\phi_{i})^{2}\right\\},$
(3)

where $\lambda$ is the interface width and $\gamma$ is the surface tension
that characterizes the cell elasticity Palmieri  _et al._ (2015). The double-
well potential in the first term of eqn (3) makes $\phi_{i}=1$ for regions
inside the cells and $\phi_{i}=0$ for the outside regions. The cell boundary
occurs at $\phi_{i}=0.5$. This formulation is guided by simplicity and models
the cell boundary as a diffuse interface, but could be easily extended to
model sharp interfaces Helfrich (1973); Shao  _et al._ (2010); Biben and
Misbah (2003).

The second term in eqn (2), $\mathcal{F}_{\text{area}}$, defines a soft area
constraint for each cell and is given by

$\mathcal{F}_{\text{area}}=\sum_{i}\mu\Big{(}1-\frac{1}{\pi
R^{2}}\int\textup{d}\vec{x}\,\phi_{i}^{2}\Big{)}^{2},$ (4)

where $\mu$ controls the relaxation of area changes and $R$ is the cell
radius. $\mathcal{F}_{\text{area}}$ ensures that cell areas,
$A_{i}=\int\textup{d}\vec{x}\phi_{i}^{2}$, are close to circle, i.e. $\pi
R^{2}$. Note that area conservation is not implemented at the level of eqn
(1), i.e. even though cells are mostly incompressible, the cell area can be
dramatically altered, for example, squeezed by neighbouring cells or stretched
on the substrate Peyret  _et al._ (2019).

Finally, the repulsion term, $\mathcal{F}_{\text{rep.}}$, penalizes regions
where two cells overlap:

$\mathcal{F}_{\text{rep.}}=\sum_{i}\sum_{j\neq
i}\frac{\kappa}{\lambda}\int\textup{d}\vec{x}\,\phi_{i}^{2}\phi_{j}^{2}\,.$
(5)

Here, $\kappa$ denotes the strength of the repulsion. Normalisation has been
chosen such that the width of the interfaces at equilibrium is given by
$\lambda$ and that the properties of the cells are roughly preserved when
$\lambda$ is rescaled.

The above mentioned free energies describe the cellular interfaces, as shown
in Fig. 1a, and specify the equilibrium properties of the interfaces including
surface tension and interface width. To drive the system away from the
equilibrium we define the cell velocities $\vec{v}_{i}$ in eqn (1) based on
the force balance for an overdamped system:

$\xi\vec{v}_{i}=\vec{F}^{\text{\text{int.}}}_{i},$ (6)

where $\xi$ is the substrate friction coefficient and
$\vec{F}^{\text{\text{int.}}}_{i}$ is the total force acting on the interface
of the cell $i$. When modelling epithelia from a material science perspective,
these microscopic interface forces are usually defined in terms of a
macroscopic tissue stress tensor, $\bm{\sigma}_{\text{tissue}}$. In our phase-
field formulation, this corresponds to

$\vec{F}^{\text{\text{int.}}}_{i}=\int\mathrm{d}\vec{x}\,\phi_{i}\;\vec{\nabla}\cdot\bm{\sigma}_{\text{tissue}}=-\int\mathrm{d}\vec{x}\,\bm{\sigma}_{\text{tissue}}\cdot\vec{\nabla}\phi_{i}.$
(7)

These expressions can be interpreted as follows: the first is the integral of
the local force, $\vec{\nabla}\cdot\bm{\sigma}_{\text{tissue}}$, weighted by
the phase-field, $\phi_{i}$, and the second is the integral of the force
exerted by the stress tensor on the vector $-\vec{\nabla}\phi_{i}$ normal to
the interface. The eqn (7) connects the local forces at the level of the
individual cells with the properties of the cell monolayer. In the limit of
sharp interfaces, the integral tends to a contour integral over the cell
boundary and leads to the usual expressions for the local force in terms of
the stress tensor, see Section III.2 below.

Figure 1: (a) Interfaces between cells defined as the overlap $\sum_{i\neq
j}\phi_{i}\phi_{j}$ and (b) contours $\phi_{i}=1/2$ with an illustration of
the contractile force exerted by the cell in the centre on its neighbours.

In order to reproduce the phenomenology of active nematic liquid crystals, we
include local dipolar interactions into $\bm{\sigma}_{\text{tissue}}$. To do
so, we associate internal nematic tensor, $\bm{Q}_{i}$, to each cell
describing the orientational alignment of the cell. The tissue stress tensor
is defined as follows,

$\bm{\sigma}_{\text{tissue}}=-P\bm{I}-\zeta\bm{Q},$ (8)

where $\bm{Q}=\sum_{i}\bm{Q}_{i}\phi_{i}$ is tissue’s nematic field,
$P=\sum_{i}\partial\mathcal{F}/\partial\phi_{i}$ is the pressure field, and
$\bm{I}$ is the identity tensor. The first term in eqn (8) induces steric
repulsion between cells, while the second term mirrors the active term found
in continuum theories of active liquid crystals Marchetti  _et al._ (2013). It
can be interpreted as a dipolar force density distributed along the cells
interfaces such that each cell pushes or pulls its neighbours depending on the
direction of their contact area with respect to the nematic tensor, as
demonstrated schematically in Fig. 1b.

We now introduce the dynamics of the internal nematic degree of freedom,
$\bm{Q}_{i}$, of each cell. In the continuum theory, most of the phenomenology
of active liquid crystals depends on a balance between the restoring elastic
forces and flow alignment of the nematic tensor. By writing
$\bm{Q}_{i}=2\left(\vec{n}_{i}\otimes\vec{n}_{i}-\vec{n}_{i}^{2}\,\bm{I}/2\right)$
with $\vec{n}_{i}=(\cos\theta_{i},\sin\theta_{i})$, we can mirror these two
components in our model and define the following dynamics of the angle,
$\theta_{i}$:

$\partial_{t}\theta_{i}=K\tau_{i}+J\omega_{i},$ (9)

where the torques are given by

$\displaystyle\tau_{i}=\frac{1}{\lambda}\int\textup{d}\vec{x}\;\bm{Q}_{i}\phi_{i}\wedge\bm{Q},\quad\omega_{i}=\int\textup{d}\vec{x}\;\vec{v}\wedge\vec{\nabla}\phi_{i},$
(10)

and $K$ is the nematic elastic parameter and $J$ is the nematic flow alignment
strength. We have defined the tissue velocity as
$\vec{v}=\sum_{i}\vec{v}_{i}\phi_{i}$ and
$\bm{A}\wedge\bm{B}=A_{xx}B_{xy}-A_{xy}B_{xx}$ for symmetric traceless
matrices $\bm{A}$ and $\bm{B}$, while $\vec{v}\wedge\vec{\nabla}\phi_{i}$ is
the usual wedge product. The first torque, $\tau_{i}$, aligns $\bm{Q}_{i}$ to
the tissue nematic tensor, $\bm{Q}$, and induces an elastic restoring force
favouring the homogeneous state. The second torque,
$\bm{\sigma}_{\text{tissue}}$, rotates $\bm{Q}_{i}$ with the local vorticity
computed as the integral of the neighbouring cells velocity projected on the
interface. Note that our definitions are local in the sense that each torque
is given by an integral along the cell interface with its neighbours.

We simulate eqns (1) and (9) using an open-source software Celadro (available
on GitHub: https://github.com/rhomu/celadro), which utilizes a finite
difference scheme on a square lattice with a predictor-corrector step (for
details of the implementation, see Mueller et al. Mueller  _et al._ (2019)).
For all the simulations presented here, we initialize each cell with a
director along the $x$-direction. We simulate the dynamics of 896 cells of
radius $R=8$ in a $400\times 400$ periodic domain. Cells are placed randomly
and relaxed for 200 time steps under passive dynamics to reach confluence.
Simulations are run for 15000 time steps with the following
parameters/parameter ranges (unless stated otherwise): $\mu=0.03$,
$\lambda=8$, $\kappa=0.2$, $\xi=1$, $K=2\times 10^{-5}$, $J=4\times 10^{-3}$,
$\zeta\in\pm[3\times 10^{-4},8\times 10^{-2}]$, $\gamma\in[0.01,0.1]$, chosen
within the range that was previously shown to reproduce defect flow fields in
epithelial layers Balasubramaniam  _et al._ (2021).

##  III Results

To demonstrate the new features of the developed model and highlight the
consistency with the continuum model, we begin by numerically investigating
the mechanism and dynamics of the system starting from a quiescent state. In
our analysis, we explore the behaviour of the model for a range of activity
parameter values, $\zeta$, (both positive and negative) and consider different
values of the cell elasticity $\gamma$.

###  III.1 Bend and splay instabilities

Continuum theory shows that when a nematic system is active, it is prone to
bend or splay instability, which depends on the sign of the activity parameter
Simha and Ramaswamy (2002); Ramaswamy and Rao (2007); Thampi  _et al._ (2014).
Let us define bend and splay measure as follows:

$b=\int\mathrm{d}\vec{x}\,(\vec{n}\wedge\vec{\nabla}\wedge\vec{n})^{2}/2,\,\,\,s=\int\mathrm{d}\vec{x}\,(\vec{\nabla}\cdot\vec{n})^{2}/2.$
(11)

Initializing a system from a quiescent state, we observe that after some
transient relaxation time, the tissue nematic field, $\bm{Q}$, develops non-
zero bend (splay) when the cell-cell interaction is extensile (contractile),
as demonstrated in Fig. 2a,b. This initial instability then leads to the
formation of walls (Fig. 2a) – elongated distortions in the director field –
which are also predicted by the continuum theory Ramaswamy and Rao (2007);
Giomi  _et al._ (2011); Marchetti  _et al._ (2013).

Figure 2: Bend-splay instabilities of a cell monolayer. (a) Tissue nematic
field, $\bm{Q}=\sum_{i}\bm{Q}_{i}\phi_{i}$, showing typical bend (left) and
splay (right) instabilities of an initially homogeneous state for extensile
($\zeta>0$) and contractile ($\zeta<0$) activities, respectively. Parameters
used here: $\zeta=\pm 0.01$, $\gamma=0.05$. (b) Evolution of the total bend,
$b$, and splay, $s$, given by eqn (11) of the tissue nematic field, $\bm{Q}$,
as a function of time for extensile (top) and contractile activities (bottom).
Averaged over 5 simulations with $\zeta=\pm 0.02$ and $\gamma=0.05$. (c)
Number of time steps (rescaled by $1/RK$) before creation of the first defect
pair for different extensile (top) and contractile activities (bottom). The
dashed line is the exponential fit. Mean$\pm$std from 5 simulations with
$\gamma=0.1$. (d) Snapshots of simulation demonstrating unzipping of walls by
topological defects during bend instability for an extensile activity. Green
circle: $+1/2$-defect, blue triangle: $-1/2$-defect. The simulation was
performed with $\zeta=0.01$ and $\gamma=0.05$.

At longer times, we observe the spontaneous creation of pairs of topological
defects that lead to the unzipping of walls in the orientation field for both
bend and splay instabilities, see Fig. 2d. In order to quantify the timescale
associated with the formation of these instabilities, we measure the time
before the creation of the first defect pair for different values of the
activity, $\zeta$ (Fig. 2c). We observe that as the activity increases, the
time to the formation of first defect pair decreases exponentially, which is
consistent with the predictions from continuum theory on the exponential
growth of the instabilities in time with the time constant proportional to
$1/|\zeta|$, where $\zeta$ is the activity coefficient Giomi  _et al._ (2011);
Thampi  _et al._ (2014).

###  III.2 Analysis of the instabilities at linear order

In this section, we show how the dynamics presented in eqns (9) and (10)
capture the bend and splay instabilities observed in continuum systems of
active nematics by performing a linear stability analysis. To keep the
description analytically tractable, we make several assumptions that allow us
to derive a continuum approximation of our multi-agent model in the double
limit where both interfaces and cells are small. We then consider the
evolution of this continuum description from a quiescent state at linear order
and obtain stability conditions for the system.

Our procedure is as follows. First, note that in the absence of activity, i.e.
$\zeta=0$, the epithelium relaxes to a honeycomb lattice where all cells take
an approximate hexagonal form (Fig. 3).

Figure 3: A cell and its neighbours in a honeycomb configuration showing
approximate hexagonal shape. Additional surrounding cells are not shown. The
locations of the centre of the first cell $\textbf{{x}}_{0}$ and of two of its
neighbours $\vec{x}_{i}=\vec{x}_{0}+\delta\vec{m}_{i-1}$ are indicated.

Consider then a cell with phase-field $\phi_{0}$ centred at $\vec{x}$ with
neighbouring cells located at $\vec{x}_{i}=\vec{x}_{0}+\delta\vec{m}_{i-1}$
for $i=1,\ldots,6$, where $\delta$ is the distance between the centres of two
neighbouring cells, and
$\vec{m}_{i}=\left(\cos\left(i\pi/3\right),\sin\left(i\pi/3\right)\right)$. We
then assume that the widths of the cells interfaces are much smaller than the
cells sizes, i.e. $\lambda\ll\delta$, such that the cells take an hexagonal
form with $\delta\approx\sqrt{3}R$. In this limit, all the integrals over the
phase fields can easily be performed by separating the directions
perpendicular and parallel to the interfaces. Assuming constant field along
the parallel direction and equilibrium interface profiles
$\phi_{\pm}(\vec{x}-\vec{x}_{0})=\frac{1}{2}\left(1\pm\tanh\left(\vec{x}/l\right)\right)$
along the perpendicular direction, we can use the explicit integrals

$\int\mathrm{d}\vec{x}\,\phi_{\pm}(\vec{x})\phi_{\mp}(\vec{x})=l/2,\quad\int\mathrm{d}\vec{x}\,\phi_{\pm}(\vec{x})\partial_{x}\phi_{\mp}(\vec{x})=\mp
1/2,$

to obtain the following expressions

$\int\mathrm{d}\vec{x}\,\phi_{i}\phi_{0}\approx\frac{\delta
l}{2\sqrt{3}},\qquad\text{and}\quad\int\mathrm{d}\vec{x}\,\phi_{i}\vec{\nabla}\phi_{0}\approx-\frac{\delta\vec{m}_{i}}{2\sqrt{3}},$

where $\delta/\sqrt{3}$ is the length of the interface between cells $0$ and
$i$. Additionally, we have used the fact that the gradients are perpendicular
to the interfaces, i.e.
$\phi_{i}\vec{\nabla}\phi_{0}\approx\phi_{i}\vec{m}_{i}(\vec{m}_{i}\cdot\vec{\nabla})\phi_{0}$.
Note that even though we assume that the interfaces take their equilibrium
profiles, the interface length, $l$, will be in principle different from
$\lambda$ because we consider a confluent state where the steric interactions
between cells cannot be ignored. Once such integrals over the individual cells
have been computed, we can consider the limit where cells are small, i.e.
$\lambda\ll\delta\ll 1$. We rescale all the physical parameters that multiply
terms proportional to integrals over the cells interfaces by the interface
length, i.e. we replace $K\to K/\delta$, $J\to J/\delta$, and
$\xi\to\xi\delta$, and rewrite $\bm{Q}_{i}$ and $\vec{v}_{i}$ as fields.
Performing an expansion to first order in $\delta$ finally allows us to obtain
expressions for the tissue level fields in the continuum limit and perform a
linear stability analysis.

We first apply this procedure to better understand the effect of the torques
in eqn (10). Considering the cell $i=0$ located at $\vec{x}$, we can write

$\displaystyle\tau_{0}$
$\displaystyle=\frac{1}{\lambda}\int\mathrm{d}\vec{x}\,\bm{Q}_{0}\phi_{0}\wedge\bm{Q}\approx\frac{l}{\lambda}\frac{\delta}{2\sqrt{3}}\sum_{i=1}^{6}\bm{Q}_{0}\wedge\bm{Q}_{i}$
$\displaystyle=\frac{l}{\lambda}\frac{\delta}{2\sqrt{3}}\sum_{i=1}^{6}\bm{Q}(\vec{x})\wedge\bm{Q}(\vec{x}+\delta\vec{m}_{i})$
$\displaystyle=\frac{l}{\lambda}\frac{\sqrt{3}}{2}\Delta\theta(\vec{x})\delta^{3}+\mathcal{O}(\delta^{5}),$

where we have set $Q_{xx}=\cos 2\theta$ and $Q_{xy}=\sin 2\theta$. This form
is similar to the elastic term of continuum theories of liquid crystals in the
one-elastic constant approximation, and shows that the torque $\tau_{0}$ does
indeed have the effect of an elastic restoring force. We can similarly rewrite
the vorticity torque as

$\displaystyle\omega_{0}$
$\displaystyle=\int\mathrm{d}\vec{x}\,\vec{v}\wedge\vec{\nabla}\phi_{0}=-\frac{\delta}{2\sqrt{3}}\sum_{i=1}^{6}\vec{v}_{i}\wedge\vec{m}_{i}$
$\displaystyle\approx-\frac{\delta}{2\sqrt{3}}\sum_{i=1}^{6}\vec{v}(\vec{x}+\delta\vec{m}_{i})\wedge\vec{m}_{i}$
$\displaystyle=\frac{\sqrt{3}}{2}(\partial_{x}v_{y}(\vec{x})-\partial_{y}v_{x}(\vec{x})\delta^{2}+\mathcal{O}(\delta^{4}).$

This shows that the torque $\omega_{i}$ indeed measures the local tissue
vorticity around a cell. Finally, we consider the total force on the interface
and, neglecting the pressure term, obtain the following expansion

$\displaystyle\vec{F}_{0}^{\text{int.}}$
$\displaystyle=\zeta\int\mathrm{d}\vec{x}\,\bm{Q}\cdot\vec{\nabla}\phi_{0}\approx-\frac{\zeta\delta}{2\sqrt{3}}\sum_{i=1}^{6}\bm{Q}_{i}\cdot\vec{m}_{i}$
$\displaystyle=-\frac{\zeta\delta}{2\sqrt{3}}\sum_{i=1}^{6}\bm{Q}(\vec{x}+\delta\vec{m}_{i})\cdot\vec{m}_{i}$
$\displaystyle=-\frac{\sqrt{3}\zeta}{2}\,\begin{pmatrix}\partial_{x}Q_{xx}+\partial_{y}Q_{xy}\\\
\partial_{x}Q_{xy}-\partial_{y}Q_{xx}\end{pmatrix}\delta^{2}+\mathcal{O}(\delta^{4}),$

which is similar to the local force arising from the active stress term in
continuum theories of liquid crystals.

We can now obtain the linear stability condition for the homogeneous state
$\theta=0$ by using the force balance equation
$\vec{v}=\vec{F}_{0}^{\text{int.}}/\xi$ and expand eqn (9) to linear order in
$\theta$. This gives us

$\partial_{t}\theta=\delta^{2}\left[\frac{l}{\lambda}\frac{\sqrt{3}K}{2}\Delta+\frac{3J\zeta}{2\xi}(\partial_{x}^{2}-\partial_{y}^{2})\right]\theta+\mathcal{O}(\delta^{4},\theta^{2}),$

where we have neglected terms of order $\delta^{4}$ and $\theta^{2}$. For
extensile (contractile) systems, it is easy to consider perturbations
homogeneous under translation in the $y$ direction ($x$ direction) and obtain
the linear stability condition, $-\zeta_{c}<\zeta<\zeta_{c}$, where the
critical activity is given by

$\zeta_{c}=\frac{l}{\lambda}\frac{\xi K}{3RJ}\propto\frac{\xi K}{J}.$

Hence, we see that the system shows bend and splay instabilities when
$\zeta>\zeta_{c}$ and $\zeta<\zeta_{c}$, respectively.

We examine the predictions of the linear stability analysis by performing
simulations for different values of the activity parameter, $\zeta$, of the
nematic elastic parameter, $K$, and of the nematic flow alignment strength,
$J$. For each of the simulations we measure the average bend and splay defined
by eqn. (11) over the last 5000 timesteps. Fig. 4a shows the average bend and
splay as a function of activity for fixed $K$ and $J$, confirming the
existence of an activity threshold. It is apparent that beyond a certain
activity value, the system develops non-zero bend and splay as predicted by
linear stability analysis, both in the extensile and contractile case.

Figure 4:  Activity threshold for the instabilities. (a) Bend and splay (eqn
(11)) after 10000 steps for different values of activity $\zeta$, with
$\gamma=0.05$, $K=2\times 10^{-5}$ and $J=4\times 10^{-3}$. Mean$\pm$std over
5 simulations. (b) Dependence of the critical activity, $\zeta_{C}$, on the
nematic elastic parameter, $K$, (for $J=4\times 10^{-3}$) and nematic flow
alignment strength, $J$, (with $K=2\times 10^{-5}$) for $\gamma=0.05$. Dashed
lines represent scaling: $\zeta_{C}\propto K$ and $\zeta_{C}\propto 1/J$, as
predicted from the linear stability analyses.

We then estimate the critical activity $\zeta_{C}$ from simulations by finding
the first value of $\zeta$ for which bend deviates substantially, i.e. $b\geq
100$, from zero at the end of the simulation (last 5000 time steps) and
examine how this critical value depends on the elastic constant $K$ and flow-
alignment constant $J$. Fig. 4b shows the critical activity $\zeta_{C}$ as a
function of $K$ (top) and $J$ (bottom). The dashed lines demonstrates the
expected scaling as predicted from the linear stability analyses, i.e.,
$\zeta_{C}\propto K$ and $\zeta_{C}\propto 1/J$. Together, these results
establish the existence of a well-defined critical activity, above which a
collection of active deformable cells exhibit bend/splay instabilities.

Figure 5: Topological defect statistics. (a) Rescaled rate of creation of
defect pairs (top) and density of defects (bottom) as a function of the
activity, $|\zeta|$, for different values of the elastic constant, $\gamma$.
Solid line – extensile, dotted – contractile (dashed lines not visible since
extensile and contractile curves overlap). Mean$\pm$std from 5 simulations.
(b) Defect speed divided by mean root-mean-square velocity of the entire
monolayer for positive (red) and negative (blue) defects for extensile (solid)
and contractile (dotted) activities. Here, $\gamma=0.1$. Mean$\pm$std from 5
simulations. (c) Mean-squared-distance (rescaled by squared cell radius,
$R^{2}$) travelled by defect from its creation point as a function of time for
positive (red) and negative (blue) defects with $\zeta=\pm 0.01$ (left) and
$\zeta=\pm 0.06$ (right), and $\gamma=0.1$. Top row: extensile, bottom row:
contractile system. The time is rescaled by $1/RK$. Mean$\pm$std from all
defects in one simulation. (d) Diffusion exponent for $\gamma=0.1$ as a
function of activity for extensile (solid) and contractile (dotted) defects.
Mean$\pm$std from 5 simulations. (e) Average properties of $+1/2$-defects
(top) and $-1/2$-defects (bottom) in an extensile and contractile systems with
$\zeta=\pm 0.04$ and $\gamma=0.1$. Average over time and over 5 simulations.

###  III.3 Active turbulence: defects characteristics

As expected from continuum theory of liquid crystals, the hydrodynamic
instability of the ordered state and the subsequent defect pair nucleation is
followed, at longer times and for high enough activity, by the emergence of
active turbulence Giomi  _et al._ (2013). The chaotic flows of active cellular
systems are often interleaved with dynamic organisation of topological defects
Saw  _et al._ (2017); Kawaguchi  _et al._ (2017); Blanch-Mercader  _et al._
(2018). As such, we next look at the rate of defect pair creation and the
density of defects for different values of activity, $\zeta$, and elasticity,
$\gamma$, shown in Figure 5a. In line with the predictions of continuum
theories of active nematics Giomi  _et al._ (2013); Thampi  _et al._ (2014),
both values of rate and density increase for larger activities in both
contractile and extensile systems until they saturate. The rate of saturation
depends on the elasticity, i.e. larger values of $\gamma$ lead to more rapid
increase in defect nucleation rate and a faster saturation of the defect
density.

It is well established that although topological defects in active nematics
appear and annihilate in pairs with positive and negative half-integer charge,
the defects with different charge show distinct dynamic behaviour. In
particular, while negative defects are three-fold symmetric and lack self-
propulsion, the positive, comet-shaped, defects in active systems are endowed
with self-propulsive velocity along their comet-head (-tail) for extensile
(contractile) activities.

To characterize the dynamics of positive and negative half-integer defects we
first measured the average speed of defects normalised by root-mean-squared
velocity of the flow in the entire domain, which demonstrate higher velocities
of positive defects at low activities (Fig. 5b). At low activities, both
positive and negative half-integer defects show larger velocities in extensile
systems compared to the contractile one. At higher activities, however, the
normalised defect velocities, for both positive and negative half-integer
defects and in both extensile and contractile systems, drop to the same level.
To study this trend more closely, we calculate the mean squared distance
travelled from the time of defect creation for each defect as a function of
time (Fig. 5c) and obtain the corresponding diffusion exponent at long times
(Fig. 5d). When the activity is relatively low, the defect density is lower
and positive defects show larger diffusion exponents, indicating the impact of
self-propulsion on their dynamics. However, when the activity is increased and
the density of the defects approaches its saturation values (Fig. 5a), the
diffusion constants of both positive and negative defects become equal and
tend to $1$. Similar behaviour is observed for different values of elasticity,
$\gamma$, and remains the same in both extensile and contractile systems. This
is in contrast to the predictions of continuum theories Thampi  _et al._
(2014) and observations in bacterial colonies Meacock  _et al._ (2021) and
subcellular filaments-motor protein mixtures Sanchez  _et al._ (2012); Kumar
_et al._ (2018), where positive and negative half-integer defects show
distinct propulsive and diffusive behaviours, respectively. Remarkably,
however, the diffusive behaviour for both defects, observed in our simulations
at high activities, is in agreement with the experimental observations on
Human Bronchial Cells (HBC) Blanch-Mercader  _et al._ (2018) and recent phase-
field modelling with distinct mechanisms of active driving Wenzel and Voigt
(2021), suggesting that such diffusive defect dynamics are generic features of
deformable active particles at high activities, where the nematic order is an
emergent property and defects motion is dominated by the interaction between
defects at long times.

Figure 6:  Global flow features of the active turbulence in the monolayer. (a)
Rescaled root-mean-square velocity, $\sqrt{v_{x}^{2}+v_{y}^{2}}$, (solid) and
sum of both components, $v_{x}+v_{y}$, (dashed) plotted for different values
of activity, $\zeta$, and elasticity, $\gamma$. Averaged over the whole domain
and time (5 realisations). (b) Dependence of the vorticity length defined as
the location of the minimum of the vorticity autocorrelation function as a
function of $\zeta$ for different values of $\gamma$. Mean$\pm$std from 5
simulations. (c) Distribution of vortex areas, defined by Okubo-Weiss
parameter, for different values of activities, $\zeta$, with $\gamma=0.1$.
Averaged over 5 realisations. The $x$-axis is rescaled by $\pi R^{2}$ and the
black solid-line marks a possible exponential scaling at intermediate vortex
sizes. (d) Absolute value of the divergence as a function of activity,
$\zeta$, rescaled by $v_{\text{rms}}/R$. Averaged over 5 simulations.

To further complement the analyses of the defect dynamics and highlight the
structural differences between the different defect types, we also calculated
the average flows around positive and negative defects since the flow patterns
around defects are crucial in determining the dynamics of active turbulence
Giomi  _et al._ (2013) and have been recently shown to play a role in the
mechanical properties of epithelial monolayers Saw  _et al._ (2017); Kawaguchi
_et al._ (2017); Doostmohammadi and Ladoux (2021). We extract the average
properties of defects in the case of extensile interaction $\zeta>0$ and
obtain the local properties for positive defects (top row) and for negative
defects (bottom row), as shown in Fig. 5e for both contractile and extensile
defects. They reproduce faithfully the theoretical predictions from continuum
models as well as the experimental observations made in epithelial monolayers
Saw  _et al._ (2017).

###  III.4 Active turbulence: flow characteristics

Next, we focus on the global features of the flow field in the monolayer as
the active turbulent state is established.

First, we measure both the root-mean-square velocity,
$\sqrt{v_{x}^{2}+v_{y}^{2}}$, and the ensemble average of the sum of both
velocity components, $v_{x}+v_{y}$. In a confluent epithelium, the total force
is approximately zero and the system does not develop any total velocity under
periodic boundary conditions. This is demonstrated in Fig. 6a (dashed line)
for both extensile and contractile systems. In contrast, the rms-velocity has
a finite value and is enhanced with increasing activity for both extensile and
contractile systems. More importantly, the generated flow is interleaved with
vortical structures. This is quantified by calculating the vorticity-vorticity
correlation function that clearly shows the existence of a characteristic
vorticity length scale for different activities. Consistent with continuum
theories Thampi  _et al._ (2014), the vorticity length scale decreases with
increasing activity (Fig. 6b).

In addition to having an intrinsic vorticity length scale that depends on the
activity of the system, active turbulence is shown to manifest exponential
distribution of vortex sizes, that was first predicted using active
nematohydrodynamic equations Giomi (2015), and was confirmed in epithelial
monolayers of Human Bronchial Cells Blanch-Mercader  _et al._ (2018), and in
dense suspensions of microtubule-kinesin motor mixtures Martínez-Prat  _et
al._ (2019). We quantify the vortex areas using Okubo-Weiss criterion Blanch-
Mercader  _et al._ (2018), $\mathcal{W}$:

$\mathcal{W}=-\textup{det}\left[\vec{\nabla}\vec{v}\right].$ (12)

The region is considered to be a vortex if $\mathcal{W}<0$. Fig. 6c
demonstrates the semi-log representation of the total vertex area probability
distribution for different values of activity with $\gamma=0.1$. While we see
a possible evidence of activity-dependent, exponential scaling at intermediate
vortex areas, we could not find a conclusive evidence of the exponential
distribution of vortex areas in our simulations.

Additionally, continuum models of epithelial layers often assume the velocity
field of the cells to be incompressible Saw  _et al._ (2017); Duclos  _et al._
(2018); Pérez-González  _et al._ (2019); Alert  _et al._ (2019). We directly
check this constraint in our cell-based model by calculating the average of
the absolute values of the flow divergence, $\vec{\nabla}\cdot\vec{v}$, for
different values of $\zeta$ and $\gamma$, as shown in Fig. 6d. Interestingly,
as the activity increases, the divergence becomes non-zero, suggesting that
the monolayer of cells deviates from an incompressible behavior. Higher values
of divergence are observed for larger values of elasticity, $\gamma$. The non-
zero divergence further indicates the emergence of dilating-shrinking domains
in the system, which we explore next.

###  III.5 Spontaneous gap formation

While physical properties of confluent monolayers have been intensely
investigated, recent modelling and in vivo experiments have begun to unravel
the importance of extracellular spaces on structural properties of the cell
layers Kim  _et al._ (2021). Interestingly, we find that the emergence of
extracellular spaces can be governed by the activity level of the cells: at a
given cell density, for sufficiently high activities the confluency condition
is broken and extracellular spaces spontaneously emerge (Fig. 7a).

Figure 7: Higher activity values lead to heterogeneity in cell densities,
eventually leading to the formation of highly dense regions and empty spaces,
i.e. gaps. (a) Plots of $\phi_{tot}(t)=\sum_{i}\phi_{i}(t)$ for $\zeta=0.01$
with $\gamma=0.1$ (left) and $\zeta=0.08$ with $\gamma=0.01$ (right) at the
end of the simulation. Black lines represent cell contours, grey – interior of
the cell, and white – exterior. (b) Fraction of empty space on the grid as a
function of time for contractile (left) and extensile (right) activities. The
empty space was evaluated at the end of simulation (15000 time steps) and
averaged over 5 simulations. (c) Cell area distribution (normalized) for
different values of activity, $\zeta$, and elasticity, $\gamma$. Dashed
vertical line corresponds to the area of hexagon, $3\sqrt{3}R^{2}/2$, and the
dotted line to the area of circle, $\pi R^{2}$. The $x$-axis is rescaled by
$\pi R^{2}$.

To quantitatively characterize the spontaneous gap formation, we measure the
fraction of empty space within the cell layer. To this end, in the phase-field
model, we use the threshold method. We calculate the total cell density by
summing over all cells:

$\phi_{tot}=\sum_{i}^{N}\phi_{i}.$ (13)

We then assume that the grid points with $\phi_{tot}<\phi_{0}$ correspond to
empty space. The value of $\phi_{0}=0.2$ is chosen by qualitatively
identifying the gaps.

Interestingly, measurement of the empty space fraction over simulation time
shows that after their formation the fraction reaches a constant value
throughout the simulations. Fig. 7b demonstrates the steady-state fraction of
empty space for different values of activity and elasticity, showing an
elasticity-dependent threshold for activity, above which extracellular spaces
are formed. Beyond this point, the empty space fraction increases as the
activity is increased. Similar gap formation behavior is observed for both
extensile and contractile systems (Fig. 7b), indicating that the spontaneous
emergence of extracellular spaces is a generic feature of an active
multicellular layer.

Furthermore, the crossover from the confluent to non-confluent state is
accompanied by a significant alteration of the individual cell features. This
is best quantified by measuring the distribution of cell areas within the
confluent and non-confluent states of the monolayer (Fig. 7c). When the
monolayer is confluent (left panel in Fig. 7c), the cells obtain more
hexagonal shape, i.e. on average closer to the dashed line, indicating the
area of a hexagon. However, when the confluency is broken (right panel in Fig.
7c), the cells have more circular shape, i.e. on average closer to the dotted
line, indicating the area of a circle, and the cell distributions are narrower
compared to the confluent case. Furthermore, as the activity increases, the
distribution broadens significantly, indicating the emergent activity-induced
heterogeneity in cell areas, and is true for both confluent case and
monolayers with gaps.

##  IV Conclusions

We have presented a phase-field model of cellular monolayers that bridges
cell-level interaction forces to the active nematohydrodynamic behavior at the
multicellular level. We do so by introducing a minimal form of local dipolar
interaction between cells and show that this simple formulation reproduces
phenomenology of the continuum active nematohydrodynamic approach. In
particular, by performing numerical simulations, we show the emergence of bend
and splay instabilities for extensile and contractile interactions. At later
times, the system develops active turbulence and the formation of topological
defects. We quantify the observed features of flows and defects for both
extensile and contractile systems, and show that they are in agreement with
the experimental observations.

Furthermore, we find that when activity and elasticity of individual cells is
sufficiently high, the monolayer flow field deviates from its incompressible
behavior. Concomitantly, we demonstrate that this deviation of the cell layer
flow field from incompressibility leads to spontaneous formation of gaps
within the confluent layers. We observe larger fraction of empty space for
higher activity and lower elasticity. Furthermore, we quantify the emergent
heterogeneity in the cell areas, showing that cells in confluent monolayer
have more hexagonal shapes and higher degree of heterogeneity compared to
cells in monolayer with gaps. The spontaneous gap formation is accompanied by
more homogeneous circular cell shapes.

Inter-cellular gap formation has been observed experimentally in epithelial
cell monolayers, where gaps can reach the length of several cell bodies Zheng
_et al._ (2017). Appearance of cell-free space has been linked to viscoelastic
properties of the substrate Xi  _et al._ (2019). More recently, Sonam et al.
Sonam  _et al._ (2022) have suggested that the spontaneous formation of holes
within cell monolayer can occur due to the heterogeneity in substrate
stiffness and using a vertex model they show that the gaps can form in the
areas close to negative half-integer topological defects, which result in
higher tensile stresses within the cell layer. Thus, extending the model to
account for substrate heterogeneity will provide a theoretical platform to
study the gap formation and its implications on collective cell movement.

## Conflicts of interest

There are no conflicts to declare.

## Acknowledgements

A. D. acknowledges funding from the Novo Nordisk Foundation (grant No.
NNF18SA0035142), Villum Fonden (Grant no. 29476). A. A. and A. D. acknowledge
support from the European Union’s Horizon 2020 research and innovation program
under the Marie Sklodowska-Curie grant agreement No. 847523 (INTERACTIONS).

## References

  * Heisenberg and Bellaïche (2013) C.-P. Heisenberg and Y. Bellaïche, Cell 153, 948 (2013).
  * Brugués  _et al._ (2014) A. Brugués, E. Anon, V. Conte, J. H. Veldhuis, M. Gupta, J. Colombelli, J. J. Muñoz, G. W. Brodland, B. Ladoux, and X. Trepat, Nature Physics 10, 683 (2014).
  * Friedl  _et al._ (2012) P. Friedl, J. Locker, E. Sahai, and J. E. Segall, Nature Cell Biology 14, 777 (2012).
  * Ladoux and Mège (2017) B. Ladoux and R.-M. Mège, Nature Reviews Molecular Cell Biology 18, 743 (2017).
  * Doostmohammadi and Ladoux (2021) A. Doostmohammadi and B. Ladoux, Trends in Cell Biology (2021).
  * Duclos  _et al._ (2018) G. Duclos, C. Blanch-Mercader, V. Yashunsky, G. Salbreux, J.-F. Joanny, J. Prost, and P. Silberzan, Nature Physics 14, 728 (2018).
  * Kawaguchi  _et al._ (2017) K. Kawaguchi, R. Kageyama, and M. Sano, Nature 545, 327 (2017).
  * Prost  _et al._ (2015) J. Prost, F. Jülicher, and J.-F. Joanny, Nature Physics 11, 111 (2015).
  * Doostmohammadi  _et al._ (2018) A. Doostmohammadi, J. Ignés-Mullol, J. M. Yeomans, and F. Sagués, Nature Communications 9, 1 (2018).
  * Maroudas-Sacks  _et al._ (2021) Y. Maroudas-Sacks, L. Garion, L. Shani-Zerbib, A. Livshits, E. Braun, and K. Keren, Nature Physics 17, 251 (2021).
  * Doostmohammadi  _et al._ (2015) A. Doostmohammadi, S. P. Thampi, T. B. Saw, C. T. Lim, B. Ladoux, and J. M. Yeomans, Soft Matter 11, 7328 (2015).
  * Rossen  _et al._ (2014) N. S. Rossen, J. M. Tarp, J. Mathiesen, M. H. Jensen, and L. B. Oddershede, Nature Communications 5, 1 (2014).
  * Doostmohammadi  _et al._ (2016) A. Doostmohammadi, M. F. Adamer, S. P. Thampi, and J. M. Yeomans, Nature Communications 7, 1 (2016).
  * Blanch-Mercader  _et al._ (2018) C. Blanch-Mercader, V. Yashunsky, S. Garcia, G. Duclos, L. Giomi, and P. Silberzan, Physical Review Letters 120, 208101 (2018).
  * Lin  _et al._ (2021) S.-Z. Lin, W.-Y. Zhang, D. Bi, B. Li, and X.-Q. Feng, Communications Physics 4, 1 (2021).
  * Saw  _et al._ (2017) T. B. Saw, A. Doostmohammadi, V. Nier, L. Kocgozlu, S. Thampi, Y. Toyama, P. Marcq, C. T. Lim, J. M. Yeomans, and B. Ladoux, Nature 544, 212 (2017).
  * Marchetti  _et al._ (2013) 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).
  * Peshkov  _et al._ (2014) A. Peshkov, E. Bertin, F. Ginelli, and H. Chaté, The European Physical Journal Special Topics 223, 1315 (2014).
  * Lee and Wolgemuth (2011a) P. Lee and C. W. Wolgemuth, PLoS Computational Biology 7, e1002007 (2011a).
  * Lee and Wolgemuth (2011b) P. Lee and C. Wolgemuth, Physical Review E 83, 061920 (2011b).
  * Gov (2009) N. Gov, HFSP Journal 3, 223 (2009).
  * Camley and Rappel (2017) B. A. Camley and W.-J. Rappel, Journal of Physics D: Applied Physics 50, 113002 (2017).
  * Bi  _et al._ (2016) D. Bi, X. Yang, M. C. Marchetti, and M. L. Manning, Physical Review X 6, 021011 (2016).
  * Henkes  _et al._ (2011) S. Henkes, Y. Fily, and M. C. Marchetti, Physical Review E 84, 040301 (2011).
  * Etournay  _et al._ (2015) 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).
  * Balasubramaniam  _et al._ (2021) L. Balasubramaniam, A. Doostmohammadi, T. B. Saw, G. H. N. S. Narayana, R. Mueller, T. Dang, M. Thomas, S. Gupta, S. Sonam, A. S. Yap, _et al._ , Nature Materials , 1 (2021).
  * Barton  _et al._ (2017) D. L. Barton, S. Henkes, C. J. Weijer, and R. Sknepnek, PLoS Computational Biology 13, e1005569 (2017).
  * Aranson (2016) I. S. Aranson,  _Physical Models of Cell Motility_ (Springer, 2016).
  * Najem and Grant (2016) S. Najem and M. Grant, Physical Review E 93, 052405 (2016).
  * Basan  _et al._ (2013) M. Basan, J. Elgeti, E. Hannezo, W.-J. Rappel, and H. Levine, Proceedings of the National Academy of Sciences 110, 2452 (2013).
  * Sepúlveda  _et al._ (2013) N. Sepúlveda, L. Petitjean, O. Cochet, E. Grasland-Mongrain, P. Silberzan, and V. Hakim, PLoS Computational Biology 9, e1002944 (2013).
  * Monfared  _et al._ (2021) S. Monfared, G. Ravichandran, J. E. Andrade, and A. Doostmohammadi, arXiv preprint arXiv:2108.07657 (2021).
  * Mueller  _et al._ (2019) R. Mueller, J. M. Yeomans, and A. Doostmohammadi, Physical Review Letters 122, 048004 (2019).
  * Zhang  _et al._ (2020) G. Zhang, R. Mueller, A. Doostmohammadi, and J. M. Yeomans, Journal of the Royal Society Interface 17, 20200312 (2020).
  * Zhang and Yeomans (2021) G. Zhang and J. M. Yeomans, arXiv preprint arXiv:2111.14401 (2021).
  * Palmieri  _et al._ (2015) B. Palmieri, Y. Bresler, D. Wirtz, and M. Grant, Scientific Reports 5, 1 (2015).
  * Helfrich (1973) W. Helfrich, Z Naturforsch C 28, 693 (1973).
  * Shao  _et al._ (2010) D. Shao, W.-J. Rappel, and H. Levine, Physical Review Letters 105, 108104 (2010).
  * Biben and Misbah (2003) T. Biben and C. Misbah, Physical Review E 67, 031908 (2003).
  * Peyret  _et al._ (2019) G. Peyret, R. Mueller, J. d’Alessandro, S. Begnaud, P. Marcq, R.-M. Mège, J. M. Yeomans, A. Doostmohammadi, and B. Ladoux, Biophysical Journal 117, 464 (2019).
  * Simha and Ramaswamy (2002) R. A. Simha and S. Ramaswamy, Physical Review Letters 89, 058101 (2002).
  * Ramaswamy and Rao (2007) S. Ramaswamy and M. Rao, New Journal of Physics 9, 423 (2007).
  * Thampi  _et al._ (2014) S. P. Thampi, R. Golestanian, and J. M. Yeomans, EPL (Europhysics Letters) 105, 18001 (2014).
  * Giomi  _et al._ (2011) L. Giomi, L. Mahadevan, B. Chakraborty, and M. Hagan, Physical Review Letters 106, 218101 (2011).
  * Giomi  _et al._ (2013) L. Giomi, M. J. Bowick, X. Ma, and M. C. Marchetti, Physical Review Letters 110, 228101 (2013).
  * Meacock  _et al._ (2021) O. J. Meacock, A. Doostmohammadi, K. R. Foster, J. M. Yeomans, and W. M. Durham, Nature Physics 17, 205 (2021).
  * Sanchez  _et al._ (2012) T. Sanchez, D. T. N. Chen, S. J. DeCamp, M. Heymann, and Z. Dogic, Nature 491, 431 (2012).
  * Kumar  _et al._ (2018) N. Kumar, R. Zhang, J. J. De Pablo, and M. L. Gardel, Science Advances 4, eaat7779 (2018).
  * Wenzel and Voigt (2021) D. Wenzel and A. Voigt, arXiv preprint arXiv:2106.10552 (2021).
  * Giomi (2015) L. Giomi, Physical Review X 5, 031003 (2015).
  * Martínez-Prat  _et al._ (2019) B. Martínez-Prat, J. Ignés-Mullol, J. Casademunt, and F. Sagués, Nature Physics 15, 362 (2019).
  * Pérez-González  _et al._ (2019) C. Pérez-González, R. Alert, C. Blanch-Mercader, M. Gómez-González, T. Kolodziej, E. Bazellieres, J. Casademunt, and X. Trepat, Nature Physics 15, 79 (2019).
  * Alert  _et al._ (2019) R. Alert, C. Blanch-Mercader, and J. Casademunt, Physical Review Letters 122, 088104 (2019).
  * Kim  _et al._ (2021) S. Kim, M. Pochitaloff, G. A. Stooke-Vaughan, and O. Campàs, Nature Physics , 1 (2021).
  * Zheng  _et al._ (2017) J. Y. Zheng, S. P. Han, Y.-J. Chiu, A. K. Yip, N. Boichat, S. W. Zhu, J. Zhong, and P. Matsudaira, Biophysical Journal 113, 1585 (2017).
  * Xi  _et al._ (2019) W. Xi, T. B. Saw, D. Delacour, C. T. Lim, and B. Ladoux, Nature Reviews Materials 4, 23 (2019).
  * Sonam  _et al._ (2022) S. Sonam, L. Balasubramaniam, S.-Z. Lin, Y. M. Y. Ivan, I. P. Jaumà, C. Jebane, M. Karnat, Y. Toyama, P. Marcq, J. Prost, _et al._ , bioRxiv (2022).