# A Mechanistic Model of the Organization of Cell Shapes in Epithelial Tissues

Kanaya Malakar 1, Rafael Rubenstein 1, Dapeng Bi 2 and Bulbul Chakraborty 1,∗

###### Abstract

The organization of cells within tissues plays a vital role in various
biological processes, including development and morphogenesis. As a result,
understanding how cells self-organize in tissues has been an active area of
research. In our study, we explore a mechanistic model of cellular
organization that represents cells as force dipoles that interact with each
other via the tissue, which we model as an elastic medium. By conducting
numerical simulations using this model, we are able to observe organizational
features that are consistent with those obtained from vertex model
simulations. This approach provides valuable insights into the underlying
mechanisms that govern cellular organization within tissues, which can help us
better understand the processes involved in development and disease.

##  1 Keywords:

tissue, elasticity, organization, Monte Carlo

##  2 Introduction

Organ surfaces are covered with confluent monolayers of epithelial cells or
endothelial cells, providing physical barriers for organs and bodies. Under
healthy conditions, the cells in these layers remain static and non-migratory,
meaning that they do not move or change position. This feature is an important
aspect of their function as a physical barrier, as any gaps or spaces between
the cells could allow harmful substances to pass through. Instead, the cells
remain securely connected to each other through a network of tight junctions,
adherens junctions, and desmosomes, creating a solid-like structure. The
cellular collective behaves as a solid-like object because of the cohesive
forces that hold the cells together. Together, these forces create a strong
and stable structure that resists deformation and maintains the integrity of
the organ surface.

Inside a confluent tissue, cells form a polygonal tiling of space without any
gaps between them. In most cases, they form an amorphous tiling without any
spatial order. In spite of major differences between their microscopic
constituents and interactions, these cellular assemblies bear a strong
resemblance to jammed granular solids. Cells within dense tissues can be
considered as particles in a jammed granular solid, where the packing and
arrangement of cells play a crucial role in determining the overall tissue
properties. This packing can be influenced by various factors, such as cell
size, shape, and cell-cell interactions, as well as the mechanical
microenvironment.

The essential features that are shared by these two distinct classes of
amorphous solids include (i) both form via out-of-equilibrium, non-thermal
processes, (ii) their elasticity develops in response to external stresses,
and (iii) the cells within tissues are instantaneously in a state of
mechanical equilibrium with each cell maintaining force and torque balance.
These properties define jammed solids: solids whose rigidity emerge in
response to external stresses Cates et al. (1998). A further consequence is
that cells in tissues can be represented by force dipoles since force-balance
eliminates the force-monopole contribution. Recent work has uncovered a
universal statistical distribution that governs cell shapes in cell monolayers
across a wide range of tissues and biological processes Atia et al. (2018).
This distribution has been observed in various contexts, including the
maturation of the pseudostratified bronchial epithelial layer in both
asthmatic and non-asthmatic donors, as well as in the formation of the ventral
furrow in the Drosophila embryo. These findings imply a relationship between
jamming and geometry that extends beyond specific system details, encompassing
both living organisms and non-living jammed systems.

Crucially, tissues differ from jammed granular solids in their display of a
complex pattern of cell attributes, including the shapes of cells. Our primary
hypothesis is that, in a solid-like tissue, the organization of cell shapes is
driven by the constraints of mechanical equilibrium. This mechanistic
perspective naturally leads to a model of interacting force dipoles Schwarz
and Safran (2013). Recently, this same mechanistic perspective has been framed
as an “emergent theory” of elasticity of jammed solids Nampoothiri et al.
(2022, 2020). From this viewpoint, the organization of cells in tissues and
the elasticity of tissues, emerges from the constraints of force and torque
balance, locally, as the collection of cells respond to externally imposed
forces, which are the natural “charges” of this emergent gauge theory
Nampoothiri et al. (2022, 2020). The emergent elasticity theory parallels that
of classical elasticity, with two crucial differences: physical displacements
or strains do not appear, and the elastic moduli are not material properties
but depend on how these non-equilibrium, jammed solids are created.

In this paper, we adopt this perspective, however, in practice, our model
reduces to the force-dipole models that have been extensively used to study
cell-cell interactions in tissues Schwarz and Safran (2013), if we assume a
set of elastic moduli for a tissue. In future work, we would like to determine
these elastic moduli, phenomenologically, through the measurement of stress-
stress correlations Vinutha et al. (2023); Nampoothiri et al. (2022).

The article is divided into the following sections. In section 2, we define
the mechanistic model, and discuss how this is related to other models of cell
organizations in tissues and the Vertex model(Farhadifar et al. (2007); Bi et
al. (2016); Li et al. (2018); Yan and Bi (2019); Das et al. (2021)). In
section 3, we present the results from numerical simulations, and discuss
relationships with the vertex model, focusing on the appearance of
orientational order and the organization of the magnitudes of cell
polarizations (aspect ratios). In section 4, we present a summary and
discussion of future work.

##  3 Theoretical Model

Each cell is modeled as a force dipole, with two equal and opposite forces
acting along the same line Schwarz and Safran (2013). These force dipoles
interact with each other by virtue of the fact that they are embedded in an
elastic medium, the tissue. The force dipole associated with a cell is given
by the tensor Bischofs et al. (2004a):

$P_{\alpha,\beta}=P\hat{n}_{\alpha}\hat{n}_{\beta}\equiv
Pq_{\alpha,\beta}~{},$ (1)

where the cell is viewed as a pair of juxtaposed forces with dipole strength
$P$ and orientation $\hat{n}$. The dipole strength is a measure of the
anisotropy of the cell shape, and can be measured via a cell’s aspect ratio
Atia et al. (2018). Cell aspect ratio is important because it can influence
various cellular processes such as cell division Gibson et al. (2006),
migration Atia et al. (2018); Mitchel et al. (2020), and differentiation. The
aspect ratio is the ratio of the length of the longest axis of a cell to the
length of the shortest axis.

Figure 1:  Interaction energy $\varepsilon$ between a pair of dipoles. The
nature of the energy landscape depends on both the relative distance and the
relative orientation of the dipoles. In the two figures shown here, one dipole
is fixed at the origin with its axis aligned in the x direction. A second
dipole, either (A) parallel to the first dipole or (B) perpendicular to the
first dipole, is placed at different points in space and energy of the
configuration is indicated by the color at that point. Blue and yellow regions
represent attractive and repulsive regions respectively. (A) When both dipoles
are parallel, the most preferred configuration is the head-to-tail line up.
Stacking along the y-axis is also a preferred configuration. (B) However when
the dipoles are mutually perpendicular, both of the previously attractive
regions become repulsive. The new attractive regions lie along $45^{\circ}$ to
the x axis.

The interaction energy between cells can be written as:

$E=\sum_{i<j}P({\bf r}_{i})q_{\alpha\gamma}({\bf
r}_{i})G_{\alpha\beta,\gamma\delta}({\bf r}_{i}-{\bf r}_{j})P({\bf
r}_{j})q_{\beta\delta}({\bf r}_{j})~{},$ (2)

where $G_{\alpha\beta,\gamma\delta}({\bf r}_{i}-{\bf r}_{j})$ is the Green’s
function encoding the interaction between ith and jth dipole. Greek letters
denote tensor components and Einstein summation is implied over repeated Greek
indices. The elastic properties of the tissue are encoded in the Green’s
function,]. In our analysis, we use the 2D elastic Green’s function obtained
from the classic Boussinesq solution for an isotropic elastic medium Bischofs
et al. (2004b):

$G_{\alpha\beta}^{2D}(r)=a_{1}\left(a_{2}\delta_{\alpha\beta}+\frac{r_{\alpha}r_{\beta}}{r^{2}}\right)\frac{1}{r}$
(3)

where the parameters, $a_{1}$ and $a_{2}$, are defined in terms of the shear
modulus, $\mu$ and the Poisson ratio, $\nu$ as: $a_{1}=\nu/2\pi\mu$ and
$a_{2}=(1-\nu)/\nu$. In this study, lengths of the dipoles are considered to
be much smaller than the distance between dipoles. Therefore, this mechanistic
model does not incorporate any steric interactions or geometrical constraints
arising from the physical shapes of cells. A distinction that will be
important to bear in mind when comparing the results of this purely
mechanistic model to the Vertex model. The effective Green’s function between
two interacting dipoles separated by distance $r$ is

$\displaystyle
G_{\alpha\beta,\gamma\delta}(r)=a_{1}[\frac{1}{r^{3}}(a_{2}\delta_{\alpha\beta}\delta_{\gamma\delta}-\delta_{\alpha\delta}\delta_{\beta\gamma}-\delta_{\beta\delta}\delta_{\alpha\gamma})+\frac{3}{r^{5}}(\delta_{\alpha\gamma}r_{\beta}r_{\delta}+\delta_{\beta\gamma}r_{\alpha}r_{\delta}$
$\displaystyle+\delta_{\gamma\delta}r_{\alpha}r_{\beta}+\delta_{\alpha\delta}r_{\beta}r_{\gamma}+\delta_{\beta\delta}r_{\alpha}r_{\gamma}-a_{2}\delta_{\alpha\beta}r_{\gamma}r_{\delta})-\frac{1}{r^{7}}r_{\alpha}r_{\beta}r_{\gamma}r_{\delta}].$
(4)

The results presented in this paper are for an isotropic and incompressible
elastic medium, characterized by a Poisson ratio, $\nu=1/2$, and all energies
are measured in terms of the shear modulus, $\mu$. For convenience, we will
often use a dimensionless form of the energy given by
$\varepsilon=\frac{El^{3}}{N\langle P\rangle^{2}d^{2}}$, where $l$ is the
lattice distance, $N$ is total number of dipoles and $\langle P\rangle$ is the
mean value of polarization magnitude, and $d(\ll l)$ is the dipole length. It
is worth noting that this mechanistic model of cell interactions is inherently
tensorial: the interaction strength depends both on the relative positions and
on relative orientations of the dipole pair. In Fig. 1, we illustrate this
feature, for two representative configurations of force dipoles.

There is increasing evidence that liquid crystalline, nematic order, and
topological defects associated with such order, play a prominent role in
cellular development and function Armengol-Collado et al. (2022); Saw et al.
(2017); Mueller et al. (2019). This liquid crystalline order is associated
with the shape of cells, and is viewed as emerging from the collection of
cells in tissues behaving as an active liquid crystal. In general, migrating
cells tend to have a long axis, and the movement direction of neighboring
cells is strongly correlated Kawaguchi et al. (2017). For example, Saw et al
Saw et al. (2017). demonstrated that defects akin to those observed in nematic
liquid crystals manifest in epithelial tissues. The distribution of stress
around the defects is comparable to that observed in nematic liquid crystals,
and the occurrence of extrusions correlates closely with +1/2 defects. The
interaction potential in Eq. 2, derived from a mechanistic perspective of a
jammed solid, has commonalities with liquid-crystal models, and as we show in
this work, orientational order can naturally emerge from the interaction
between force-dipoles in a jammed solid. This model of force dipoles in a
solid has features that are distinct from an active liquid crystal model.
Before embarking on a discussion of our numerical simulations and results, we
explore these differences.

In liquid crystals, the microscopic degrees of freedom are the orientations,
$\hat{n}$, of the objects, whether they are ellipsoidal molecules or cells
with a longitudinal axis. Active liquid crystals, are driven out of
equilibrium by active stresses and fluid flow Mueller et al. (2019). Nematic
order, in both equilibrium and active liquid crystals, is characterized by the
traceless tensor Chaikin and Lubensky (1995):
$Q_{\alpha,\beta}=(\hat{n}_{\alpha}\hat{n}_{\beta}-\delta_{\alpha\beta}/2)$,
in 2D. Rewriting the force-dipole interaction in terms of $Q$, reveals the
qualitatively different features of the mechanistic model from the energy
function in liquid crystal models (see Supplementary Information for detailed
calculation):

$\displaystyle E$ $\displaystyle=\sum_{i<j}P({\bf r}_{i})Q_{\alpha\gamma}({\bf
r}_{i})G_{\alpha\beta,\gamma\delta}({\bf r}_{i}-{\bf r}_{j})P({\bf
r}_{j})Q_{\beta\delta}({\bf r}_{j})$ $\displaystyle+\sum_{i}P({\bf
r}_{i})Q_{\alpha\gamma}({\bf r}_{i})h_{\alpha\gamma}({\bf r}_{i})~{},$ (5)

where we have defined a local field, $h_{\alpha\gamma}({\bf
r}_{i})=\sum_{i<j}G_{\alpha\beta,\gamma\beta}({\bf r}_{i}-{\bf r}_{j})P({\bf
r}_{j})$. This representation highlights two important sources of difference
between the force dipole model and a generic liquid crystal model: (1) there
is an effective field that couples to the nematic field (second set of terms
in Eq. 5), and (2) the field $P(\bf{r})$, which is believed to be distributed
broadly.

A universal statistical distribution that governs the shapes of cells in cell
monolayers across various tissues and biological processes has recently been
discovered in research Atia et al. (2018). This distribution has been observed
in different situations, such as the development of the pseudostratified
bronchial epithelial layer in both asthmatic and non-asthmatic donors, and the
formation of the ventral furrow in the Drosophila embryo. Using the aspect
ratio (defined as the ratio between the long and short axes of the cell) to
quantify cell shape elongation, Atia et al Atia et al. (2018) showed that
epithelial cells inside confluent tissues obey a universal distribution that
is well described by a $k-\Gamma$ distribution Sadhukhan and Nandi (2022).
This naturally occuring hetereogenity form the basis for an inherent
polydispersity in our model.

Figure 2:  (A) Initial and (B) final stages of a Monte Carlo simulation of
dipoles on a square lattice. Arrows indicate orientation of dipoles, and
colorbar indicates magnitude of polarization. Polarization magnitudes are
randomly drawn from a gamma distribution ($k$=2.5, $\theta$=2.0). (A) In the
initial state dipoles have random orientations and there is no spatial
correlation between polarization magnitudes of neighboring dipoles. (B) The
final state obtained by slowly anenealing the system to a low temperature,
shows formation of domains of similar polarization magnitudes and domains with
non-zero nematic order. we should give a sense of how low the temperature is
compared to the energy of the initial states (C) During the anneal, the energy
of the system decreases sharply and soon the system gets stuck in a metastable
state. (D) Average nematic order parameter $S$ shows a steady shift from zero
to $\pm 1$ indicating emergence of nematic order. Here $\phi$ is the
orientation angle of each dipole. (E) The weighted order parameter $\Omega$ is
plotted against time. In the disordered state the value of this weighted
nematic order parameter is close to zero, but as the system gets ordered its
value becomes more positive or negative. In this example, this order parameter
approaches unity. Simulation parameters :
$k=2.5,\theta=2.0,l=0.1,L=1.5,T=0.001,\nu=0.5,\mu=0.5$.

An interesting feature, which we will explore in our numerical studies, is the
possibility of different dynamics associated with the fields $P$ and
$Q_{\alpha,\beta}$. If the two fields relax on the same time scales, then we
have a system with annealed disorder arising from the polydispersity. In this
case, a generalized nematic order parameter can be defined as
$p_{\alpha,\beta}\equiv\langle PQ_{\alpha,\beta}\rangle$. This is similar to
the generalized shape function defined in Huang et al. (2022); Fielding et al.
(2022); Hernandez and Marchetti (2021); Armengol-Collado et al. (2022) We
refer to this dynamics as Model 1, and unless otherwise stated, the results in
the main text are all obtained from numerical simulations of this model. The
configurations of Model 1 can be fully characterized by the order parameter,
$p_{\alpha,\beta}$. This model is not “frustrated” in the classic sense since
the interactions are specified just by the elastic Green’s function. The local
field, $h_{\alpha\gamma}({\bf r}_{i})$, is strongly influenced by the total
magnitude of the polarization, $\sum_{i}P(\bf{r}_{i})$, and as we will show,
numerically, the energy of the configurations in Model 1 is primarily governed
by the total polarization magnitude of the sample. Since this quantity depends
on the distribution of $P(\bf{r})$, the energetics also depends on the
distribution that characterizes a particular tissue line. We will present
results, both for the nature of orientational order, and the characteristics
of the energy landscape for different realizations of the distribution of
$P(\bf{r})$.

If the field $P$ evolves on a much longer time scale, for example, because the
distribution of $P$ is fixed, then we have a situation of quenched disorder.
In this case, the effective interaction is a quenched random variable:
$J_{ij}=P({\bf r}_{i})G({\bf r}_{i}-{\bf r}_{j})P({\bf r}_{j})$. This model
can be frustrated if the effective interaction $J_{ij}$ is incompatible with
long range order, which could lead to the appearance of defects in the nematic
field. This scenario is closer to nematic ordering in porous media such as
aerogels Mertelj and Copic (1997), if we threshold $P$ and envision regions
with $P$ lower than the threshold as being “pores”. In the mechanistic model,
the local field, $h_{\alpha,\beta}$, is also a quenched variable, and this is
a difference from classical nematic models in disordered media. We refer to
this quenched-disorder model as Model 2.

In this paper, we focus primarily on the results of Model 1, obtained from
Monte Carlo simulations. We will also present analysis of data obtained from
vertex-model simulations of sheared tissues. We will briefly discuss results
from Model 2 in the Supplementary Information.

##  4 Numerical results

We present numerical results from Monte Carlo simulations of our mechanistic
model. Our primary objective is to understand the interplay between the
spatial organization of the magnitude of the force dipoles (cell aspect
ratios) and their orientations. We also explore the energy landscape and
reveal important connections with magnitudes of force dipoles. We compare
results from Monte Carlo simulations with those from Vertex model simulations.

###  4.1 Monte Carlo Simulation Results

Figure 3:  Snapshots of final states obtained from Monte Carlo annealing runs
starting with three different initial realizations of the distribution of $P$.
(a) $k=2.5$, $\theta=3.0$, (b) $k=2.5$, $\theta=2.0$, (c) $k=2.5$,
$\theta=1.0$. All three columns (A-C) of each row (a,b,c) shows the same
configuration through different lenses. Column (A) shows heatmap of
polarization magnitudes and orientations of dipoles. Column (B) represent the
nematic scalar order parameter $S=\cos 2\phi$ as heatmap. The heatmap in
column (C) shows the weighted nematic order parameter $\Omega=P\cos
2\phi/\langle P\rangle$. In columns (B) and (C), the polarization magnitudes
are represented by contour lines. It is evident from column (C) that the
domains of large polarizations align with largest absolute values of $\Omega$
and hence are coincident with ordered domains. Simulation parameters:
$\nu=0.5$, $\mu=0.5$, $L=1.5$, $l=0.1$, $T=0.001$ .

We perform Monte Carlo simulation to explore the energy landscape and the
spatial organization of the force dipoles emerging from the interactions
between force dipoles (Eq. 2). We consider a square lattice of size $L\times
L$ where each lattice site ($r_{i}$) contains a point force dipole
characterized by a polarization magnitude $P(r_{i})$ and an orientation angle
$\phi_{i}$. Several biological cell lines are known to have their aspect
ratios distributed according to a $k-\Gamma$ distribution. We replicate this
in our simulations by drawing polarization magnitudes ($P$) from a
distribution characterized by parameters $k$ and $\theta$, and described by
probability density function

$f(P)=\frac{1}{\Gamma(k)\theta^{k}}P^{k-1}e^{-P/\theta}$ (6)

The orientation angles ($\phi$) are randomly chosen between $0$ and $\pi$. The
interaction between force dipoles is given by equation 2, where a dipole’s
nematic tensor $q_{\alpha\beta}(r_{i})=\hat{n}_{\alpha}\hat{n}_{\beta}$ can be
obtained by writing the dipole orientation vector as
$\hat{n}=(\cos\phi,\sin\phi)$. The total energy is obtained by summing the
long range interactions between all pairs of dipoles with periodic boundary
condition on all sides.

To analyze the energy landscape, we access a multitude of metastable states by
thermalizing the system at a high temperature and then running the simulation
at a temperature that is low compared to the energy of the system:
$\frac{k_{B}T}{E}\approx 10^{-6}$ (Fig. 7). During each Monte Carlo step, the
two properties of each dipole, polarization and orientation, are updated
following two different rules. For the orientation, we employ model A (non-
conserved) dynamics [Hohenberg and Halperin (1977)], where a new orientation
angle is proposed randomly between 0 and $\pi$. The polarization magnitudes
are updated following model B (conserved) dynamics [Hohenberg and Halperin
(1977)], where the polarization magnitudes of two randomly chosen dipoles are
exchanged, so that the overall distribution of the polarization does not
change during the simulation, however, the spatial organization evolves along
with the orientation of the dipoles. Acceptance rate of a proposed update is
determined by the Boltzmann factor $\exp(-\frac{\Delta E}{K_{B}T})$, $\Delta
E$ being the difference of energy of the proposed state and the present state,
and $T$ is the temperature.

The simulation is run for a long enough time to ensure that the system reaches
a metastable, local minimum, of the energy function: the energy fluctuates
around a well-defined average value. The time evolution of the system is
monitored by observing three different quantities - (i) the total energy of
the system, $E$, which undergoes a sharp decrease from the initial high-energy
state, and then fluctuates around a much lower value, (ii) the average nematic
scalar order parameter, $S=\langle\cos 2\phi\rangle$, which evolves from
$\approx 0$ in the initial state to a value close to unity in the low energy
state, (iii) the average weighted order parameter, $\Omega=\langle P_{i}\cos
2\phi\rangle/\langle P\rangle$, which shows a similar trend, growing from
$\approx 0$ in the initial state to a large positive value ($\geq 1$) in the
final state. The time evolution of these three quantities, for a sample run,
are plotted in fig 2.

###  4.2 Orientational Order

Figure 2 shows initial and final configurations of a sample run. The initial
state (A) shows dipoles that are placed on the sites of the square lattice,
with random orientations and polarization magnitudes drawn from a $k-\Gamma$
distribution. The final state (B), which has a much lower energy, has a
majority of the dipoles aligned along the horizontal direction. We also
observe that the dipoles are sorted into well-defined domains by their
polarization magnitudes.

In fig 3, we take a closer look at the spatial organization of the order
parameters $S$ and $\Omega$ in the final, low-energy states. The three rows
represent results from simulations with three different values of $\theta$
(Eq. 6). Each row represents adifferent representation of the same end state,
in order to highlight different measures of the ordering. The leftmost figure
in each row shows the two defining properties of every dipole - orientation
(arrows) and polarization magnitude (heatmap). We see clear formation of
domains of large polarization in all three sets of data. The middle figure in
each row shows local values of scalar order parameter $S$ which again
organizes the system into domains of nematic order, oriented in different
directions. In the rightmost figure of each row, the heatmap shows local
values of the weighted nematic order parameter $\Omega$. From the figures it
is clear that $\Omega$ provides a better characterization of the “order” than
the pure nematic order parameter, $S$, and clearly identifies regions where
the nematically ordered domains overlap with domains of large polarization
magnitudes. For clarity we have also plotted the polarization magnitudes using
contour lines.

These figures clearly demonstrate, one of our primary findings, that domains
of largest polarization magnitude overlap with domains of highest nematic
order. Although not surprising in retrospect, this aspect of the mechanistic
model is missing from pure liquid crystalline models of tissues, and points to
the need for characterizing cell-polarity in tissues by both the magnitude of
the polarization and its orientation. Our results imply that because of
mechanistic interactions between cells, defects in orientational order would
likely be localized to regions where the cells are not significantly distorted
from an isotropic shape. The consequences of this for tissue development and
morphogenesis needs to be explored further.

Figure 4: Panel showing results from voronoi Vertex model simulations where a
constant shear in $xy$ direction is applied to a tissue which undergoes
multiple rearrangements in response to the external shear. Snapshots from the
vertex model simulations from (A, B) initial and (C, D) final time step. (A,
C) shows orientation (arrows) and magnitude of polarization (heatmap), and (B,
D) shows weighted nematic order parameter (heatmap) and polarization magnitude
(contours). In the initial state, the system is random with no spatial or
orientational order, but as time progresses we observe emergence of domains
with large nematic order which also coincides with domains of large
polarization. This matches with our observation from the Monte Carlo
simulations above. (E) Distribution of polarization changes throughout the
process unlike the Monte Carlo simulations where polarization distribution is
held constant throughout an annealing run. Both (F) nematic scalar order
parameter $S$ and (G) weighted nematic order parameter $\Omega$ goes from near
zero in the random phase to a large positive value in the ordered phase.

Therefore, further exploration of the mechanisms that regulate cell
orientation and the consequences of defects in cell orientation for tissue
development and disease is crucial. This could involve using advanced imaging
and computational techniques to study the dynamics of cell orientation in
tissues, as well as investigating the role of mechanical and biochemical
factors in regulating cell orientation. Ultimately, a better understanding of
these processes could lead to new therapeutic strategies for treating a wide
range of diseases and disorders.

In a separate set of simulations, of Model 2, we keep the dipoles fixed at the
same positions as in the initial disordered state. Dipoles do not exchange
polarization magnitudes and are allowed to update only their orientation
angles. Please see supplemental data for results.

The dynamics of cells within solid-like biological tissues is, of course, much
more complex than can be captured by Model 1 or 2, discussed in this work. To
connect to observations in a much more realistic model of tissues, we present
results from simulations of the vertex model, subjected to external shear.
Most biological tissues undergoing shape change due to development, growth,
morphogenesis, wound healing etc. are under external stresses. We have not
extended our mechanistic model to explore the effect of external stresses.
Below, we compare the order-parameter correlations that develop in a vertex
model under shear to results from Model 1.

###  4.3 Vertex Model Simulations

The Voronoi-based implementation Bi et al. (2016) of the vertex model
Farhadifar et al. (2007); Li et al. (2019, 2018); Yan and Bi (2019); Mitchel
et al. (2020); Das et al. (2020) is employed to model a 2D cell layer, with
the cell centers ${\mathbf{r_{i}}}$ serving as the degrees of freedom and
their Voronoi tessellation dictating the cellular structureBi et al. (2016).
The mechanics of the cell layer are described by the energy function Staple et
al. (2010)

$E=\sum_{i=1}^{N}\left[K_{A}(A_{i}-A_{0})^{2}+K_{P}(P_{i}-P_{0})^{2}\right]$
(7)

, where $N$ is the number of cells, $K_{A}$ and $K_{P}$ are the area and
perimeter elastic moduli, respectively, and $A_{i}$ and $P_{i}$ are the area
and perimeter of the $i$-th cell. $A_{0}$ and $P_{0}$ are their corresponding
equilibrium values. The origin of the first term in the expression, which is
quadratic in the cell areas ${A_{i}}$, can be attributed to the
incompressibility of the cell volume. As a result, a 2D area elasticity
constant $K_{A}$ and a preferred area $A_{0}$ are produced, as discussed in
Farhadifar et al. (2007); Staple et al. (2010). The second term in the
expression, which is quadratic in the cell perimeters ${P_{i}}$, is a result
of the contractile nature of the cell cortex and is described by an elastic
constant $K_{P}$ Farhadifar et al. (2007). The target cell perimeter $P_{0}$,
which represents the interfacial tension between adjacent cells arising from
the competition between cortical tension and adhesion Farhadifar et al.
(2007); Staple et al. (2010). Here, we focus on the case where all cells share
the same single cell parameters: $K_{A},K_{P},P_{0},A_{0}$. We also choose
$A_{0}=\bar{A}$, the mean cell area, which defines the unit of length. For all
results presented in this work, we used $N=400$ cells.

Figure 5:  Correlation of polarization magnitudes $C_{P}$ and that of order
parameter $C_{\Omega}$ as a function of distance $r$. The correlation is
measured along two directions, $\hat{r}_{\parallel}$, which is parallel to the
global orientation of the dipoles, and its perpendicular direction
$\hat{r}_{\perp}$, as indicated in figures (A, B). (C) and (E) show
correlation functions in the initial and final stages of the Monte Carlo
simulations, disorder averaged over 100 runs. The system is completely random
in the initial state (A), with the correlation length dying out just beyond
$r=0$. In the final state (E), we see emergence of long range order and this
is reflected in the correlation lengths which span almost half the system
size. Both correlations $C_{P}$ and $C_{\Omega}$ are much stronger in the
parallel direction than their counterparts in the perpendicular direction. (D)
and (F) shows correlations in the initial and final frames respectively of the
vertex model simulation, averaged over 84 runs. In the initial state (D), the
correlations are very short ranged. However in the final state (F), the
correlation lengths are longer. Like in the Monte Carlo case, the correlations
in the parallel direction are stronger than those in the perpendicular
direction. The Monte Carlo simulations (C, E) have the same parameters as the
vertex model simulations (D, F): $k=34$, $\theta=0.05$, $L=20$, average
nearest neighbor distance $l=1$, Poisson ratio $\nu=0.5$.

To study tissue mechanical response in the vertex model, we subject the model
tissue to quasistatic simple shear using Lees-Edwards boundary conditions
Allen and Tildesley (1989); Huang et al. (2022). We start from a strain-free
state, the strain $\gamma$ is increased in increments of $\Delta\gamma=2\times
10^{-3}$, while cells are subject to an affine deformation given by
$\Delta\mathbf{r_{i}}=\Delta\gamma\ y_{i}\ \hat{x}$. At each strain step, the
tissue energy is minimized to find a mechanically stable configuration.

We plot the data obtained from vertex model simulations in fig 4. Fig 4(A)
shows initial state ($\gamma=0$) of the tissue with the cellular aspect ratios
indicated by the heatmap and cell orientations indicated by arrows. Fig 4(B)
shows a configuration of the tissue in the quasistatic plastic flow regime at
larger strain values. Here the cell shapes are oriented at $\sim~{}45^{\circ}$
to the $x$ axis due to shear. We observe alignment of cells and formation of
domains of large polarization magnitudes parallel to the direction of the
external shear. Figs 4(C,D) shows the corresponding heat maps for the order
parameter $\Omega$, and contour lines indicating magnitude of polarization.
Here again we observe that the domains of large polarization coincides with
domains of large nematic order, similar to the Monte Carlo simulations. This
system reaches steady state following a very different mechanism from the
Monte Carlo simulations. Here the cells undergo elongation and plastic failure
Huang et al. (2022) multiple times as indicated in the time series plot of the
mean polarization magnitude. Still we find end states that are qualitatively
similar to those of annealed Monte Carlo runs. The emergence of order in the
system is indicated by the time series plots of $S$ and $\Omega$, both of
which starts at a value close to $0$ in the initial state and plateaus to a
large positive value as time progresses.

###  4.4 Comparison between Monte Carlo and Vertex model

Figure 6:  Energy versus average mean field. The $y$ axis represents
dimensionless energy in units of $\mu$, and $x$ axis represents components of
the dimensionless field tensor $h$. Each data point represents final time
frame of a different experiment. The color of each point is the $\theta$ value
of the corresponding polarization distribution. (A,B) Results from Monte
Carlo, (C,D) Vertex model simulations. Energy is calculated using eq 2 for
both models. In case of the vertex model, we have assumed the Poisson ratio of
the tissue to be $0.5$. The imposed shear in the vertex model, in the $xy$
direction creates a strong effective field in that direction, which is absent
in the Monte Carlo simulations without shear.

Both Monte Carlo and vertex model simulations show the common features of (i)
formation of ordered domains, (ii) formation of domains of large polarization,
and (iii) overlap of these two types of domains. In order to quantify these
observations we calculate spatial correlations of (a) the polarization
magnitudes and the (b) order parameter $\Omega$ for the two models (fig 5). In
the initial states, the correlation of both $P$ and $\Omega$ die out very
fast, indicating that there is no correlation in these variables to start
with. In the final state since we observe formation of anisotropic domains
with a long and a short axis , we calculate correlations in two independent
directions - (i) parallel to the largest domain and (ii) perpedicular to it.
Because of the externally imposed shear, the “parallel” direction is always at
$45^{\circ}$ to the $x$ axis in case of the vertex model simulation. In the
case of Monte Carlo simulations, for each final state, we align the largest
domain along the $x$ axis, for convenience. In the final state, both $C_{P}$
and $C_{\Omega}$ decay simultaneously indicating a strong correlation between
the two quantities. For the Monte Carlo runs we see a significant increase in
correlation lengths of both quantities to almost $1/5$-th of the system size
in perpendicular direction and system size in parallel direction. For the
vertex model simulations too, the correlations decay at the same rate for both
$P$ and $\Omega$ as a function of distance. The correlations in this case too
are stronger in the parallel direction than those in the perpendicular
direction. But the correlation lengths are systematically smaller in the vertx
model.

Figure 7:  (A) In this figure we plot the time dependence of energy of the
system for 10 independent MC runs. Each run is characterised by an effective
value of $k$ and $\theta$. Each system reaches a different steady state energy
depending on its specific $k,\theta$ value. (B) Distribution of normalized
energies $\varepsilon$ obtained from the last configuration of each Monte
Carlo run.

Lastly, we analyze the energy landscape explored at low temperatures by
looking at the energies of the final states in the simulations. In the Monte
Carlo simulations, independent runs achieve different values of the final
energy. We observe, however, that if the simulation is run with a fixed set of
$P$ but with different spatial realizations in the initial state, they attain
final states with energy values that are virtually indistinguishable from each
other. This observation indicates that the energy landscape at low
temperatures is not sensitive to the initial configurations but is controlled
by the particular realization of $P$, which are drawn from a $k-\Gamma$
distribution. To quantify this feature, we plot the energy as a function of
the average field $h_{\alpha\beta}$ (eqn 5) in the final states. Each point in
fig 6 represents one configuration, with its $y$ coordinate indicating total
energy and $x$ coordinate indicating value of average field
$h_{\alpha\gamma}=\frac{1}{N}\sum_{i}\sum_{i<j}G_{\alpha\beta,\gamma\beta}({\bf
r}_{i}-{\bf r}_{j})P({\bf r}_{j})$. We see a strong correlation between energy
and the average local field $h_{xx}+h_{yy}$ in both Monte Carlo and vertex
model simulations. As seen from Eq.5, the average local field is controlled by
the total polarization, $\sum_{i}P_{i}$. We have checked that replacing
average field $h_{\alpha\beta}$ by the total polarization produces plots that
capture the same correlation as seen in Fig. 6. The vertex model runs also
show a strong correlation between final energy and the transverse field
$h_{xy}$ because of the shearing in that direction, but this trend is absent
in the Monte Carlo simulations, where the values of $h_{xy}$ are also much
smaller because of the isotropy of the system.

###  4.5 Energy Landscape

Figure 8:  (A) Scatter plot of absolute values of final energy of each run
against the respective $k,\theta$ values obtained from the particular
distribution of P. The three groups of data correspond to simulations run at
three sets of parameters $k=2.5$ and $\theta=1.0,2.0,3.0$. The actual values
of $k,\theta$ in these finite size systems vary about these mean values. The
scaled energy values in all three groups of points have the same range. We
divide the data by the values of final scaled energy into three sections
$\varepsilon_{1},\varepsilon_{2},\varepsilon_{3}$, such that each section has
equal number of data points. (B) Schematic diagram showing
$\hat{r}_{\parallel}$ and $\hat{r}_{\perp}$ directions with respect to the
domain of largest polarization. (C) Correlation functions $C_{P}$ ad
$C_{\Omega}$ in the parallel direction averaged over points belonging to the
energy sectors $\varepsilon_{1}<\varepsilon_{2}<\varepsilon_{3}$. Both $p$ and
$\Omega$ have the longest correlation length corresponding to
$\varepsilon_{3}$, and shortest correlation length corresponding to
$\varepsilon_{1}$. (D) Correlation functions calculated in the same way in the
perpendicular direction. There is no significant effect of energy on the
correlation length in this direction.

The other aspect of the model that we explore in this study is the energy
landscape explored in Model 1. Fig 7(A) shows the time dependence of energy
for ten independent Monte Carlo runs. Each of our simulations is initiated
with a set of polarization magnitudes drawn from a $k-\Gamma$ distribution
which is then kept fixed throughout that run. The initial configuration is
completely random and hence is a high energy state. As the simulation
progresses, the system approaches a low energy configuration. It is evident
from the plot that each system approaches a different value of the final
energy. The energy of the final state depends on the average value of
polarization magnitudes of the dipoles. We performed three different sets of
simulations fixing the value of $k$ to be $2.5$ but with $\theta=1.0,2.0,3.0$.
Since the systems are of finite size, each set of polarization magnitudes
drawn from, nominally, the same $k-\Gamma$ distribution is characterized by
$k,\theta$ values that are slightly different, leading to different values of
the final energy.

If we normalize the final energies from all three sets of runs by the square
of the average polarization of the respective run, then they fall in a pretty
narrow range of values as shown by the distribution in fig 7(B).

In fig 8(A), we plot the normalized final energies of each system against
their respective $k,\theta$ values. In all three sets of data, the range of
the normalized energy is approximately the same. We divide the data points
into three sections by their energy values and plot correlations in two
mutually perpendicular directions $\hat{r}_{\parallel}$ and $\hat{r}_{\perp}$.
$\hat{r}_{\parallel}$ is defined as the direction in which the largest
polarization domain is aligned and correlation $C_{\parallel}$ is calculated
by taking pairs of points lying on thin strips along this direction. Similarly
for $C_{\perp}$ we take pairs of points on thin strips along the perpendicular
direction. We observe that both $P$ and $\Omega$ have longer-ranged
correlations in the parallel direction, as compared to the perpendicular one.
In addition, as we go to higher values of absolute energies (lower in the
energy landscape) the systems become more strongly correlated in the parallel
direction whereas energy has no significant effect on the correlations in the
perpendicular direction. Unsurprisingly, the conclusion is that lower energy
states are more ordered.

##  5 Conclusion

Using a combination of numerical simulations and theoretical analysis, our
study investigates the interplay between the spatial distributions of cell
aspect ratios and orientational (nematic) order, as well as the complexity of
the resulting energy landscape.

Our main finding reveals a strong correlation between nematic order and the
magnitude of cell polarizations (aspect ratios), which are known to follow
$k-\Gamma$ distributions that vary across different tissue types. Therefore,
it is critical to analyze the spatial organization of both cell aspect ratios
and orientations in order to fully understand tissue morphogenesis and
development.

These results have significant implications for the field of tissue
engineering and regenerative medicine, as they highlight the importance of
understanding how cell shape and orientation impact tissue function and
organization. By better understanding the complex interplay between these
factors, we can develop more effective strategies for engineering functional
tissues in vitro and for promoting tissue regeneration in vivo.

The mechanistic model used in the study has not been extended to explore the
effects of external stresses on tissue shape changes. This means that the
model does not take into account the influence of external factors that can
affect tissue shape and organization. In reality, the shape changes that occur
in biological tissues are often the result of complex interactions between
internal cellular processes and external mechanical and biochemical factors.
For example, during embryonic development, the shape and positioning of organs
are influenced by mechanical forces generated by the surrounding tissues and
organs, as well as by the biochemical signals that regulate cell behavior.

Therefore, it is important to extend the mechanistic model used in the study
to include the effects of external stresses on tissue shape changes. This
would require incorporating additional parameters and variables into the
model, such as the magnitude and direction of external forces, the stiffness
and elasticity of the surrounding tissues, and the presence of biochemical
signals that modulate cell behavior. Such an extended model would provide a
more realistic representation of tissue shape changes and could help to
uncover the underlying mechanisms that regulate these processes. This could
ultimately lead to a better understanding of how tissues develop, grow, and
repair, as well as how they respond to various physiological and pathological
conditions.

## Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any
commercial or financial relationships that could be construed as a potential
conflict of interest.

## Author Contributions

KM and BC conceptualized the project. KM, BC and DB designed the study. KM
designed the simulation protocol. KM and RR performed simulations. KM analyzed
data from Monte Carlo simulations and the Vertex model. KM, BC and DB wrote
different sections of the manuscript. All authors contributed to manuscript
revision, and have read and approved the submitted version.

## Funding

KM and BC acknowledge financial support from Brandeis Bioinspired MRSEC
through grant no : NSF-DMR 2011846. DB was supported in part by NSF
DMR-2046683, Alfred P. Sloan Foundation and The Human Frontier Science
Program.

## Acknowledgments

KM and BC acknowlege many useful discussions with Micheal D’Eon. The authors
would like to thank Anh Nguyen and Junxiang Huang for providing vertex model
simulation data. KM thanks Christopher Amey and Minu Varghese for useful
discussions.

## References

  * Allen and Tildesley (1989) Allen, M. P. and Tildesley, D. J. (1989).  _Computer Simulation of Liquids_ (New York, NY, USA: Clarendon Press)
  * Armengol-Collado et al. (2022) Armengol-Collado, J.-M., Carenza, L. N., Eckert, J., Krommydas, D., and Giomi, L. (2022).  _Epithelia are multiscale active liquid crystals_.  Tech. rep., bioRxiv.  10.1101/2022.02.01.478692.  Section: New Results Type: article
  * Atia et al. (2018) Atia, L., Bi, D., Sharma, Y., Mitchel, J. A., Gweon, B., A. Koehler, S., et al. (2018).  Geometric constraints during epithelial jamming.  _Nature Physics_ 14, 613–620.  10.1038/s41567-018-0089-9.  Number: 6 Publisher: Nature Publishing Group
  * Bi et al. (2016) Bi, D., Yang, X., Marchetti, M. C., and Manning, M. L. (2016).  Motility-Driven Glass and Jamming Transitions in Biological Tissues.  _Physical Review X_ 6, 021011.  10.1103/PhysRevX.6.021011.  Publisher: American Physical Society
  * Bischofs et al. (2004a) Bischofs, I. B., Safran, S. A., and Schwarz, U. S. (2004a).  Elastic interactions of active cells with soft materials.  _Physical Review E_ 69, 021911.  10.1103/PhysRevE.69.021911.  Publisher: American Physical Society
  * Bischofs et al. (2004b) Bischofs, I. B., Safran, S. A., and Schwarz, U. S. (2004b).  Elastic interactions of active cells with soft materials.  _Phys. Rev. E_ 69, 021911.  10.1103/PhysRevE.69.021911
  * Cates et al. (1998) Cates, M. E., Wittmer, J. P., Bouchaud, J.-P., and Claudin, P. (1998).  Jamming, force chains, and fragile matter.  _Phys. Rev. Lett._ 81, 1841–1844.  10.1103/PhysRevLett.81.1841
  * Chaikin and Lubensky (1995) Chaikin, P. M. and Lubensky, T. C. (1995).  _Principles of Condensed Matter Physics_ (Cambridge Universisty Press)
  * Das et al. (2020) [Dataset] Das, A., Sastry, S., and Bi, D. (2020).  Controlled neighbor exchanges drive glassy behavior, intermittency and cell streaming in epithelial tissues.  10.1101/2020.02.28.970541.  Pages: 2020.02.28.970541 Section: New Results
  * Das et al. (2021) Das, A., Sastry, S., and Bi, D. (2021).  Controlled Neighbor Exchanges Drive Glassy Behavior, Intermittency, and Cell Streaming in Epithelial Tissues.  _Physical Review X_ 11, 041037.  10.1103/PhysRevX.11.041037.  Publisher: American Physical Society
  * Farhadifar et al. (2007) Farhadifar, R., Röper, J.-C., Aigouy, B., Eaton, S., and Jülicher, F. (2007).  The influence of cell mechanics, cell-cell interactions, and proliferation on epithelial packing.  _Current biology: CB_ 17, 2095–2104.  10.1016/j.cub.2007.11.049
  * Fielding et al. (2022) [Dataset] Fielding, S. M., Cochran, J. O., Huang, J., Bi, D., and Marchetti, M. C. (2022).  Constitutive model for the rheology of biological tissue
  * Gibson et al. (2006) Gibson, M. C., Patel, A. B., Nagpal, R., and Perrimon, N. (2006).  The emergence of geometric order in proliferating metazoan epithelia.  _Nature_ 442, 1038–1041
  * Hernandez and Marchetti (2021) Hernandez, A. and Marchetti, M. C. (2021).  Poisson-bracket formulation of the dynamics of fluids of deformable particles.  _Phys. Rev. E_ 103, 032612.  10.1103/PhysRevE.103.032612
  * Hohenberg and Halperin (1977) Hohenberg, P. C. and Halperin, B. I. (1977).  Theory of dynamic critical phenomena.  _Reviews of Modern Physics_ 49, 435–479.  10.1103/RevModPhys.49.435.  Publisher: American Physical Society
  * Huang et al. (2022) Huang, J., Cochran, J. O., Fielding, S. M., Marchetti, M. C., and Bi, D. (2022).  Shear-driven solidification and nonlinear elasticity in epithelial tissues.  _Phys. Rev. Lett._ 128, 178001.  10.1103/PhysRevLett.128.178001
  * Kawaguchi et al. (2017) Kawaguchi, K., Kageyama, R., and Sano, M. (2017).  Topological defects control collective dynamics in neural progenitor cell cultures.  _Nature_ 545, 327–331
  * Li et al. (2018) Li, X., Das, A., and Bi, D. (2018).  Biological tissue-inspired tunable photonic fluid.  _Proceedings of the National Academy of Sciences_ 115, 6650–6655.  10.1073/pnas.1715810115.  Publisher: Proceedings of the National Academy of Sciences
  * Li et al. (2019) Li, X., Das, A., and Bi, D. (2019).  Mechanical Heterogeneity in Tissues Promotes Rigidity and Controls Cellular Invasion.  _Physical Review Letters_ 123, 058101.  10.1103/PhysRevLett.123.058101.  Publisher: American Physical Society
  * Mertelj and Copic (1997) Mertelj, A. and Copic, M. (1997).  Evidence of dynamic long-range correlations in a nematic-liquid-crystal–aerogel system.  _Phys. Rev. E_ 55, 504–507.  10.1103/PhysRevE.55.504
  * Mitchel et al. (2020) Mitchel, J. A., Das, A., O’Sullivan, M. J., Stancil, I. T., DeCamp, S. J., Koehler, S., et al. (2020).  In primary airway epithelial cells, the unjamming transition is distinct from the epithelial-to-mesenchymal transition.  _Nature Communications_ 11, 5053.  10.1038/s41467-020-18841-7.  Number: 1 Publisher: Nature Publishing Group
  * Mueller et al. (2019) Mueller, R., Yeomans, J. M., and Doostmohammadi, A. (2019).  Emergence of active nematic behavior in monolayers of isotropic cells.  _Phys. Rev. Lett._ 122, 048004.  10.1103/PhysRevLett.122.048004
  * Nampoothiri et al. (2022) Nampoothiri, J. N., D’Eon, M., Ramola, K., Chakraborty, B., and Bhattacharjee, S. (2022).  Tensor electromagnetism and emergent elasticity in jammed solids.  _Phys Rev E_ 106, 065004.  10.1103/PhysRevE.106.065004
  * Nampoothiri et al. (2020) Nampoothiri, J. N., Wang, Y., Ramola, K., Zhang, J., Bhattacharjee, S., and Chakraborty, B. (2020).  Emergent elasticity in amorphous solids.  _Physical review letters_ 125, 118002
  * Sadhukhan and Nandi (2022) Sadhukhan, S. and Nandi, S. K. (2022).  On the origin of universal cell shape variability in confluent epithelial monolayers.  _eLife_ 11, e76406.  10.7554/eLife.76406.  Publisher: eLife Sciences Publications, Ltd
  * Saw et al. (2017) Saw, T. B., Doostmohammadi, A., Nier, V., Kocgozlu, L., Thampi, S., Toyama, Y., et al. (2017).  Topological defects in epithelia govern cell death and extrusion.  _Nature_ 544, 212–216.  10.1038/nature21718.  Number: 7649 Publisher: Nature Publishing Group
  * Schwarz and Safran (2013) Schwarz, U. S. and Safran, S. A. (2013).  Physics of adherent cells.  _Reviews of Modern Physics_ 85, 1327–1381.  10.1103/RevModPhys.85.1327.  Publisher: American Physical Society
  * Staple et al. (2010) Staple, D. B., Farhadifar, R., Röper, J. C., Aigouy, B., Eaton, S., and Jülicher, F. (2010).  Mechanics and remodelling of cell packings in epithelia.  _The European Physical Journal E_ 33, 117–127.  10.1140/epje/i2010-10677-0
  * Vinutha et al. (2023) Vinutha, H. A., Diaz Ruiz, F. D., Mao, X., Chakraborty, B., and Del Gado, E. (2023).  Stress-stress correlations reveal force chains in gels.  _J Chem Phys_ 158, 114104.  10.1063/5.0131473
  * Yan and Bi (2019) Yan, L. and Bi, D. (2019).  Multicellular Rosettes Drive Fluid-solid Transition in Epithelial Tissues.  _Physical Review X_ 9, 011029.  10.1103/PhysRevX.9.011029.  Publisher: American Physical Society

## Supplementary Data

### Derivation of equation 5 in the main text

Total energy of a system is given by

$\displaystyle E$
$\displaystyle=\sum_{i<j}P(r_{i})q_{\alpha\gamma}(r_{i})G_{\alpha\beta,\gamma\delta}(r_{i}-r_{j})P(r_{j})q_{\beta\delta}(r_{j}).$
(8)

Rewriting $q_{\alpha\beta}=Q_{\alpha\beta}+\frac{\delta_{\alpha\beta}}{2}$ we
get

$\displaystyle E$
$\displaystyle=\sum_{i<j}P(r_{i})[Q_{\alpha\gamma}(r_{i})+\frac{\delta_{\alpha\gamma}}{2}]G_{\alpha\beta,\gamma\delta}(r_{i}-r_{j})P(r_{j})[Q_{\beta\delta}(r_{j})+\frac{\delta_{\beta\delta}}{2}]$
(9)
$\displaystyle=\sum_{i<j}P(r_{i})Q_{\alpha\gamma}(r_{i})G_{\alpha\beta,\gamma\delta}(r_{i}-r_{j})P(r_{j})Q_{\beta\delta}(r_{j})$
(10)
$\displaystyle+P(r_{i})\frac{\delta_{\alpha\gamma}}{2}G_{\alpha\beta,\gamma\delta}(r_{i}-r_{j})P(r_{j})\frac{\delta_{\beta\delta}}{2}$
(11)
$\displaystyle+P(r_{i})Q_{\alpha\gamma}(r_{i})G_{\alpha\beta,\gamma\delta}(r_{i}-r_{j})P(r_{j})\frac{\delta_{\beta\delta}}{2}$
(12)
$\displaystyle+P(r_{i})\frac{\delta_{\alpha\gamma}}{2}G_{\alpha\beta,\gamma\delta}(r_{i}-r_{j})P(r_{j})Q_{\beta\delta}(r_{j})$
(13)

The sum of the the third and the fourth term is

$\displaystyle T_{3}+T_{4}$
$\displaystyle=P(r_{i})Q_{\alpha\gamma}(r_{i})G_{\alpha\beta,\gamma\delta}(r_{i}-r_{j})P(r_{j})\frac{\delta_{\beta\delta}}{2}+P(r_{i})\frac{\delta_{\alpha\gamma}}{2}G_{\alpha\beta,\gamma\delta}(r_{i}-r_{j})P(r_{j})Q_{\beta\delta}(r_{j})$
(14)
$\displaystyle=\frac{1}{2}\left[P(r_{i})Q_{\alpha\gamma}(r_{i})G_{\alpha\beta,\gamma\beta}(r_{i}-r_{j})P(r_{j})+P(r_{j})Q_{\alpha\gamma}(r_{j})G_{\alpha\beta,\alpha\delta}(r_{j}-r_{i})P(r_{i})\right]\text{change
indices}$ (15)
$\displaystyle=P(r_{i})Q_{\alpha\gamma}(r_{i})G_{\alpha\beta,\gamma\beta}(r_{i}-r_{j})P(r_{j})$
(16)

The energy expression in 9 without the constant term is, therefore:

$\displaystyle
E=\sum_{i<j}P(r_{i})Q_{\alpha\gamma}(r_{i})G_{\alpha\beta,\gamma\delta}(r_{i}-r_{j})P(r_{j})Q_{\beta\delta}(r_{j})$
(17)
$\displaystyle+\sum_{i}P(r_{i})Q_{\alpha\gamma}(r_{i})\sum_{i<j}G_{\alpha\beta,\gamma\beta}(r_{i}-r_{j})P(r_{j})$
(18)

The above expression for the energy shows that the last term can be expressed
in terms of a field coupling to the order parameter, as discussed in the main
text.

$\displaystyle
E=\sum_{i<j}P(r_{i})Q_{\alpha\gamma}(r_{i})G_{\alpha\beta,\gamma\delta}(r_{i}-r_{j})P(r_{j})Q_{\beta\delta}(r_{j})$
(19)
$\displaystyle+\sum_{i}P(r_{i})Q_{\alpha\gamma}(r_{i})h_{\alpha\gamma}(r_{i})$
(20)

### Results of Model 2

Fig 9 shows the end result of Model 2 simulations (see main text). In this
model, the dipoles are fixed in space and are allowed to only update their
orientations using model A dynamics. Due to the absence of the exchnage
dynamics, dipoles of large polarization magnitudes cannot cluster together
like in Model 1. Instead we find small domains of local orientational order.

Figure 9:  Results from Model 2. Rows (a) and (b) shows end state of two
different runs with parameters (a)$k=2.5,\theta=1.0$, and (b)
$k=2.5,\theta=3.0$. The first figure (A) of each panel shows the two
properties of a dipole - polarization (heatmap) and orientation (arrows). The
second figures (B) show local values of $S$, and the third figures (C) show
local values of $\Omega$, In both (B) and (C), polarization magnitudes are
indicated by contours. Simulation parameters :
$L=1.5,l=0.1,T=0.001,\nu=0.5,\mu=0.5$.

### Calculation of energy and order parameter in vertex model

From the vertex model simulations we obtain four quantities for each dipole at
every time step - x and y position coordinates, orientation angle $\phi$, and
aspect ratio. From this information we calculate total energy $E$ of the
system at each time step using eq 2 (in the main text), where $P_{i}$ is the
cellular aspect ratio, $q_{\alpha\beta}=\hat{n}_{\alpha}\hat{n}_{\beta}$ is
the nematic tensor, $\hat{n}=(\cos\phi,\sin\phi)$ describes the unit vector of
the dipole. Tissues are incompressible and hence we assume their Poisson ratio
$\nu=0.5$ for calculating energy. To get the dimensionless form of the energy
we multiply by $l^{3}/N\langle P\rangle^{2}d^{2}$, where $l$ is the average
nearest neighbor distance, $N$ is the number of cells, $\langle P\rangle$ is
the average polarization and the $d$ is the dipole length which is taken as an
arbitrary small number ($d\ll l$).

In the vertex model simulations, a shear is applied in the $xy$ direction,
resulting in the cells to be aligned in that direction. In order to calculate
order parameter $\Omega$ we need to measure orientation of cells with respect
to the global orientation $\theta_{g}$ which can be obtained by diagonalizing
the global nematic tensor

$\displaystyle
Q_{g}=\frac{1}{N}\sum_{i}\begin{pmatrix}\cos^{2}\theta_{i}-1/2&\cos\theta_{i}\sin\theta_{i}\\\
\cos\theta_{i}\sin\theta_{i}&\sin^{2}\theta_{i}-1/2\end{pmatrix}.$ (21)

The largest eigenvalue ($\Lambda$) of $Q_{g}$ is related to the global
orientation $\theta_{g}$ as $\Lambda=\cos\theta_{g}\sin\theta_{g}$. A few
lines of algebra leads us to the relation
$\cos^{2}\theta_{g}=1/2\left(1\pm\sqrt{1-4\Lambda^{2}}\right)$. The local
nematic scalar order parameter

$\displaystyle S_{i}$ $\displaystyle=\cos 2(\theta_{i}-\theta_{g})$
$\displaystyle=2\cos^{2}(\theta_{i}-\theta_{g})-1$
$\displaystyle=2(\cos^{2}\theta_{g}\cos^{2}\theta_{i}+\sin^{2}\theta_{g}\sin^{2}\theta_{i}+2\cos\theta_{g}\sin\theta_{g}\cos\theta_{i}\sin\theta_{i})-1$
$\displaystyle=2(\cos^{2}\theta_{g}\cos
2\theta_{i}+\sin^{2}\theta_{i}+\lambda\sin 2\theta_{i})-1,$

and similarly $\Omega_{i}=\frac{P_{i}}{\langle
P\rangle}[2(\cos^{2}\theta_{g}\cos 2\theta_{i}+\sin^{2}\theta_{i}+\lambda\sin
2\theta_{i})-1]$.