
2
tissue, the energy can be written as
E=1
2X
cκA(Ac−Ac0)2+κP(Pc−Pc0)2,(1)
where each cell experiences an energy cost for deviation
of its area Acand perimeter Pcfrom target values Ac0
and Pc0, with area and perimeter stiffness κAand κP.
The first term on the RHS models 3D cell volume incom-
pressibility via an effective 2D area elasticity [19,37].
The second describes the competition between cell corti-
cal contractility and adhesion between neighbouring cells
in controlling cell-edge tension and perimeter [7,19,37].
We denote by ⃗
Fn=−δE
δ⃗xnthe total force on the nth
vertex of the tiling at position ⃗xndue to interactions with
all other vertices. In an applied shear of rate ˙γ, with flow
direction xand shear gradient y, we assume over-damped
dynamics with drag ζ,d⃗xn
dt =ζ−1⃗
Fn+ ˙γyn⃗
ˆx, with Lees-
Edwards periodic boundary conditions. The cells also
undergo T1 topological neighbor exchanges that allow
the tissue to plastically relax stresses [9,38–40].
To focus on amorphous tissue structures, we simulate a
50 : 50 bidisperse tiling of Nc= 4096 cells of target areas
A0= 1,1.4, which sets our length unit. We adjust Pc0for
the two cell populations to maintain the target cell shape
p0=Pc0/√Ac0the same for all cells. We choose units
in which κA= 1 and ζ= 1 and set κp= 1.0 throughout.
We vary p0and the imposed shear rate ˙γ. As an initial
condition, we seed a planar Voronoi tiling then evolve
the above dynamics to steady state in zero shear. At
time t= 0, we switch on shear and measure the shear
stress Σij (t) = 1
NPN
n=1 Fnixnj , where the sum is over
all Nvertices in the tiling, and the mean cell perimeter
p(t) = 1
NcPNc
c=1 pc. Denoting by ⃗
tncthe unit vector along
the edge of length lncbetween the ncth and (nc+ 1)th
vertices of cell c, we define a single-cell shape tensor σc
ij =
1
νcPνc
n=1 lnctnc
itnc
j, where the sum is over the νcvertices
of the c-th cell, and the tissue-scale averaged orientation
tensor σij =1
NcPNc
c=1 σc
ij . We use the same notation
Σij , σij , p for the counterpart coarse-grained quantities
in our constitutive model below.
In the absence of external stress, the vertex model ex-
hibits a liquid solid transition as a function of the target
shape p0[7,41]. For p0< p∗
0the energy barriers to T1
transitions are finite and the system is a solid with a
finite zero-frequency linear shear modulus. At the criti-
cal value p∗
0, the mean energy barrier for T1 transitions
vanishes, giving liquid response for p0> p∗
0. For our
bidisperse tiling, p∗
0= 3.85. For monodisperse disordered
polygons p∗
0≃3.81, a value close to that of a regular
pentagon [7]. This value is renormalized by motility [9]
and by cell alignment with local spontaneous shear [42].
It was recently realized that this transition has a geo-
metric origin associated with the underconstrained na-
ture of the energy in Eq. 1[20,22,43]. For regular
hexagons the transition occurs at the isoperimetric value
piso =p8√3≃3.722. Below this value it is not possible
to satisfy both target area and perimeter and the ground
state has p=p∗
0and finite energy. This is the solid or
incompatible state. For p0> piso there is a family of
zero energy area and perimeter preserving ground states,
with p=p0. The system can accommodate an exter-
nally applied linear shear by adjusting its shape within
this degenerate manifold [22]. The compatible system is
therefore a liquid with zero shear modulus, although it
stiffens and acquires rigidity at finite strains [27].
Constitutive model — We now construct a continuum
model that accounts for the mean-field liquid-solid transi-
tion, and also captures the key rheological features of the
vertex model: (i) reversibility of linear response to small
strains, (ii) strain stiffening at intermediate strains, (iii)
plastic relaxation at larger strains, due to T1 cell rear-
rangements, and (iv) a yield stress in the steady state
flow curve Σ( ˙γ), as obtained in Ref. [27]. Although our
model below is cast in frame invariant form, capable of
addressing any flow, we focus on response to simple shear,
to compare with our vertex model simulations.
We assume dynamics of the cell perimeter governed by:
˙p+vk∇kp= ˙γ−1
τp
(p−p0)(p−p∗
0−ασij σij ),(2)
with αand τpconstants and invariant strain rate ˙γ=
p2Dij Dij . In the absence of shear, prelaxes on a
timescale τpto a steady state that displays a transcriti-
cal bifurcation as a function of the target cell perimeter
p0, with p=p∗
0in the solid phase p0< p∗
0and p=p0
in the liquid phase p0> p∗
0, capturing the liquid-solid
transition [7]. The same transcritical structure emerges
by writing exact equations for the relaxation of a single
cell modeled as a regular n−sided polygon according to
the vertex model dynamics prescribed above.
In shear, the perimeter is advected by flow and
stretched by the shear rate ˙γ. In addition, the coupling
ασij σij captures a key intuition of our approach: that a
shear-induced global cell orientation σij provides an ef-
fective mean field that distorts the individual cell’s shape
paway from its zero-shear value. As a result, in the solid
phase pincreases relative to its zero shear value p=p∗
0
from the outset of straining. In the liquid phase, pin-
creases relative to its zero shear value p=p0only after a
critical strain amplitude γc, capturing the strain-induced
stiffening transition [27]. The behavior introduced by the
coupling of single-cell shape, as quantified by the mean
perimeter p, to the tissue-scale cell shape σij is analogous
to the influence of cell alignment due to internally gener-
ated stresses in Drosophila germband extension [42]. In-
deed, the form of coupling of pto σij in Eqn. 2is justified
both by experiment [42] and mean field theory [22,27].
The cell orientation tensor is taken to obey an evolu-
tion equation of the widely used Maxwellian form,
˙σij +vk∇kσij =σikKkj +Kkiσkj + 2Dij −aσij ,(3)
where Kij =∂jviis the strain rate tensor and Dij =
1
2(Kij +Kji). The last term in Eq. 3describes plastic
relaxation. It vanishes in linear response (small strains),