#  Alignment interactions drive structural transitions in biological tissues

Matteo Paoluzzi<EMAIL_ADDRESS>Departament de Física de la Matèria
Condensada, Universitat de Barcelona, C. Martí Franquès 1, 08028 Barcelona,
Spain.  Luca Angelani  ISC-CNR, Institute for Complex Systems, Piazzale A.
Moro 2, I-00185 Rome, Italy.  Dipartimento di Fisica, Sapienza Università di
Roma Piazzale A. Moro 2, I-00185 Rome, Italy.  Giorgio Gosti  Center for Life
Nano Science, Istituto Italiano di Tecnologia, Viale Regina Elena 291,
I-00161, Rome, Italy.  M Cristina Marchetti  Department of Physics, University
of California Santa Barbara, Santa Barbara, CA 93106, USA  Ignacio
Pagonabarraga  CECAM Centre Européen de Calcul Atomique et Moléculaire, École
Polytechnique Fédérale de Lausanne (EPFL), Batochime, Avenue Forel 2, 1015
Lausanne, Switzerland.  UBICS University of Barcelona Institute of Complex
Systems, Martí i Franquès 1, E08028 Barcelona, Spain.  Departament de Física
de la Matèria Condensada, Universitat de Barcelona, C. Martí Franquès 1, 08028
Barcelona, Spain.  Giancarlo Ruocco  Center for Life Nano Science, Istituto
Italiano di Tecnologia, Viale Regina Elena 291, I-00161, Rome, Italy
Dipartimento di Fisica, Sapienza Università di Roma Piazzale A. Moro 2,
I-00185 Rome, Italy.

###### Abstract

Experimental evidence shows that there is a feedback between cell shape and
cell motion. How this feedback impacts the collective behavior of dense cell
monolayers remains an open question. We investigate the effect of a feedback
that tends to align the cell crawling direction with cell elongation in a
biological tissue model. We find that the alignment interaction promotes
nematic patterns in the fluid phase that eventually undergo a non-equilibrium
phase transition into a quasi-hexagonal solid. Meanwhile, highly asymmetric
cells do not undergo the liquid-to-solid transition for any value of the
alignment coupling. In this regime, the dynamics of cell centers and shape
fluctuation show features typical of glassy systems.

††preprint: APS/123-QED

##  I Introduction

Eukaryotic cells at high packing fraction organize themselves into confluent
monolayers, develop collective motion, and trigger a variety of patterns that
play a fundamental role in complex biological processes ranging from wound
healing to metastasis invasion [1, 2]. Pattern formation in biological tissues
involves length scales that are much bigger than the typical cell length. This
observation suggests that a coarse-grained model of biological tissues needs
to take into account only a few key ingredients of the single-cell dynamics
[3]. Different approaches have been developed during the last few decades to
capture the large-scale behavior of biological tissues [3, 4, 5, 6].

Experimental studies discovered that biological tissues show glassy dynamics,
support viscoelastic response, and behave as a disordered soft material in the
vicinity of jamming or glassy transition [7, 8, 9, 10, 11, 12, 13, 14, 15,
16]. However, differently from particulate systems, cell shape anisotropy is
the driver of the jamming transition in confluent monolayers [17, 18, 19].
Thus, cell shape and its fluctuations are important ingredients that have to
be taken into account in a mesoscopic description.Shape fluctuations can be
introduced in different ways [4, 5, 20]. Among the other alternatives, Vertex
and Voronoi models are successful coarse-grained descriptions that have been
tested against different experiments in the last few years [17, 21, 19, 22,
23, 24, 25, 26, 27, 28, 29, 30].

Figure 1:  Pictorial representation of the model. The blue line represents the
direction $\hat{\bf{n}}_{i}$ of the largest eigenvalue of the shape tensor
(the direction in the lab frame is parametrized by the angle $\alpha_{i}$).
The red arrow is the self-propulsion direction $\hat{\bf{e}}_{i}$ (that is
parametrized by the angle $\phi_{i}$). The alignment interaction tends to
reduce the distance $\alpha_{i}-\phi_{i}$.  Figure 2:  Representative
snapshots of steady-state configurations taken in the weakly nematic liquid
phase (a), disordered solid phase (b), and in the glassy regime (c). The color
indicates the modulus of the velocity in the lab frame from zero (dark) to its
maximum value (yellow). The alignment interaction acts as an inverse effective
temperature: As the strength of the interaction increases, velocity
fluctuations become strongly inhibited promoting solidification. In the glassy
regime, we observe the formation of dynamically heterogeneous regions
(highlighted by the dashed white circles).

Because cells can move autonomously, biological tissue can be seen as a soft
and active material [31]. It has been shown that feedback mechanisms at the
single-cell level can trigger collective motion [22, 32, 33, 34]. Besides
promoting collective migration, alignment interactions can also change the
structural properties of a biological tissue [22]. Structural changes and
morphological transitions play a fundamental role in morphogenesis and
organogenesis [35], however, the key ingredients responsible for self-
organization and spatial differentiation in organoids are still poorly
understood. Consequently, isolating the few fundamental ingredients that play
the role of control parameters for the emergent structural organization allows
us to gain insight into complex biological processes.

Cells become elongated during motion and tend to move along the direction of
their long axis. As a consequence, cell motility is correlated with cell
anisotropy [36, 37]. Recent studies on phase-field models show that minimal
dipolar interactions in monolayers of isotropic cells promote spontaneous
symmetry breaking and nematic order [38].

In this paper, we introduce a generalization of the Voronoi model of
biological tissues where we consider a minimal alignment interaction between
cell shape and cell displacement. We show that the feedback between shape and
displacement triggers morphological transitions in the confluent monolayer.
The alignment interaction acts as an inverse effective temperature that cools
down the system as the intensity of the interaction increases. Starting from
fluid configurations, the liquid becomes weakly nematic as the interaction is
turned on. For higher values of the alignment interaction, the system falls
into a hexagonal disordered solid [39, 40]. We observe that the alignment
interaction promotes the formation of cooperative clusters that tend to slow
down the dynamics and trigger the proliferation of dynamical heterogeneities
typical of glassy systems. Glassy dynamics involve both, the correlation of
density fluctuations, as in the case of supercooled liquid, and shape
fluctuations.

##  II Kinetic Monte Carlo Voronoi Model

We implement a kinetic Monte Carlo (KMC) dynamics based on a Voronoi model of
biological tissues [41]. The confluent monolayer is represented through the
Voronoi tessellation of the $N$ cell centers (labeled by $i\\!=\\!1,...,N$, of
coordinates $\mathbf{r}_{i}\\!=\\!(x_{i},y_{i})$, in a two dimensional square
box of side $L\\!=\\!\sqrt{N}$ with periodic boundary conditions). Let
$\mathbf{r}\\!\equiv\\!(\mathbf{r}_{1},...,\mathbf{r}_{N})$ be a configuration
of the system. The dynamics is governed by the following configurational
energy [41, 42, 27, 18, 19, 43]

$E[{\mathbf{r}}]=\sum_{i}\left[K_{a}(a_{i}-a_{0})^{2}+K_{p}(p_{i}-p_{0})^{2}\right]\;,$
(1)

where the function $a_{i}$ and $p_{i}$ return the value of the area and the
perimeter of the $i-$th polygon of the Voronoi tessellation. Cell area and
cell perimeter fluctuate around the preferred (or target) values $a_{0}$ and
$p_{0}$, their fluctuations are regulated by the stiffnesses $K_{a}$ and
$K_{p}$. The square deviation from $a_{0}$ enforces the constraint of
incompressibility in three dimensions. The square deviation from $p_{0}$
encodes the competition between cell-cell adhesion and active contractility in
the actomyosin cortex [19]. In the following we set $a_{0}\\!=\\!1$,
$K_{p}\\!=\\!K_{a}\\!=\\!1$, and we express (1) in terms of the target shape
index $s_{0}\\!=\\!p_{0}/\sqrt{a_{0}}$ [43]. We sample stationary
configurations of the energy functional (1) numerically using a Monte Carlo
algorithm in which we propose time-correlated trial moves for the cell
centers. This algorithm is general enough to capture two key features of the
real cell movement dynamics: (i) Cells usually move at a velocity that can
fluctuate in magnitude [44], this is ensured by the noise parameter $T$ that
enters in the Monte Carlo algorithm, and (ii) Cells displace positively
correlated steps on a microscopic time scale $\tau$ [45, 44], which is the
second parameter of the algorithm. Voronoi and Vertex models develop a rough
energy landscape where energy barriers separate local minima [18, 19, 43]. MC
algorithms are particularly suitable for reaching steady-state configurations
in these situations [46, 47].

Here the persistent motion typical of active systems is modeled using time-
correlated trial displacements. We also consider an alignment interaction
acting between the direction where the cell is elongated and the crawling
direction. Indicating with $\bm{\delta}_{i,t}$ the displacement performed by
the cell $i$ at the time step $t$, and adopting polar coordinates, we can
write $\bm{\delta}_{i,t}=\delta_{i,t}(\cos\phi_{i,t},\cos\phi_{i,t})$, with
$\delta_{i,t}=|\bm{\delta}_{i,t}|$. The angle $\phi_{i,t}$ determines the
direction of the displacement. The shape tensor
$\mathbf{Q}_{i}\\!=\\!(\sum_{l=1}^{N_{v}^{i}}\Delta\mathbf{r}_{l,i}\otimes\Delta\mathbf{r}_{l,i})/N_{v}^{i}$
encodes information on cell shape, where $N_{v}^{i}$ is the number of cell $i$
vertices, $\Delta\mathbf{r}_{l,i}=\mathbf{r}_{l,i}-\mathbf{r}_{CM,i}$ is cell
$i$’s $l-$th vertex position, $\mathbf{r}_{CM,i}$ indicates the center of
mass, and the symbol $\\!\otimes\\!$ indicates the standard diadic product. In
our case, $\mathbf{Q}_{i}$ is a $2\times 2$ symmetric matrix. The eigenvector
corresponding to the largest eigenvalue defines the direction of maximum cell
elongation. We explore the effect of a nematic alignment interaction between
the principal axis of elongated cell shape (parametrized through the angle
$\alpha_{i}$) and the crawling direction (parametrized by the angle
$\phi_{i}$). The alignment interaction is sketched in Fig. (1). For enforcing
the alignment interaction, at each MC step, we update the new direction
$\phi_{i}$ using the following rule

$\phi_{i,t+1}=\phi_{i,t}-J\sin 2(\phi_{i,t}-\alpha_{i,t})\;.$ (2)

At the beginning of the time step $t\\!+\\!1$, the alignment interaction in
Eq. (2) tends to align the (trial) displacement performed during the previous
time step $\bm{\delta}_{i,t}$. Once updated the displacement direction, we
propose the trial move

$\displaystyle\mathbf{r}_{i,t+1}\\!=\\!\mathbf{r}_{i,t}\\!+\\!\bm{\delta}_{i,t}\;,$
(3)

and then we evolve the displacement

$\displaystyle\bm{\delta}_{i,t}\\!=\\!\bm{\delta}_{i,t-1}\\!+\\!\delta_{1}\bm{\eta}_{i}$
(4)

with the condition $\bm{\delta}_{i,0}\\!=\\!\delta_{0}\bm{\eta}_{i}$ [48], i.
e., in this way, the trial moves are correlated on a time scale
$\tau\\!=\\!(\delta_{0}/\delta_{1})^{2}\tau_{MC}$ [49, 48] (the Monte Carlo
time step $\tau_{MC}$ is defined as the succession of $N$ elementary moves
[47]). The components of the random vector $\bm{\eta}_{i}$ are extracted from
a uniform distribution independently at each time step. The distribution is
centered around zero and has unit variance. Moreover, following Refs. [50,
48], the displacements are constrained to be
$|\bm{\delta}_{i,t}|\\!\leq\\!\delta_{0}$ and $\delta_{0}\\!\geq\\!\delta_{1}$
(in our simulations $\delta_{0}=0.25$).

The time evolution of the displacement $\bm{\delta}_{i,t}$ introduces a
correlation on the time scale $\tau$ so that
$\langle\delta_{i,t}\delta_{j,s}\rangle\sim\delta_{i,j}e^{-|t-s|/\tau}$ [49,
48], as well as in the case of Active Ornstein-Uhlenbeck particles [51, 52,
53, 54, 55]. The model interpolates between an equilibrium relaxation dynamics
for $\tau=0$, representing cells that perform a random crawling, to a
persistent non-equilibrium dynamics characterized by a ballistic regime on
short time scales, which is the hallmark of self-propelled motion at low
Reynolds numbers. It is important to stress that, although the time evolution
of $\phi_{i}$ given by Eq. (2) is deterministic, once we evolve
$\bm{\delta}_{i,t}$ with Eq. (4), the stochastic term
$\delta_{1}\bm{\eta}_{i}$ introduces a rotational noise on $\phi_{i}$ that
makes it to diffuse with a rotational diffusion constant
$D_{r}\propto\tau^{-1}$. The trial move is thus accepted with probability
$P_{acc}\\!\propto\\!\exp{(-\Delta E/T)}$, with $\Delta E\equiv
E[\mathbf{r}_{t+1}]-E[\mathbf{r}_{t}]$.

Figure 3:  Structural properties. (a) phase diagram of the model using
$D_{eff}$ as a dynamical order parameter. The dashed black line indicates the
solid/liquid transition. Crossing the liquid-solid transition the mean-squared
displacement changes discontinuously, as shown in panel (b) ($J\in[0,4]$ from
violet to green, and $s_{0}=3.0$). Approaching the glassy region, the mean-
squared displacement shows a subdiffusive regime on intermediate times, as
shown in (c) for $s_{0}=3.8$ (same values of $J$ shown in (b)). The dashed red
line is the diffusive scaling $\Delta r^{2}\sim t$, dashed blue line the
ballistic scaling $\Delta r^{2}\sim t^{2}$. (d) $D_{eff}$ as a function of the
nematic coupling $J\in[0,4]$ for different values of $s_{0}\in[3.0,3.9]$,
increasing values from violet to yellow. (e) Phase diagram using $\psi_{6}$ as
a structural order parameter (increasing values of $\psi_{6}$ from blue to
green). The dashed black line represents the transition line to the hexagonal
solid. The dotted white line corresponds to the value $s_{0}=s_{0}^{hex}\sim
3.722$. (f) Order parameter $\psi_{6}$ as a function of $J$ for different
values of $s_{0}\in[3.0,3.9]$, increasing values from violet to yellow. Static
structure factor $S(q_{x},q_{y})$ for $s_{0}=3.0$ (g) and $s_{0}=3.6$ (h). (i)
Nematic order parameter $S$ as a function of $J$ (increasing values of
$s_{0}\in[3.0,3.9]$ from violet to yellow). (j-l) Radial distribution function
$g_{\parallel,\perp}(r)$ evaluated along $g_{\parallel}(r)=g(r_{\parallel},0)$
and perpendicular $g_{\perp}(r)=g(0,r_{\perp})$ the nematic director
$\mathbf{n}$ ((j) liquid, (k) solid, and (l) glassy).

We perform simulations of a tissue composed of $N\\!=\\!100,400,1600,6400$
cells with $T\\!=\\!0.002,0.05$ and $\tau\\!=\\!200$ (in Monte Carlo unit
$\tau_{MC}$).

##  III Phase diagram

We study the phase diagram of the tissue using as control parameters $s_{0}$
and $J$. The target shape index $s_{0}$ tunes the typical cell asphericity,
i.e., the larger is $s_{0}$ the more elongated is the cell [19]. We anticipate
that the system shows a liquid-solid transition that is driven by the
alignment coupling $J$. We perform numerical simulations in a region of the
phase diagram where the system at finite temperature $T$ behaves as a fluid
for any values of $\tau$ at $J\\!=\\!0$ (the phase diagram in the $T$ vs
$\tau$ plane for $J\\!=\\!0$ is shown in Appendix (B)).

In Fig. (2) we report three representative snapshots taken in the liquid,
solid, and glassy regime. We observe an increase in velocity fluctuations as
$s_{0}$ increases (see the Appendix (C)). On the other hand, $J$ triggers the
formation of islands of slow-moving cells, as shown in Fig. (2). For
$J\\!\neq\\!0$, the system develops nematic patterns as signaled by a non-
vanishing value of the nematic order parameter $S$. We monitored the average
cell anisotropy via the parameter
$\Delta=\langle\frac{(\lambda_{1}^{i}-\lambda_{2}^{i})^{2}}{(\lambda_{1}^{i}+\lambda_{2}^{i})^{2}}\rangle$,
where the eigenvalues of the shape tensor $\lambda_{1,2}^{i}$ are sorted in a
way so that $\lambda_{1}^{i}>\lambda_{2}^{i}$. Our analysis reveals a jump to
lower values of $\Delta$ crossing the liquid-to-solid transition. We
anticipate that, in the solid phase, the system arranges into hexagonal
patches with $\Delta\neq 0$ but small. Moreover, from the study of the
relaxation dynamics of the principal axis, we obtain that axis fluctuations
decorrelate on a finite time scale and no flipping dynamics between the two
principal directions occurs (see the Appendix (F)). These findings show that
the alignment interaction is always well defined.

##  IV Structural Properties

We start our quantitative discussion from the phase diagram which is shown in
Fig. (3)-(a). The phase diagram has been obtained considering the long time
behavior of the mean-squared displacement $\Delta r^{2}$ (defined in the
Appendix (A)) as a dynamic order parameter [19], i. e., looking at the
effective diffusion constant $D_{eff}\\!\equiv\\!\lim_{t\to\infty}\Delta
r^{2}/(4t)$ [19, 22]. The typical behavior of $\Delta r^{2}$ is shown in (b)
and in (c) for $s_{0}\\!=\\!3.0,3.8$. As one can see, $\Delta r^{2}$ undergoes
a crossover from a ballistic regime on short time scales, i.e., $\Delta
r^{2}\\!\sim\\!t^{2}$, to a diffusive regime, i.e., $\Delta
r^{2}\\!\sim\\!4D_{eff}t$, for longer times. As $J$ increases in intensity, we
observe two different behaviors in $\Delta r^{2}$. For small $s_{0}$ values
(b), $\Delta r^{2}$ discontinuously develops a plateau right after the
ballistic regime. This fact signals a liquid-solid transition. For larger
$s_{0}$ values (c), the plateau is replaced by a subdiffusive regime. Panel
(d) shows $D_{eff}$ as a function of $J$. Through this analysis, we identify
three regimes in the phase diagram (see panel (a)): a liquid phase for small
$J$ values, a solid-state, at larger $J$ values and small $s_{0}$, and a
glassy regime, where $\Delta r^{2}$ develops a subdiffusive behavior at
intermediate times.

To gain insight into the structural properties of the system, we take a look
at the positional and orientational order. We start our discussion with the
order parameter $\psi_{6}$ (see the Appendix (A) for its definition) for
revealing the presence of sixfold order. $\psi_{6}$ leads to the phase diagram
that is shown in Fig. (3-e), where the color map indicates the magnitude of
the order parameter. The behavior of $\psi_{6}$ indicates that the solid phase
is characterized by hexagonal order (see Fig. (3-f))). In agreement with early
studies on Vertex models [27, 29], for higher $J$ values, the transition
between a glassy fluid and a disordered hexatic solid phase matches the
critical value of a regular hexagon $s_{0}\\!=\\!s_{0}^{hex}\\!\sim\\!3.722$.
Complementary information about the positional order is provided by the static
structure factor $S(q_{x},q_{y})$ which allows us allows to visualize the
emerging ordered patches. In Fig. (3)-(g,h) we report $S(q_{x},q_{y})$ in the
glassy (g) and in the solid (h) phase. The solid phase shows hexatic patches
that are compatible with a disordered hexagonal solid. The increase in
positional order in the solid regime is signaled by marked damped oscillations
in the $g(r)$ that imply the lack of a true crystalline structure ($g(r)$ is
reported in the Appendix (D)). The order parameter $\psi_{6}$ jumps almost
discontinuously at the transition (see panel (f)) providing evidence for an
increase of hexatic order in the solid phase rather than in the liquid [39,
40].

The control parameter that triggers the transition between liquid and solid is
the alignment coupling $J$ which plays the role of an (inverse) effective
temperature. For rationalizing this effect, we consider the simplest case
where each cell is represented by a self-propelled spheroid undergoing an
active Brownian dynamics with self-propulsion velocity $v_{0}$, and rotational
diffusion $\tau^{-1}$. During the dynamics, each particle tends to (i)
minimize the mechanical energy, and (ii) align towards the direction given by
$\alpha_{i}$. We indicate with
$\mathbf{r}^{0}=({\mathbf{r}}_{1},\ldots,{\mathbf{r}}_{N})$ the inherent state
configuration that minimizes $E[\mathbf{r}]$, i.e., $\left.\nabla
E\right|_{{\mathbf{r}}={\mathbf{r}}^{0}}=0$, and we linearize the dynamics
around those minima [56, 21]. The equations of motion for the fluctuations are
$\dot{\delta\mathbf{r}_{i}}\\!=\\!v_{0}\mathbf{e}_{i}-\mu\mathbf{M}_{ij}\delta\mathbf{r}_{j}$,
with $\delta\mathbf{r}_{i}\\!=\\!\mathbf{r}_{i}-\mathbf{r}_{i}^{0}$, $\mu$ the
mobility, and $\mathbf{M}_{ij}$ a $2\\!\times\\!2$ block of the dynamical
matrix 111We adopt Einstein summation convention. The orientation
$\mathbf{e}_{i}\\!=\\!(\cos\theta_{i},\sin\theta_{i})$ follows the linearized
equation $\dot{\theta}_{i}\\!=\\!-J(\theta_{i}-\alpha_{i})\\!+\\!\eta_{i}$,
with $\langle\eta_{i}\rangle\\!=\\!0$ and
$\langle\eta_{i}(t)\eta_{j}(s)\rangle\\!=\\!2\tau^{-1}\delta_{ij}\delta(t-s)$.
Performing the replacement $\theta_{i}\to\theta_{i}\\!-\\!\alpha_{i}$ and
projecting the equations for $\delta\mathbf{r}_{i}$ onto the normal modes (see
Refs. [20, 21] and the Appendix (E) for details), we obtain that the mean
energy per mode can be written as
$e_{\nu}\\!=\\!\frac{1}{2}T_{eff}^{0}I(\tau,J)$ with
$T_{eff}^{0}\\!=\\!v_{0}^{2}\tau/2$. For $J\\!=\\!0$ and $\tau\\!\neq\\!0$, it
follows a generalization of the equipartition theorem [58]. Using the
expression of the energy per mode, we can thus define
$T_{eff}(\tau,J)\\!=\\!T_{eff}^{0}I(\tau,J)$. In general, it is not possible
to compute analytically $T_{eff}(\tau,J)$, however, it turns out that it is a
decreasing function of $J$ bounded above by $T_{eff}^{0}$ and below by
$T_{eff}^{0}/\lambda_{\nu}\tau$.

We conclude our analysis of the structural properties of the tissue by
studying the features of the liquid state. As shown in Fig. (3)-(i), where we
report the nematic order parameter $S$ (defined in Appendix (A)), in the
liquid phase the system develops weak nematic order. We can thus define two
preferred global directions that are individuated by the average direction of
the director fields $\mathbf{n}\\!=\\!\frac{1}{N}\sum_{i}(\cos
2\alpha_{i},\sin 2\alpha_{i})$ computed at a given time step. We indicate with
$x_{\parallel}$ and $x_{\perp}$ respectively the directions parallel and
orthogonal to $\mathbf{n}$. In the nematic phase, the positional order of the
cell centers is different along these two directions, as it is shown in Fig.
(3)-(j) where we report the radial distribution functions $g_{\parallel}(r)$
and $g_{\perp}(r)$. The space isotropy is restored in the glassy regime (see
panel (l)). The sixfold orientational order replaces the twofold orientational
order in the solid phase (see (k) and (h)).

Figure 4:  Dynamical slowing down. (a) The self-part of the intermediate
scattering function $F(q,t)$ ($J\in[2.3,3.2]$, increasing values from violet
to yellow) and corresponding four-point susceptibility $\chi_{4}(q,t)$ (c).
The shape parameters is $s_{0}\\!=\\!3.73$ and $T=0.002$. (b) The correlation
function $C_{\psi}(t)$. Inset in (c): Map of displacements for $J\\!=\\!3.1$
and $t\\!\sim\\!10^{2}$. (d) Structural relaxation time as a function of the
nematic coupling $J$, the dashed black curve is a fit to Vogel-Fulcher-Tammann
law with $T\to J^{-1}$.

##  V Relaxation dynamics

We now probe the region of the phase diagram where the tissue develops a
subdiffusive regime. Fig. (4)-(a) shows the behavior of the intermediate
scattering function $F_{s}(q,t)$ for $s_{0}\\!=\\!3.73$, $T=0.002$, and
$q\\!=\\!q_{peak}$ (with $q_{peak}$ the position of the first peak of the
static structure factor $S(q)$). We have also measured the time correlation
function $C_{\psi}(t)$ of the hexatic order parameter $\psi_{6}$ [59, 60] (see
the Appendix (A) for the definition). The behavior of $C_{\psi}$ is shown in
panel (b) of the same figure. As one can see, $C_{\psi}$ undergoes a dynamical
slowing down as $J$ increases similar to that observed in $F_{s}(q,t)$. Since
$\psi_{6}(t)$ reflects the local structure that is determined by the number of
cell sides at time $t$, a non-vanishing correlation $C_{\psi}(t)$ signals a
dynamical slowing down of shape fluctuations. The dynamical slowing down is
usually due to the presence of relaxation dynamics on different time scales.
The emerging of complex and heterogeneous relaxation dynamics becomes more
evident probing the dynamical susceptibility $\chi_{4}(q,t)$ defined as the
sample-to-sample fluctuations of $F_{s}(q,t)$, shown in (c) [50].
$\chi_{4}(q,t)$ shows a broad peak, due to the presence of dynamical
heterogeneity (see the displacement field, inset in (c)), that grows in height
and shifts towards longer times as $J$ increases, the typical feature of
glassy systems approaching the glass transition [50, 61]. We can provide a
quantitative measure of the dynamical slowing down using a characteristic
relaxation time $\tau_{\alpha}$ defined as $C(\tau_{\alpha})\\!=\\!e^{-1}$,
with $C(t)$ a time-correlation function. In panel (d) we show the behavior of
the relaxation time for $F_{s}(q_{peak},t)$, and $C_{\psi}(t)$. Another
estimate of $\tau_{\alpha}$ is provided by the position of the peak of
$\chi_{4}$ (shown in the same panel). As it has been observed in the case of
flocking transition in biological tissues [22], the behavior of
$\tau_{\alpha}$ as a function of $J$ proves that the alignment interaction
acts as an inverse temperature, causing a cooling down of the system as $J$
increases. This is confirmed by the fact that a Vogel-Fulcher-Tamman formula
$\tau_{\alpha}\propto\exp(B/(J^{-1}-J^{-1}_{c}))$ (with $T\to J^{-1}$ and
$J_{c}\sim 2.9$) well captures the behavior of the relaxation time.

##  VI Discussion and Conclusions

The collective behavior of biological tissues shows features remarkably
similar to those of active nematics, disordered solids, and supercooled
liquids. Some of these facts can be rationalized in the framework of dense
active matter [56]. However, to capture the collective properties of
biological tissues in an opportune coarse-grained description one needs to
take into account all the relevant ingredients of single-cell dynamics. Since
moving cells assume an asymmetric configuration that spontaneously breaks
spatial symmetries, we have studied numerically how feedback between cell
shape and displacement changes the structural properties of the tissue. We
performed our study within the framework of the Voronoi models. We focused our
attention on alignment interactions tending to couple the direction of cell
motion with its elongation. Besides the experimental evidence at the single-
cell level that highlights the importance of feedback in cell locomotion [62],
the impact of these interactions on the large-scale behavior of confluent
monolayers remains poorly understood.

We have explored the phase diagram of the tissue using as a control parameter
the target shape index $s_{0}$, which is experimentally accessible [17, 63].
The second control parameter is the strength of the alignment interaction
$J>0$. We documented that the interplay of these control parameters triggers
structural changes giving rise to a rich phase diagram characterized by
liquid-to-solid transitions and glassy dynamics. For large enough $J$ values,
the tissue undergoes a phase transition between a disordered state and a quasi
hexagonal lattice at $s_{0}\sim s_{0}^{hex}$ [27].

For larger $s_{0}$, the system remains in a disordered liquid state, showing
typical features of glassy dynamics as the strength of the alignment force
increases [7, 19, 64, 65]. In particular, we observed the proliferation of
dynamical heterogeneities, subdiffusive dynamics, broad peaks in the dynamical
susceptibilities, and dynamical slowing down of density and shape
fluctuations. We showed that the solidification of the system for increasing
values of $J$ can be generally understood in terms of an effective temperature
$T_{eff}$ that scales with the inverse of $J$, similarly to what has been
observed in the case of the self-propelled Voronoi model with polar
interactions [22].

In conclusion, the intensity of the shape-displacement feedback at the single-
cell level can trigger structural transitions in confluent monolayers that
impact dramatically the collective behavior of the biological tissue. Our
results suggest that, if the coupling between cell elongation and displacement
is small, and thus $J$ assumes small values, cells tend to rearrange as in a
weakly nematic fluid rather than form a confluent monolayer. In other words,
biochemical mechanisms that tend to alter cell polarization might impact the
collective properties of the biological tissue. There are biological
implications for a weak coupling between shape and locomotion. For instance,
the interplay between these two properties might become weak for cells that
suddenly lose a polarized shape, as in the case of metastasis invasion, where
the metastatic cell does not show epithelial polarity [66]. Our analysis
suggests that, in such a condition, the tissue tends to melt into a fluid
phase. More in general and in agreement with recent studies that revealed the
crucial role of cell symmetry [67], our results show that alignment
interactions might provide an additional control parameter for the epithelial-
mesenchymal transition.

### Acknowledgments

M.P. has received funding from the European Union’s Horizon 2020 research and
innovation programme under the MSCA grant agreement No 801370 and by the
Secretary of Universities and Research of the Government of Catalonia through
Beatriu de Pinós program Grant No. BP 00088 (2018). I.P. acknowledges MICINN,
DURSI and SNSF for financial support under Projects No. PGC2018-098373-B-I00,
No. 2017SGR-884, and No. 200021-175719, respectively. M.P. acknowledges
Sapienza University of Rome, Physics Department, for kind hospitality during
the preparation of this manuscript. MCM was supported by the US National
Science Foundation though award No. DMR-2041459.

##  Appendix A Observables

We indicate with ${\mathbf{r}}_{i}(t)$ the position of the $i-$th cell center
in the lab reference frame and with ${\mathbf{r}}_{i}^{\prime}(t)$ the
position in the center-of-mass reference frame. In the following, we indicate
with $\langle\mathcal{O}\rangle_{s}$ the average of the observable
$\mathcal{O}$ with respect to independent runs, i.e., the subscript $s$
indicates sample averages. We indicate with $\langle\mathcal{O}\rangle_{t}$
time-averaging in the stationary state.

For studying the single-cell diffusion and the solid-to-liquid transition, we
compute the mean-squared displacement $\Delta r^{2}$ that is

$\Delta
r^{2}=\frac{1}{N}\left\langle\sum_{i}\left[{\mathbf{r}}^{\prime}_{i}(t)-{\mathbf{r}}_{i}^{\prime}(0)\right]^{2}\right\rangle_{s}\;.$
(5)

We quantify the emerging of hexatic order through the complex field
$\Psi_{i}(t)$ defined for each cell

$\Psi_{i}(t)=\frac{1}{n}\sum_{j\in n.n.}^{n}e^{i6\theta_{ij}(t)}$ (6)

with $n$ the number of Voronoi neighbors to the cell $i$. The angle
$\theta_{ij}$ is individuated by the two cell centers $i$ and $j$. The hexatic
order parameter at the time step $t$ reads

$\psi_{6}(t)=\frac{1}{N}\left|\sum_{i}\Psi_{i}(t)\right|$ (7)

and we indicate with $\psi_{6}$ its time average. We obtain additional and
complementary information on the positional order by measuring the static
structure factor $\tilde{S}(q_{x},q_{y})$ that is

$\tilde{S}(q_{x},q_{y})=\frac{1}{N}\left\langle\sum_{j,k}e^{i\mathbf{q}\cdot({\mathbf{r}}_{j}^{\prime}-{\mathbf{r}}_{k}^{\prime})}\right\rangle_{t,s}$
(8)

where the wave vector $\mathbf{q}=(q_{x},q_{y})$ satisfies the periodic
conditions imposed to the dynamics, i.e.,
$q_{x,y}=\frac{2\pi}{L}(n_{x},n_{y})$, with $n_{x,y}=0,\pm 1,\pm 2,...$ (and
avoiding the combination $n_{x}=n_{y}=0$).

We detect the presence of nematic order measuring the order parameter

$S=2\langle\frac{1}{N}\sum_{i}\cos^{2}\alpha_{i}\rangle-1\;.$ (9)

We obtain additional information on the nematic phase measuring the pair
distribution function

$g(\mathbf{r})=\frac{1}{N}\left\langle\sum_{i,j\neq
i}\delta(\mathbf{r}-\mathbf{r}_{j}+\mathbf{r}_{i})\right\rangle_{t,s}\;.$ (10)

In particular, we compute $g_{\parallel,\perp}(r)$, where the subscription
indicates that the observable is computed along the principal directions of
the nematic director, i.e., $g_{\perp}(r)\equiv g(0,r_{\perp})$ and
$g_{\parallel}(r)\equiv g(r_{\parallel},0)$.

As dynamical observables, we measure the self-part of the Intermediate
Scattering Function $F_{s}(q,t)$, and the time-correlation function of the
hexatic order parameter $C_{\psi}(t)$ [68, 50, 59, 60]. The intermediate
scattering function is

$F_{s}(q,t)=\frac{1}{N}\left\langle\sum_{i}e^{-i\mathbf{q}\cdot({\mathbf{r}}^{\prime}_{i}(t)-{\mathbf{r}}^{\prime}_{i}(0))}\right\rangle_{s}\;,$
(11)

where $\mathbf{q}$ follows the same prescription used to compute
$S(q_{x},q_{y})$. The sample-to-sample fluctuations of $F_{s}(q,t)$ provides a
measure of dynamical heterogeneity through the four-point dynamical
susceptibility $\chi_{4}(q,t)$. The position of the peak $\chi_{4}(q,t)$,
i.e., $t=\tau_{4}$, individuates the typical time scale of dynamical
heterogeneity. We thus compute the displacement field
$\Delta{\mathbf{r}}(x,y,\tau_{4})$[50]. Furthermore, we measure the relaxation
time of shape fluctuations using $C_{\psi}(t)$ defined through the correlation
function

$C_{\psi}(t)=\frac{1}{C_{\psi}(0)}\left\langle\sum_{i}\Psi_{i}(t)\Psi_{i}^{*}(0)\right\rangle_{s}$
(12)

##  Appendix B Phase diagram for $J=0$

Figure 5:  Phase diagram of the model for $J=0$ and $\tau=20,200,2000$, (a),
(b), and (c) panel, respectively. The dashed black line indicates the
transition between solid and fluid using the shape parameter $q$ as a
structural order parameter. The dashed green line in panel (b) is the region
explored in the main text for $J\neq 0$.

In the main text, we have used as control parameters $s_{0}$ and $J$. For
$J=0$, the model reduces to a self-propelled Voronoi model where the self-
propulsion is due to correlated noise. We have thus probed different regions
of the phase diagram for $J=0$ using as control parameters the strength of the
noise $T$ and the shape index $s_{0}$. The resulting phase diagrams for
$\tau=20,200,2000$ and $N=256$ are shown in Fig. (5). As a structural
parameter for discriminating the solid from the fluid phase, we adopt the
shape parameter $q$ defined as

$q=\left\langle\sum_{i}\frac{p_{i}}{\sqrt{a_{i}}}\;.\right\rangle$ (13)

Following [21, 19, 28], we use the criterium $q>3.81$ for defining the fluid
regime and $q<3.81$ for the solid regime. We explored noise values
$T\in[10^{-4},10^{-1}]$ and $s_{0}\in[3.0,3.7]$. As one can see, for
$\tau=200$ the system is always in the fluid state. The region of parameters
explored in the main text is highlighted in panel (b). The phase diagram shows
the same qualitative features as for the self-propelled Voronoi model [21].

##  Appendix C Collective motion

The presence of migratory patterns has been evaluated using global quantities
as the Vicsek order parameter $\Phi$ [69] and the nematic order parameter
$\Omega$ of the displacements. The order parameter $\Phi$ captures collective
cell migration and it is defined as follows

$\Phi=\left\langle\frac{1}{N}\left|\sum_{k}e^{i\phi_{k}(t)}\right|\right\rangle_{t}$
(14)

The order parameter $\Omega$ captures the emerging of nematic order in the
cell displacements. It is defined as follows

$\Omega=\left\langle\frac{1}{N}\left|\sum_{k}e^{i2\phi_{k}(t)}\right|\right\rangle_{t}$
(15)

The resulting phase diagrams are shown in Fig. (6).

Figure 6:  Collective motion. Global nematic order (a) and polar order (b) in
the velocity field. The dashed black line indicates the transition to solid.

The system develops a weak nematic phase in the velocity field in the liquid
phase for $0<J<3$, as shown in panel (a). It is worth noting that the order
parameter does not overcome the value of $\Omega\sim 0.3$, indicating that
only in a small system fraction nematic order is appreciable. In panel (b) we
report the behavior of $\Phi$. Around the liquid-solid transition, the system
develops weak flocking patterns for $J\sim 3$. In this case, the Vicsek order
parameter does not overcome $\Phi\sim 0.4$.

Figure 7:  Order parameters $\Omega$ and $\Phi$ as a function of $J$ for
increasing values of $s_{0}$ (from violet to yellow).

The behavior of $\Omega$ and $\Phi$ as a function of $J$ for different values
of $s_{0}$ is shown in Fig. (7). The parameter $\Omega$ turns out to be
different from zero in the liquid state, almost independently by $s_{0}$
(panel a). The presence of migratory patterns characterized by polar order
($\Phi\neq 0$), is more evident for small values of $s_{0}$ (panel b).

Figure 8:  Representative snapshots for $p_{0}=3.2$ for different values of
$J$. In the first row (a-g), the color indicates the modulus of the velocity
of each cell. In the second (h-n) and in the third (o-u) row the color
indicates nematic and polar angle, respectively, obtained from the velocity of
the cell and calculated with respect to the $x$-axis. The fourth row (1-2)
shows the distribution of the velocity for $J=0.0$ (panel (1)) and $J=3.6$
(panel (2)).

The snapshots of steady-state configurations for $s_{0}=3.2$ are shown in Fig.
(8). In the first row, cells are colored according to their velocity. In the
second row, the color indicates the angle of the nematic and polar order,
respectively. These parameters have been obtained considering the velocity
$\mathbf{v}_{i}$ of the cell $i$ that can be written as
$\mathbf{v}_{i}=v_{i}(\cos\theta_{i},\sin\theta_{i})$. Region of the same
color indicates local nematic/polar order. In the fourth row, we report the
probability distribution function of the velocity $\mathcal{P}(v)$.
$\mathcal{P}(v)$ turns out to be strongly peaked around zero in the solid
phase and develops a long tile towards higher values in the liquid state.

##  Appendix D Static Structure Factor $S(q_{x},q_{y})$ and Radial
Distribution Function $g(r)$

In Fig. (9), we show the static structure factor $S(q_{x},q_{y})$. Panels (a),
(b), and (c) report results for $s_{0}=3.0$ as the intensity of the alignment
interaction grows. Crossing the critical value $J\sim 3$, the heat map
develops patterns peculiar to the hexatic phase. For a larger value of
$s_{0}=3.8$, the system does not undergo a liquid-solid transition anymore. In
this situation, the structure factor does not develop regular peaks (panels
(d), (e), and (f)).

Figure 9:  Static Structure factor $S(q_{x},q_{y})$ for $s_{0}=3.0$ (a,b,c)
and $s_{0}=3.8$ (d,e,f). Figure 10:  Radial distribution function $g(r)$ for
$s_{0}=3.0$ (a) and $s_{0}=3.7$ (b). Increasing values of
$J=0.0,0.5,1.0,1.5,2.0,2.5,3.0,3.5,5.0$ from violet to yellow, respectively.
For clarity curves have been shifted vertically. The grey area in (a)
indicates configurations in the solid phase.

In Fig. (10) we report the radial distribution function $g(r)$ for $s_{0}=3.0$
and $s_{0}=3.7$. In both cases, $g(r)$ becomes more structured as the nematic
coupling increases. However, while at high $s_{0}$ values the liquid becomes
progressively more structured, at low $s_{0}$ values when the system crosses
from a liquid-like to a solid-like phase, $g(r)$ changes suddenly. The
presence of damped oscillations indicates a disordered rather than a
crystalline structure.

##  Appendix E Effective temperature with alignment interactions

We consider a system composed of $N$ self-propelled particles where each of
them tends to align towards a given direction defined by the angle
$\alpha_{i}$. The self-propulsion of magnitude $v_{0}$ acts along
$\mathbf{e}_{i}=(\cos\theta_{i},\sin\theta_{i})$ and changes direction with a
rate $\tau^{-1}$. We indicate with $\mathbf{r}_{i}^{0}$ the inherent state
configuration that minimizes the mechanical energy of the system. We indicate
with $\delta\mathbf{r}_{i}=\mathbf{r}_{i}-\mathbf{r}_{i}^{0}$ a small
displacement around the equilibrium configuration. Linearizing the potential
around the minimum of the mechanical energy we obtain the following equations
of motion

$\displaystyle\delta\dot{\mathbf{r}}_{i}$ $\displaystyle=$ $\displaystyle
v_{0}\mathbf{e}_{i}+\mu\mathbf{M}_{ij}\delta\mathbf{r}_{i}\;,$ (16)
$\displaystyle\dot{\theta}_{i}$ $\displaystyle=$
$\displaystyle-J(\theta_{i}+\alpha_{i})+\eta_{i}\;,$ (17)

where the random force satisfies $\langle\eta_{i}\rangle=0$ and
$\langle\eta_{i}(t)\eta_{j}(s)\rangle=2\tau^{-1}\delta_{ij}\delta(t-s)$.
$\mathbf{M}_{ij}$ is the $2\times 2$ block of the dynamical matrix (hereafter
we adopt the Einstein summation convention). Expanding the perturbation in
terms of the normal modes $\mathbf{u}_{i}^{\lambda}$ of the dynamical matrix,
one has

$\delta\mathbf{r}_{i}=a_{\nu}(t)\,\mathbf{u}_{i}^{\nu}$ (18)

and the amplitude $a_{\nu}$ follows the equation of motion

$\dot{a}_{\nu}=-\mu\lambda_{\nu}a_{\nu}+\tilde{\eta}_{\nu}\;.$ (19)

The noise $\tilde{\eta}_{\nu}$ satisfies

$\displaystyle\langle\tilde{\eta}_{\nu}(t)\rangle$ $\displaystyle=$
$\displaystyle M(t)$ (20)
$\displaystyle\langle\tilde{\eta}_{\nu}(t)\tilde{\eta}_{\mu}(s)\rangle$
$\displaystyle=$ $\displaystyle\delta_{\mu,\nu}C(t-s)$ (21) $\displaystyle
C(t-s)$ $\displaystyle\equiv$
$\displaystyle\frac{v_{0}^{2}}{2}\langle\cos(\theta(t)-\theta(s))\rangle$ (22)

In terms of the probability density function $\rho(\theta,t)$, the solution of
the Fokker-Planck equation generated by (17) with the initial condition
$\rho(\theta,t=0)=\delta(\theta-\theta_{0})$ reads

$\displaystyle\rho(\theta,t)$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2\pi\sigma^{2}(t)}}e^{-\frac{(\theta-M(t))^{2}}{2\sigma(t)^{2}}}$
(23) $\displaystyle\sigma^{2}(t)$ $\displaystyle\equiv$
$\displaystyle\frac{1}{\tau J}(1-e^{-2Jt})$ (24) $\displaystyle M(t)$
$\displaystyle\equiv$ $\displaystyle\theta_{0}e^{-Jt}$ (25)

and thus we can write

$\langle\cos\Delta\theta\rangle=\cos M(t)e^{-\frac{1}{2}\sigma^{2}(t)}$ (26)

where we have defined $\Delta\theta\equiv\theta(t)-\theta(0)$ and, without
loss of generality we have set $s=0$ and $\theta_{0}=0$. Through (19) we can
compute the average potential energy stored in each mode

$\displaystyle e_{\nu}=\langle\frac{1}{2}\lambda_{\nu}a_{\nu}^{2}\rangle$
$\displaystyle=$ $\displaystyle\frac{1}{2}T_{eff}I(\tau,J)$ (27)
$\displaystyle T_{eff}^{0}$ $\displaystyle\equiv$
$\displaystyle\frac{v_{0}^{2}\tau}{2\mu}$ (28) $\displaystyle I(\tau,J)$
$\displaystyle\equiv$
$\displaystyle\int_{0}^{\infty}\frac{dt}{\tau}e^{-\frac{1}{2}\sigma^{2}(t)-\mu\lambda_{\nu}t}\;.$
(29)

Where we have introduced the effective temperaure $T_{eff}^{0}$ that is one of
the control parameters of the KMC algorithm. In this way, the equilibrium
equipartition theorem is recovered in the limit $\tau=0$. For $\tau\neq 0$ and
$J=0$, we recover a generalization of the equipartition theorem, as shown in
Ref. [58]

$e_{\nu}(J=0)=\frac{1}{2}\frac{T_{eff}^{0}}{1+\mu\lambda_{\nu}\tau}\;.$ (30)

Another limiting case is obtained for $J\to\infty$ (and equivalently
$\tau\to\infty$), for which

$e_{\nu}(J\to\infty)=\frac{1}{2}\frac{T_{eff}^{0}}{\lambda_{\nu}\tau}\;.$ (31)

As a consequence, the effective temperature $T_{eff}(\tau,J)$ turns out to be
bounded above by $T_{eff}^{0}$ and decreases towards
$T_{eff}^{0}/\lambda_{\nu}\tau$ as $J$ increases.

##  Appendix F Cell anisotropy

Indicating with $\lambda_{1,2}^{i}$ the eigenvalues of the shape tensor of the
cell $i$ and using the convention $\lambda_{1}^{i}>\lambda_{2}^{i}$, we define
the cell asphericity $\Delta_{i}$ as

$\displaystyle\Delta_{i}=\frac{(\lambda_{1}^{i}-\lambda_{2}^{i})^{2}}{(\lambda_{1}^{i}+\lambda_{2}^{i})^{2}}\;.$
(32)

Figure 11:  Asphericity $\Delta$ for $s_{0}=3.0$ as a function of $J$.

The asphericity provides a quantitative measure of cell anisotropy, in the
case of highly symmetric cells one has $\lambda_{1}^{i}\simeq\lambda_{2}^{i}$
and thus $\Delta_{i}\simeq 0$. Contrarily, for rod-like cells, $\Delta_{i}\to
1$, i.e., cells strongly elongated towards a given direction. As shown in Fig.
(11), the parameter $\Delta=\langle\Delta_{i}\rangle$, where the angular
parentheses indicate both averages, over cells and steady-state
configurations, jumps from higher to lower values crossing the liquid-to-solid
transition. We notice that even in the solid phase, where the system arranges
in hexagonal patterns, $\Delta\neq 0$, indicating that cells are not displaced
in a perfectly regular hexagonal lattice.

Figure 12:  Relaxation time of the direction of maximum elongation for
$s_{0}=1$ as a function of $J$.

Moreover, shape fluctuations do not flip the principal axes, making the
alignment interaction always well defined. For proving this, we measure the
relaxation time $\tau_{\alpha_{i}}$ of the eigenvector corresponding to
$\lambda_{1}^{i}$. As shown in Fig. (12) the relaxation time is of the order
of hundreds of $\tau_{MC}$ in the liquid phase, and it becomes larger in the
solid state, indicating that, although the asphericity is small, the direction
corresponding to the larger eigenvalue decays on longer times.

Figure 13:  Asphericity $\Delta$ for $J=1$ as a function of $s_{0}$.

Finally, in Fig. (13) we report the behavior of $\Delta$ for $J=1$ as a
function of $s_{0}$. As one can see, although the system does not undergo to
structural changes, cells become more elongated making the system more fluid.

## References

  * Sunyer _et al._ [2016] R. Sunyer, V. Conte, J. Escribano, A. Elosegui-Artola, A. Labernadie, L. Valon, D. Navajas, J. M. García-Aznar, J. J. Muñoz, P. Roca-Cusachs _et al._ , _Science_ , 2016, 353, 1157–1161.
  * Friedl and Mayor [2017] P. Friedl and R. Mayor, _Cold Spring Harbor perspectives in biology_ , 2017, 9 (4), a029199.
  * Trepat and Sahai [2018] X. Trepat and E. Sahai, _Nature Physics_ , 2018, 14, 671–682.
  * Alert and Trepat [2020] R. Alert and X. Trepat, _Annual Review of Condensed Matter Physics_ , 2020, 11, 77–101.
  * Camley and Rappel [2017] B. A. Camley and W.-J. Rappel, _Journal of physics D: Applied physics_ , 2017, 50, 113002.
  * Banerjee and Marchetti [2019] S. Banerjee and M. C. Marchetti, in _Cell Migrations: Causes and Functions_ , Springer, 2019, pp. 45–66.
  * Angelini _et al._ [2011] T. E. Angelini, E. Hannezo, X. Trepat, M. Marquez, J. J. Fredberg and D. A. Weitz, _Proceedings of the National Academy of Sciences_ , 2011, 108, 4714–4719.
  * Manning _et al._ [2010] M. L. Manning, R. A. Foty, M. S. Steinberg and E.-M. Schoetz, _Proceedings of the National Academy of Sciences_ , 2010, 107, 12517–12522.
  * Garcia _et al._ [2015] S. Garcia, E. Hannezo, J. Elgeti, J.-F. Joanny, P. Silberzan and N. S. Gov, _Proceedings of the National Academy of Sciences_ , 2015, 112, 15314–15319.
  * Kasza _et al._ [2007] K. E. Kasza, A. C. Rowat, J. Liu, T. E. Angelini, C. P. Brangwynne, G. H. Koenderink and D. A. Weitz, _Current opinion in cell biology_ , 2007, 19, 101–107.
  * Trepat _et al._ [2009] X. Trepat, M. R. Wasserman, T. E. Angelini, E. Millet, D. A. Weitz, J. P. Butler and J. J. Fredberg, _Nature physics_ , 2009, 5, 426–430.
  * Angelini _et al._ [2010] T. E. Angelini, E. Hannezo, X. Trepat, J. J. Fredberg and D. A. Weitz, _Physical review letters_ , 2010, 104, 168104.
  * Pawlizak _et al._ [2015] S. Pawlizak, A. W. Fritsch, S. Grosser, D. Ahrens, T. Thalheim, S. Riedel, T. R. Kießling, L. Oswald, M. Zink, M. L. Manning _et al._ , _New Journal of Physics_ , 2015, 17, 083049.
  * Schötz _et al._ [2008] E.-M. Schötz, R. D. Burdine, F. Jülicher, M. S. Steinberg, C.-P. Heisenberg and R. A. Foty, _HFSP journal_ , 2008, 2, 42–56.
  * Saraswathibhatla and Notbohm [2020] A. Saraswathibhatla and J. Notbohm, _Phys. Rev. X_ , 2020, 10, 011016\.
  * Saraswathibhatla _et al._ [2021] A. Saraswathibhatla, S. Henkes, E. E. Galles, R. Sknepnek and J. Notbohm, _Extreme Mechanics Letters_ , 2021, 48, 101438.
  * Park _et al._ [2015] J.-A. Park, J. H. Kim, D. Bi, J. A. Mitchel, N. T. Qazvini, K. Tantisira, C. Y. Park, M. McGill, S.-H. Kim, B. Gweon, R. S. J. Jacob Notbohm, S. Burger, S. H. Randell, A. T. Kho, D. T. Tambe, C. Hardin, S. A. Shore, E. Israel, D. A. Weitz, D. J. Tschumperlin, E. P. Henske, S. T. Weiss, M. L. Manning, J. P. Butler, J. M. Drazen and J. J. Fredberg, _Nature Materials_ , 2015, 14, 1040–1048.
  * Bi _et al._ [2014] D. Bi, J. Lopez, J. Schwarz and M. L. Manning, _Soft Matter_ , 2014, 10, 18885–1890.
  * Bi _et al._ [2016] D. Bi, X. Yang, M. C. Marchetti and M. L. Manning, _Phys. Rev. X_ , 2016, 6, 021011.
  * Henkes _et al._ [2020] S. Henkes, K. Kostanjevec, J. M. Collinson, R. Sknepnek and E. Bertin, _Nature communications_ , 2020, 11, 1–9.
  * Bi _et al._ [2016] D. Bi, X. Yang, M. C. Marchetti and M. L. Manning, _Phys. Rev. X_ , 2016, 6, 021011.
  * Giavazzi _et al._ [2018] F. Giavazzi, M. Paoluzzi, M. Macchi, D. Bi, G. Scita, M. L. Manning, R. Cerbino and M. C. Marchetti, _Soft Matter_ , 2018, 14, 3471–3477.
  * Wang _et al._ [2020] X. Wang, M. Merkel, L. B. Sutter, G. Erdemci-Tandogan, M. L. Manning and K. E. Kasza, _Proceedings of the National Academy of Sciences_ , 2020.
  * Sharp _et al._ [2019] T. A. Sharp, M. Merkel, M. L. Manning and A. J. Liu, _PloS one_ , 2019, 14, e0209892.
  * Sussman _et al._ [2018] D. M. Sussman, J. M. Schwarz, M. C. Marchetti and M. L. Manning, _Physical review letters_ , 2018, 120, 058001.
  * Merkel and Manning [2018] M. Merkel and M. L. Manning, _New Journal of Physics_ , 2018, 20, 022002\.
  * Staple _et al._ [2010] D. B. Staple, R. Farhadifar, J.-C. Röper, B. Aigouy, S. Eaton and F. Jülicher, _The European Physical Journal E_ , 2010, 33, 117–127.
  * Sussman _et al._ [2018] D. M. Sussman, M. Paoluzzi, M. C. Marchetti and M. L. Manning, _EPL (Europhysics Letters)_ , 2018, 121, 36001.
  * Farhadifar _et al._ [2007] R. Farhadifar, J.-C. Röper, B. Aigouy, S. Eaton and F. Jülicher, _Current Biology_ , 2007, 17, 2095–2104.
  * Li _et al._ [2021] Y.-W. Li, L. L. Y. Wei, M. Paoluzzi and M. P. Ciamarra, _Phys. Rev. E_ , 2021, 103, 022607.
  * Marchetti _et al._ [2013] M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao and R. A. Simha, _Rev. Mod. Phys._ , 2013, 85, 1143–1189.
  * Henkes _et al._ [2011] S. Henkes, Y. Fily and M. C. Marchetti, _Phys. Rev. E_ , 2011, 84, 040301\.
  * Petrolli _et al._ [2019] V. Petrolli, M. Le Goff, M. Tadrous, K. Martens, C. Allier, O. Mandula, L. Hervé, S. Henkes, R. Sknepnek, T. Boudou _et al._ , _Physical review letters_ , 2019, 122, 168101.
  * Czajkowski _et al._ [2018] M. Czajkowski, D. Bi, M. L. Manning and M. C. Marchetti, _Soft matter_ , 2018, 14, 5628–5642.
  * Lancaster and Knoblich [2014] M. A. Lancaster and J. A. Knoblich, _Science_ , 2014, 345, 283\.
  * Maeda _et al._ [2008] Y. T. Maeda, J. Inose, M. Y. Matsuo, S. Iwaya and M. Sano, _PloS one_ , 2008, 3, e3734.
  * Lauffenburger and Horwitz [1996] D. A. Lauffenburger and A. F. Horwitz, _cell_ , 1996, 84, 359–369.
  * Mueller _et al._ [2019] R. Mueller, J. M. Yeomans and A. Doostmohammadi, _Physical review letters_ , 2019, 122, 048004.
  * Li and Ciamarra [2018] Y.-W. Li and M. P. Ciamarra, _Phys. Rev. Materials_ , 2018, 2, 045602\.
  * Pasupalak _et al._ [2020] A. Pasupalak, L. Yan-Wei, R. Ni and M. P. Ciamarra, _Soft Matter_ , 2020, 16, 3914–3920.
  * Honda [1978] H. Honda, _Journal of theoretical biology_ , 1978, 72(3), 523–43.
  * Nagai and Honda [2001] T. Nagai and H. Honda, _Philos. Mag. B_ , 2001, 81, 669.
  * Bi _et al._ [2015] D. Bi, J. H. Lopez, J. M. Schwarz and M. L. Manning, _Nature Phyics_ , 2015, 11, 1074.
  * Vedel _et al._ [2013] S. Vedel, S. Tay, D. M. Johnston, H. Bruus and S. R. Quake, _Proceedings of the National Academy of Sciences_ , 2013, 110, 129–134.
  * Dieterich _et al._ [2008] P. Dieterich, R. Klages, R. Preuss and A. Schwab, _Proceedings of the National Academy of Sciences_ , 2008, 105, 459–463.
  * Jorgensen and Tirado-Rives [1996] W. L. Jorgensen and J. Tirado-Rives, _The Journal of Physical Chemistry_ , 1996, 100, 14508–14513.
  * Binder [1995] K. Binder, _Monte Carlo and molecular dynamics simulations in polymer science_ , Oxford University Press, 1995.
  * Levis and Berthier [2014] D. Levis and L. Berthier, _Phys. Rev. E_ , 2014, 89, 062301.
  * Berthier [2014] L. Berthier, _Phys. Rev. Lett._ , 2014, 112, 220602.
  * Berthier and Biroli [2011] L. Berthier and G. Biroli, _Rev. Mod. Phys._ , 2011, 83, 587–645.
  * Maggi _et al._ [2015] C. Maggi, U. M. B. Marconi, N. Gnan and R. Di Leonardo, _Scientific reports_ , 2015, 5, .
  * Szamel _et al._ [2015] G. Szamel, E. Flenner and L. Berthier, _Physical Review E_ , 2015, 91, 062304.
  * Fily and Marchetti [2012] Y. Fily and M. C. Marchetti, _Phys. Rev. Lett._ , 2012, 108, 235702\.
  * Fodor _et al._ [2016] E. Fodor, C. Nardini, M. E. Cates, J. Tailleur, P. Visco and F. van Wijland, _Phys. Rev. Lett._ , 2016, 117, 038103.
  * Paoluzzi _et al._ [2016] M. Paoluzzi, C. Maggi, U. Marini Bettolo Marconi and N. Gnan, _Phys. Rev. E_ , 2016, 94, 052602.
  * Henkes _et al._ [2018] S. Henkes, M. C. Marchetti and R. Sknepnek, _Phys. Rev. E_ , 2018, 97, 042605.
  * Note [1] We adopt Einstein summation convention.
  * Maggi _et al._ [2014] C. Maggi, M. Paoluzzi, N. Pellicciotta, A. Lepore, L. Angelani and R. Di Leonardo, _Phys. Rev. Lett._ , 2014, 113, 238303.
  * Flenner and Szamel [2015] E. Flenner and G. Szamel, _Nature communications_ , 2015, 6, 1–6.
  * Massana-Cid _et al._ [2018] H. Massana-Cid, J. Codina, I. Pagonabarraga and P. Tierno, _Proceedings of the National Academy of Sciences_ , 2018, 115, 10618–10623.
  * Kob _et al._ [1997] W. Kob, C. Donati, S. J. Plimpton, P. H. Poole and S. C. Glotzer, _Phys. Rev. Lett._ , 1997, 79, 2827–2830.
  * Trepat and Fredberg [2011] X. Trepat and J. J. Fredberg, _Trends in cell biology_ , 2011, 21, 638–646.
  * Malinverno _et al._ [2017] C. Malinverno, S. Corallino, F. Giavazzi, M. Bergert, Q. Li, M. Leoni, A. Disanza, E. Frittoli, A. Oldani, E. Martini, T. Lendenmann, G. Deflorian, G. V. Beznoussenko, D. Poulikakos, K. H. Ong, M. Uroz, X. Trepat, D. Parazzoli, P. Maiuri, W. Yu, A. Ferrari, R. Cerbino and G. Scita, _Nature Materials_ , 2017, 16, 587–596.
  * Nandi _et al._ [2018] S. K. Nandi, R. Mandal, P. J. Bhuyan, C. Dasgupta, M. Rao and N. S. Gov, _Proceedings of the National Academy of Sciences_ , 2018, 115, 7688–7693.
  * Sadhukhan and Nandi [2021] S. Sadhukhan and S. K. Nandi, _Phys. Rev. E_ , 2021, 103, 062403\.
  * Paul _et al._ [2017] C. D. Paul, P. Mistriotis and K. Konstantopoulos, _Nature Reviews Cancer_ , 2017, 17, 131.
  * Jung _et al._ [2019] H.-Y. Jung, L. Fattet, J. H. Tsai, T. Kajimoto, Q. Chang, A. C. Newton and J. Yang, _Nature cell biology_ , 2019, 21, 359–371.
  * Lačević _et al._ [2003] N. Lačević, F. W. Starr, T. B. Schrøder and S. C. Glotzer, _The Journal of Chemical Physics_ , 2003, 119, 7372–7387.
  * Vicsek and Zafeiris [2012] T. Vicsek and A. Zafeiris, _Physics Reports_ , 2012, 517, 71 – 140\.