|  Active nematics with deformable particles†
---|---
|  Ioannis Hadjifrangiskou, Liam J. Ruske, and Julia M. Yeomans
|  The hydrodynamic theory of active nematics has been often used to describe
the spatio-temporal dynamics of cell flows and motile topological defects
within soft confluent tissues. Those theories, however, often rely on the
assumption that tissues consist of cells with a fixed, anisotropic shape and
do not resolve dynamical cell shape changes due to flow gradients. In this
paper we extend the continuum theory of active nematics to include cell shape
deformability. We find that circular cells in tissues must generate sufficient
active stress to overcome an elastic barrier to deforming their shape in order
to drive tissue-scale flows. Above this threshold the systems enter a
dynamical steady-state with regions of elongated cells and strong flows
coexisting with quiescent regions of isotropic cells.

##

††footnotetext: The Rudolf Peierls Centre for Theoretical Physics, Beecroft
Building, Parks Road, Oxford, OX1 3PU, UK. E-mail:
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
##  1 Introduction

Living systems are inherently active as they use chemical energy from their
surroundings to do work. It is becoming increasingly clear that many features
seen in cell layers and tissues can be understood using the continuum theory
of active nematics, which describes hydrodynamic interactions between active
anisotropic particles. Examples include biofilm initiation, topological
defects in cell monolayers and epithelial expansion.1, 2, 3, 4, 5, 6, 7, 8, 9

While the theory of active nematics predicts many qualitative features such as
short-range orientational order, active turbulence and motile topological
defects,10, 11, 12, 13, 14 it is based on liquid crystal hydrodynamics which
assumes that the active particles are nematogens with a fixed aspect ratio.15
While this assumption may be a reasonable description for systems consisting
of rod-shaped cells such as fibroblasts or Escherichia coli,16 many particles,
such as MDCK cells or soft colloids, can undergo large shape changes and their
aspect ratio can vary significantly due to active forces.17, 18, 19
Simulations have shown that intercellular stresses in monolayers are enhanced
by cell deformation, creating a positive feedback loop that affects the
collective behaviour of the layer.20

In this paper we extend the continuum equations which describe active nematic
fluids by considering not only the magnitude and direction of the nematic
order, but also the shape dynamics of the underlying particles. We consider an
equation of motion for the aspect ratio of the deformable nematogens which
incorporates flow-driven stretching and compression of particles, as well as
elastic restoring forces. The aspect ratio of particles in turn affects the
orientational dynamics by modifying thermodynamic interactions between
particles and the way they align in shear flows.

In section 2 we present the equations of motion for particle shapes,
characterised by the aspect ratio $\omega$, the nematic order tensor
$\mathbf{Q}$ of the particle orientational distribution, and the associated
velocity field $\mathbf{u}$. The following sections, 3.1 and 3.2, are
analytical and numerical investigations of systems consisting of elastic,
active particles which are isotropic in the absence of flow. We report a shape
instability in which extensile active stress generated by particles drives the
formation of regions with highly anisotropic particles with nematic order.
Finally, in section 4 we summarize our results.

##  2 Equations of motion

We introduce coarse-grained equations of motion to describe the hydrodynamics
of elastic, nematic particles. To this end we approximate the shape of
particles as ellipsoidal, and define a field $\omega>1$ to describe their
local aspect ratio, Fig. 1 (a). The time evolution of the shape distribution
follows 21, 22, 23

$\displaystyle\left(\partial_{t}+\mathbf{u}\cdot\mathbf{\nabla}\right)\omega-2\omega
E_{\parallel}=-\Gamma_{\omega}\dfrac{\delta\mathcal{F}}{\delta\omega}.$ (1)

The first term in Eq. (1) describes advection by a macroscopic flow field
$\mathbf{u}$, while the second term accounts for shape changes due to flow
gradients. $E_{\parallel}=n_{i}E_{ij}n_{j}$ is the projection of the strain
rate tensor $E_{ij}=(\partial_{i}u_{j}+\partial_{j}u_{i})/2$ along the long
axis of the particles defined by the unit vector $\mathbf{n}$. Depending on
the particle orientation with respect to the extensional flow axis, particles
are either stretched or compressed by flows (Fig. 1 (b), (c)). The final term
in Eq. (1) models particle shape elasticity, where $\Gamma_{\omega}$ controls
the relaxation towards the minimum of a free energy $\cal{F}$ as defined
below.

Nematic ordering of anisotropic particles is described by a tensorial order
parameter, $Q_{ij}=3S/2\>\left(n_{i}n_{j}-\delta_{ij}/3\right)$, where $S$
quantifies the magnitude of the nematic alignment, and we restrict the
directors to lie within a plane so that $n_{z}=0$ and $u_{z}=0$. The time
evolution of $\mathbf{Q}$ follows the Beris-Edwards equation,24

$\left(\partial_{t}+\mathbf{u}\cdot\mathbf{\nabla}\right)\mathbf{Q}-\mathcal{W}=\Gamma_{LC}\mathbf{H},\>$
(2)

where $\Gamma_{LC}$ controls the relaxation towards an equilibrium quantified
by the molecular field
$\mathbf{H}=-\delta\mathcal{F}/\delta\mathbf{Q}+(\mathbf{I}/3)\text{Tr}(\delta\mathcal{F}/\delta\mathbf{Q})$.
Nematogens are not only advected by the fluid, but also rotated by gradients
in the flow field, which gives rise to the the co-rotational term 24

$\displaystyle\mathcal{W}_{ij}=$ $\displaystyle\left(\xi
E_{ik}+\Omega_{ik}\right)\left(Q_{kj}+\frac{\delta_{kj}}{3}\right)+\left(Q_{ik}+\frac{\delta_{ik}}{3}\right)\left(\xi
E_{kj}-\Omega_{kj}\right)$
$\displaystyle-2\xi\left(Q_{ij}+\frac{\delta_{ij}}{3}\right)Q_{kl}W_{lk},$ (3)

where $\Omega_{ij}=(\partial_{j}u_{i}-\partial_{i}u_{j})/2$ is the vorticity
tensor, the antisymmetric part of the velocity gradient tensor
$W_{ij}=\partial_{i}u_{j}$. The flow-alignment parameter $\xi$ quantifies the
rotation of anisotropic particles in shear flows (Fig. 1(d)). Its value
depends on the shape of the particles and it is zero for isotropic
particles.25 Therefore we assume $\xi=\xi_{0}(\omega-1)$ where $\xi_{0}$ sets
the flow-alignment scale.

Fig. 1: Particle deformations driven by (a) free energy minimisation, (b) flow
stretching, (c) flow compression, (d) flow alignment. Fig. 2: Left: Phase
space trajectories of the evolution of $\omega$ subject to an extensional flow
at an angle $\theta$ to the director. Right: Rotation and squeezing of a
particle in an extensional flow.

The total free energy of the system,
$\mathcal{F}=\int(f_{\omega}+f_{LC})\>d\mathbf{V}$, consists of two
contributions. The first contribution is an elastic energy associated with
particle shape deformations,

$\displaystyle
f_{\omega}=A_{\omega}\frac{1}{2}(\omega-\omega_{0})^{2}+A_{\omega}^{\star}\frac{1}{4}(\omega-\omega_{0})^{4}\>,$
(4)

where $\omega_{0}$ is the particle aspect ratio at equilibrium and
$A_{\omega}$, $A_{\omega}^{\star}$ are elastic deformation parameters. The
second represents the liquid crystal free energy which emerges from particle-
particle interactions,

$\displaystyle f_{LC}=$
$\displaystyle\dfrac{A_{LC}}{2}\left(S_{0}^{2}\>\left(\frac{\omega-1}{\omega}\right)^{2}-\dfrac{2}{3}\mbox{tr}(\boldsymbol{Q}^{2})\right)^{2}+\frac{1}{2}K_{LC}\left(\mathbf{\nabla}\mathbf{Q}\right)^{2},$
(5)

where $K_{LC}$ penalizes distortions in the director field in the one-elastic-
constant approximation and $A_{LC}$ sets the bulk-energy scale. Highly
anisotropic particles tend to align nematically due to excluded volume effects
whereas there is no alignment interaction between circular particles. This is
reflected in the free energy which has a minimum at a reference value,
$S=S_{0}$, for $\omega\rightarrow\infty$ while favouring $S=0$ for $\omega=1$.

Finally, the velocity field $\mathbf{u}$ of the fluid obeys the incompressible
Navier-Stokes equations,

$\nabla\cdot\mathbf{u}=0\>,$ (6)
$\rho\left(\partial_{t}+\mathbf{u}\cdot\mathbf{\nabla}\right)\mathbf{u}=\mathbf{\nabla}\cdot\Pi,\>$
(7)

where $\rho$ is the fluid density and the stress tensor
$\Pi=\Pi^{passive}+\Pi^{active}$ includes passive and active contributions.
The passive stress tensor consists of viscous dissipation and elastic stress
arising from liquid crystal hydrodynamics,26

$\Pi^{viscous}=2\eta\mathbf{E}\>,$ (8) $\displaystyle\Pi^{elastic}=$
$\displaystyle-p\mathbf{I}-\xi[\mathbf{H}\tilde{\mathbf{Q}}+\tilde{\mathbf{Q}}\mathbf{H}-2\tilde{\mathbf{Q}}\mbox{tr}(\mathbf{Q}\mathbf{H})]+\mathbf{Q}\mathbf{H}$
$\displaystyle-\mathbf{H}\mathbf{Q}-\mathbf{\nabla}\mathbf{Q}\left(\frac{\partial
f_{LC}}{\partial(\mathbf{\nabla}\mathbf{Q})}\right),$ (9)

where $\eta$ is the viscosity, $p$ is the bulk pressure and
$\tilde{\mathbf{Q}}=\left(\mathbf{Q}+\frac{1}{3}\mathbf{I}\right)$. Active
dipolar forces produced by individual particles give rise to an active
stress,27

$\Pi^{active}=-\zeta\mathbf{Q}\>,$ (10)

where $\zeta$ quantifies the magnitude of active forces. $\zeta>0$,
corresponds to extensile activity, where fluid is pushed outwards from the
particles along their direction of elongation, and pulled inwards along the
perpendicular axis. $\zeta<0$ corresponds to contractile activity, where the
flow direction is reversed. We consider extensile systems in what follows.

We use a hybrid lattice Boltzmann–finite difference method to solve the
equations of motion (1, 2, 6, 7).28 The lattice size is $200\times 200$ with a
lattice spacing $\Delta x=1$, LB timestep $\Delta t=1$ and periodic boundary
conditions. Unless stated otherwise, the simulation parameters are
$\Gamma_{LC}=0.1,A_{LC}=0.1,S_{0}=1,K_{LC}=0.015,\xi_{0}=0,\Gamma_{\omega}=0.01,A_{\omega}=0.04,A_{\omega}^{\star}=0.003,\omega_{0}=1,\zeta=0.001,\rho=1,\eta=2/3$.
The relevant inverse time-scales are as follows: $\Gamma_{\omega}A_{\omega}$,
the rate of relaxation of the particles to their preferred shape,
$\Gamma_{LC}A_{LC}$, the rate of relaxation of the order parameter towards the
free energy minimum, and $\zeta/\eta$, the rate of active stresses injected
into the system. In the absence of flow gradients, the regime where
$\Gamma_{\omega}A_{\omega}\gg\Gamma_{LC}A_{LC}$ has nematic alignment
responding slowly to changes in shape, creating a time delay associated with
how particles will align after their shape changes. In the inverse case,
nematic alignment responds quickly to shape changes, which is the regime we
are considering. The relevant length-scales are the active length-scale
$L_{act}\sim\sqrt{K_{LC}/\zeta}$ and a passive length-scale
$L_{nem}\sim\sqrt{K_{LC}/A_{\omega}}$, which sets the scale of shape
gradients. Simulations are initialised with random noise in $\omega$ around a
finite value of $\omega=1.3$ and $S=\left(\omega-1\right)/\omega$. We use a
runtime of $2\times 10^{5}$ LB timesteps. Initial conditions become irrelevant
after a runtime of $\approx 10^{5}$ timesteps. If $\omega$ falls below unity,
we rotate the director field to point towards the extensile flow axis. This is
physically motivated by considering a particle that is instantaneously
circular; any subsequent elongation will occur along the local extensile axis
of the active flow field.

##  3 Results

To illustrate the coupled time evolution of particle orientation $\mathbf{n}$
and aspect ratio $\omega$ we consider a passive system in an external
extensional flow and solve Eq. (2) for the director field in the Ericksen-
Leslie limit (constant $S$) 24 for a finite flow-aligning parameter which
scales linearly with the aspect ratio, coupled to Eq. (1) with
$\Gamma_{\omega}=0$. Depending on the position in phase space, there are two
distinct re-orientation mechanisms: If the angle $\theta$ between the
compressive flow axis and particle orientation is sufficiently large,
particles undergo rotation where their shape stays anisotropic and their
orientation rotates until it aligns with the extensional flow axis (red
trajectory in Fig. 2). If $\theta$ is small, however, elongated particles
cannot rotate fast enough and get squeezed by the flow until they become
isotropic, $\omega=1$. Subsequently their orientation aligns along the
extensional flow axis, $Q_{ij}\sim E_{ij}$, and their aspect ratio increases
again (blue trajectory in Fig. 2).

###  3.1  Linear Stability Analysis

To investigate the behaviour of an active system consisting of nearly circular
particles, we perform a stability analysis around a quiescent isotropic system
by adding small perturbations to the aspect ratio of particles,
$\omega=1+\delta\omega$. The shape perturbation
$\delta\omega\sim\exp{(i\mathbf{q}\cdot\mathbf{r}+\lambda t)}$ is applied
along a wavevector $\mathbf{q}=q\left(\cos{\theta},\sin{\theta},0\right)$ and
we assume particles deform along their long axis. The flow fields which result
from activity give rise to the dispersion relation,

$\displaystyle
q^{2}\left\\{\left(\eta+\dfrac{3}{2}\Gamma_{\omega}K_{LC}\rho\right)\lambda+\Gamma_{\omega}A_{\omega}\eta-\dfrac{3}{4}\zeta\sin^{2}{2\theta}\right\\}$
$\displaystyle+\dfrac{3}{2}\Gamma_{\omega}K_{LC}\eta
q^{4}+\rho\lambda^{2}+\Gamma_{\omega}A_{\omega}\rho\lambda=0\>.$ (11)

Fig. 3: Dispersion relation of $\lambda(q)$ from Eq. (3.1) plotted for
$\theta=\pi/4$ for different values of {${\zeta,K_{LC}}$}. Red/Blue curves
correspond to activity values smaller/greater than the threshold activity,
$\zeta_{c}$ respectively. Fig. 4: Fully developed active turbulence with
deformable particles. (a) Aspect ratio, $\omega$ of particles. (b) Velocity
field, $\mathbf{u}$. Colour scale indicates the magnitude of the velocity. (c)
Nematic director field $\mathbf{n}$ and scalar order parameter $S$. Colour
scale indicates the magnitude of $S$. Red/blue symbols correspond to the core
of $\pm 1/2$ defects respectively. $+1/2$ defects are oriented such that the
tail is along the line shown.

For a solution with $\lambda>0$ to exist for some value of $q$, we require the
coefficient of the quadratic term to be negative. Thus the critical value of
the activity below which an isotropic system will be stable is given by:

$\displaystyle\zeta_{c}=\dfrac{4\Gamma_{\omega}A_{\omega}\eta}{3\sin^{2}{2\theta}}.$
(12)

$\zeta_{c}$ is clearly minimised when $\theta=\pi/4$, (Fig. 3). At this angle,
the parallel and perpendicular components of the flow work in tandem to create
the most efficient elongation of the particles, see Appendix A. Previously,
activity was found to create nematic order due to flow aligning effects as
long as $\xi\zeta>0$.29 Here we provide a mechanism where activity may induce
nematic order without the need for flow alignment. Starting from an isotropic
system of circular particles, activity causes elongation and then
thermodynamic forces described by the free energy may allow for nematic order
to be established.
The competition becomes apparent if we explicitly write down the wavevector
below which the instability grows:

$\displaystyle
q_{0}^{2}=\dfrac{3\zeta\sin^{2}{2\theta}-4\Gamma_{\omega}A_{\omega}\eta}{6K_{LC}\Gamma_{\omega}\eta}.$
(13)

This depends on a balance between the active length-scale,
$L_{act}\sim\sqrt{K_{LC}/\zeta}$, and the passive length-scale,
$L_{nem}\sim\sqrt{K_{LC}/A_{\omega}}$.

###  3.2 Active turbulence

For values of the extensile activity above which flows are maintained the
system enters an active turbulent regime. Fig. 4 shows the dynamical steady
state corresponding to fully developed active turbulence in the system of
deformable particles. At any given point in time the system is characterised
by distinct regions of elongated particles separated by areas where there is
little or no particle extension. There is strong nematic ordering in the
regions where the particles are extended, and it is possible to identify the
motile topological defects characteristic of active turbulence. In active
turbulence with particles of fixed length, nematic domains are broken up by
the active instability and by the passage of topological defects. Here there
is an additional mechanism which results in regions of circular particles.
Within a nematic domain the gradients in $\mathbf{Q}$ are small, and thus
active forces are diminished. Consequently, the magnitude of $E_{\parallel}$
drops, and it is no longer able to overcome the free energy cost of
deformation, Eq. (4). The particles within these domains will therefore begin
to contract. Fig. 5 shows an idealised sketch of this process.

Fig. 5: Creation and contraction of nematic domains. Nematic elasticity
smooths out gradients in the director field, resulting in a domain of aligned
directors. The active force vanishes and particles contract due to the elastic
free energy cost.

To investigate how physical parameters affect the dynamical state of the
system of deformable particles, we measure the time and space averaged aspect
ratio, $\langle\omega\rangle$ as the activity parameter $\zeta$, the elastic
deformation parameter $A_{\omega}$, the elastic constant, $K_{LC}$ and the
flow alignment parameter, $\xi_{0}$ are varied. In all cases, activity was
high enough to keep the system in active turbulence.

Fig. 6: Average aspect ratio as a function of (a) the activity, $\zeta$, (b)
elastic deformation parameter, $A_{\omega}$, (c) the elastic constant,
$K_{LC}$ and (d) the flow alignment parameter $\xi_{0}$. A linear relationship
using a least-squares fit has been plotted in each case. Error bars are given
by the standard deviation of the spatial average of the aspect ratio over
$10^{5}$ timesteps.

Fig. 6 (a) shows that $\langle\omega\rangle$ increases linearly with the
activity for extensile systems. In active turbulence where the evolution of
the flow field is dominated by driving from the active stress, the magnitude
of the velocity, and consequently the strain rate, increase with the activity
parameter, resulting in a higher degree of elongation. Conversely, Fig. 6 (b)
shows that $\langle\omega\rangle$ decreases linearly with $A_{\omega}$. As the
elastic energy cost of deforming particles increases, the average aspect ratio
will decrease. Next, as $K_{LC}$ is increased, gradients in the director field
are smoothed out, reducing the overall active stress and consequently reducing
elongation. This effect is relatively small compared to the other parameters,
as Fig. 6 (c) shows. Finally, $\langle\omega\rangle$ also increases linearly
with the value of the flow aligning parameter, $\xi_{0}$, shown in Fig. 6 (d).
As $\xi_{0}$ increases, particles will align more efficiently with the
extensional flow axis, thus increasing their elongation.
Next, we look at any emergent length-scales in the system. As we have already
shown, in a system of fully developed active turbulence, distinct domains of
elongated, anisotropic particles coexist with domains of isotropic particles.
To quantify the size of these domains, we calculate the normalised correlation
function defined as
$C_{\omega}\left(r\right)=\langle\left(\omega\left(0,t\right)-\bar{\omega}\left(t\right)\right)\left(\omega\left(\mathbf{r},t\right)-\bar{\omega}\left(t\right)\right)\rangle/\langle\left(\omega\left(0,t\right)-\bar{\omega}\left(t\right)\right)^{2}\rangle$
where $\bar{\omega}$ is the spatial average of the aspect ratio at each point
in time and $\langle\cdot\rangle$ denotes a space and time average over all
points separated by a distance $r$. Fig. 7 shows how $C_{\omega}$ varies with
distance, scaled by the active length-scale, $L_{act}\sim\sqrt{K_{LC}/\zeta}$
for different choices of the activity and elasticity parameters $\zeta$ and
$K_{LC}$. The collapse of the data to a single curve shows that the size of
the domains is governed by $L_{act}$, as in active turbulence with fixed
length nematogens, reflecting that flow gradients, which occur on this length
scale are responsible for the correlations in the aspect ratio of particles.

Fig. 7: Aspect ratio correlation, $C_{\omega}$ as a function of distance,
scaled by the active length-scale, $L_{act}=\sqrt{K_{LC}/\zeta}$. Different
coloured curves correspond to differing values of $L_{act}$, generated by
varying $\zeta$ and $K_{LC}$ independently in systems exhibiting fully
developed active turbulence.

We also measure the density of topological defects throughout the system as
the activity $\zeta$, the elastic constant $K_{LC}$ and the cell elasticity
$A_{\omega}$ are varied. To obtain this data we track defects 30, 31 and also
set a cut-off value $\omega=1.1$ below which topological defects cannot be
observed.

As Fig. 8 shows, the defect density scales as $(L^{*}_{act}L_{nem})^{-1}$
where the modified active length-scale is defined as
$L_{act}^{*}\sim\sqrt{K_{LC}/(\zeta-\zeta^{*})}$, where $\zeta^{*}$ is the
threshold activity required to keep the system in fully developed active
turbulence - this is taken to be the activity value below which circular cells
dominate the system and consequently, topological defects cannot exist. An
estimate for this threshold activity is given below.

Fig. 8: Average number of topological defects as a function of the inverse
modified active length-scale, $L^{*}_{act}$ multiplied by $L_{nem}$. The
plotted line shows a least-squares fit through all the data points. Data
points were collected by varying the activity parameter, $\zeta$ for different
realisations of the parameters $\\{K_{LC},A_{\omega}\\}$. Error bars are given
by the standard deviation of the spatial average of the defect number over
$10^{5}$ timesteps.

This is an empirical result and is a fundamentally different scaling than
fixed length active nematics where the defect density scales with the inverse
active length-scale squared, $L_{act}^{-2}=\zeta/K_{LC}$.32 Indeed it seems
very reasonable that other length-scales should be involved. One clear
difference from the usual active nematic models is that the active stresses
depend on the local magnitude of the order parameter and therefore vary across
the system, particularly at interfaces, and this will influence the defect
count. A plausible explanation for this scaling is as follows: $L^{*}_{act}$
plays a similar role to $L_{act}$ in traditional active nematics, with a
constant offset in the activity due to a minimum threshold required to achieve
active turbulence. Secondly, the deformation energy $A_{\omega}$ will tend to
reduce the magnitude of the aspect ratio, $\omega$, and therefore the order
parameter, $S$. In turn, this will reduce the effectiveness of the nematic
elasticity which depends on $S$. There is then a reduced, effective elastic
constant $\sim K_{LC}/A_{\omega}$ so the defect density scales with
$L_{nem}^{-1}\sim\sqrt{A_{\omega}/K_{LC}}$.

An estimate for the threshold activity $\zeta^{\ast}$ may be obtained by using
the velocity field around $+1/2$ defects 33 to calculate the strain rate,
$E_{\parallel}^{+}$. For a director given by
$\mathbf{n}=\left(\cos{\theta/2},\sin{\theta/2},0\right)$ and assuming
$S=\left(\omega-1\right)/\omega$ outside the core, the strain rate is given in
polar coordinates by

$\displaystyle
E_{\parallel}^{+}(r,\theta)=\dfrac{\omega-1}{\omega}\dfrac{\zeta}{16\eta}\left(\cos{2\theta}+3\right).$
(14)

By substituting (14) into Eq. (1), neglecting advection and quartic
contributions in the free energy, a simple calculation yields a lower bound to
the activity required for particles to elongate along the axis of symmetry of
$+1/2$ defects:

$\displaystyle\zeta^{*}=2\Gamma_{\omega}A_{\omega}\eta.$ (15)

This estimate ignores elasticity: it is lower than the numerical values by
between 10% for low values of the elastic constant and 50% for high values of
the elastic constant.

##  4 Discussion

We have extended the continuum theory of active nematics to allow for
deformable particles. This was achieved by introducing an equation of motion
(1) for the aspect ratio $\omega$ of particles, which couples to the equations
for the nematic order parameter (2) and velocity fields (7) through the flow
alignment parameter $\xi$ and the nematic bulk free energy $f_{LC}$.

By analysing the linear stability of isotropic systems with initially circular
nematogens, we showed that they must overcome an activity threshold in order
to extend sufficiently to drive flows. Above this threshold the system enters
a dynamical steady state characterised by coexisting regions of elongated
particles, which tend to align nematically and are associated with topological
defects and spatiotemporal chaotic flows, and quiescent regions consisting of
isotropic particles. While isotropic particles get elongated by tissue-scale
flow gradients from neighbouring nematic regions, there is little extensional
flow deep in the nematic regions, where the aspect ratio of particles
subsequently decreases due to elastic shape energy. In steady-state, particles
continuously get stretched by flow gradients and subsequently contract in the
absence of flow gradients, leading to large variations of cell shape. Similar
cell shape variations have been observed in simulations of the vertex model
for unjammed tissues 34 and interpreted in terms of a continuum model.35 Mean
field theories of cell shape variations based on the vertex cellular Potts
models have also been presented.36, 37

We have focused on extensile materials. Contractile systems which are
initially circular will not show any dynamical behaviour because active forces
just reinforce the thermodynamic forces restoring the particles to circular.
We note, however, that our model is restricted to nematic driving: polar
forces which can lead to coherent motion of circular cells are also observed
in cell monolayers, and the resultant force on a cell does not necessarily act
along its nematic axis. 38, 39

Our model was motivated by the observation of active turbulence in confluent
cell layers where cells can be deformed. In the light of our results it would
be interesting to further investigate the spatial and temporal variation of
the aspect ratio of cells within confluent layers that are on average
circular, such as the MDCK cell line.20, 8

## Author Contributions

All authors worked on defining the project, analysing data and writing the
paper. I. H. ran the simulations.

## Conflicts of interest

There are no conflicts to declare.

## Acknowledgements

This project was funded by the Gould & Watson Scholarship and the European
Commission’s Horizon 2020 research and innovation programme under the Marie
Sklodowska-Curie grant agreement No 812780.

##  Appendix A Dispersion Relation Derivation

The inverse Fourier transform for a fluctuating field $f^{\prime}$ is given by

$\displaystyle f^{\prime}(\mathbf{r},t)=\int d\lambda
dq\tilde{f}(\mathbf{q},\lambda)e^{i\mathbf{q}\cdot\mathbf{r}+\lambda t}.$ (16)

The time evolution equation for the tensorial order parameter $\mathbf{Q}$ is
eliminated by setting $\mathbf{n}=\hat{x}$ and
$S=\left(\omega-1\right)/\omega$. The relevant equations of motion to first
order in the perturbed fields are then:

$\partial_{t}\omega^{\prime}=2E_{\parallel}^{\prime}-\Gamma_{\omega}A_{\omega}\omega^{\prime}+\dfrac{3}{2}\Gamma_{\omega}K_{LC}\nabla^{2}\omega^{\prime},$
(17)
$\rho\partial_{t}u_{i}^{\prime}=\partial_{j}\left(2\eta
E_{ij}^{\prime}-p^{\prime}\delta_{ij}-\dfrac{3}{2}\zeta\omega^{\prime}(n_{i}n_{j}-\dfrac{1}{3}\delta_{ij})\right).$
(18)

In Fourier space, incompressibility reads

$\displaystyle
q_{j}\tilde{u}_{j}=0\Rightarrow\tilde{u}_{x}=-\dfrac{\sin{\theta}}{\cos{\theta}}\tilde{u}_{y}.$
(19)

In Fourier space, Eq. (18) reads:

$\displaystyle\rho\lambda\tilde{u}_{i}=\eta\left(-q_{i}\underbrace{q_{j}\tilde{u}_{j}}_{\text{=
0}}-q^{2}\tilde{u}_{i}\right)-q_{i}\tilde{p}-\dfrac{3i}{2}\zeta\tilde{\omega}q_{j}\left(n_{i}n_{j}-\dfrac{1}{3}\delta_{ij}\right).$
(20)

Taking $\mathbf{n}=\left(1,0,0\right)$ and
$\mathbf{q}=q\left(\cos{\theta},\sin{\theta},0\right)$ where $\theta$ is the
angle between $\mathbf{n}$ and $\mathbf{q}$ and considering $i=x$, we get

$\displaystyle\rho\lambda\tilde{u}_{x}=-\eta
q^{2}\tilde{u}_{x}-iq\cos{\theta}\tilde{p}-i\zeta\tilde{\omega}q\cos{\theta},$
(21)

from which we can solve for $\tilde{p}$

$\displaystyle\tilde{p}=\dfrac{i}{q\cos{\theta}}\left(\rho\lambda+\eta
q^{2}\right)\tilde{u}_{x}-\zeta\tilde{\omega}=-\zeta\tilde{\omega}-\dfrac{i\sin{\theta}}{q\cos^{2}{\theta}}(\rho\lambda+\eta
q^{2})\tilde{u}_{y},$ (22)

where we have used incompressibility (19) in the final step. Now consider
$i=y$

$\displaystyle\rho\lambda\tilde{u}_{y}=-\eta
q^{2}\tilde{u}_{y}-iq\sin{\theta}\tilde{p}-\dfrac{3i}{2}\zeta\tilde{\omega}\left(-\dfrac{1}{3}q\sin{\theta}\right).$
(23)

Substituting the expression (22) for $\tilde{p}$

$\displaystyle\left(\rho\lambda+\eta
q^{2}\right)\tilde{u}_{y}=-\dfrac{\sin^{2}{\theta}}{\cos^{2}{\theta}}\left(\rho\lambda+\eta
q^{2}\right)\tilde{u}_{y}+iq\sin{\theta}\zeta\tilde{\omega}+\dfrac{i}{2}q\sin{\theta}\zeta\tilde{\omega},$
(24)

allows us to solve for $\tilde{u}_{y}$

$\displaystyle\tilde{u}_{y}=\dfrac{3iq\cos^{2}{\theta}\sin{\theta}\zeta\tilde{\omega}}{2\left(\rho\lambda+\eta
q^{2}\right)}.$ (25)

Next, we consider Eq. (17) in Fourier space:

$\displaystyle\lambda\tilde{\omega}=2iq\cos{\theta}\tilde{u}_{x}-\Gamma_{\omega}A_{\omega}\tilde{\omega}-\dfrac{3}{2}\Gamma_{\omega}K_{LC}q^{2}\tilde{\omega}.$
(26)

Substituting for $\tilde{u}_{x}$ by using (25) for $\tilde{u}_{y}$ and the
incompressibility relation (19), we finally get

$\displaystyle\lambda\tilde{\omega}=\dfrac{3q^{2}\cos^{2}{\theta}\sin^{2}{\theta}\zeta}{\left(\rho\lambda+\eta
q^{2}\right)}\tilde{\omega}-\Gamma_{\omega}A_{\omega}\tilde{\omega}-\dfrac{3}{2}\Gamma_{\omega}K_{LC}q^{2}\tilde{\omega},$
(27)

which, after some algebraic manipulation and a double angle formula gives the
dispersion relation

$\displaystyle
q^{2}\left(\left(\eta+\dfrac{3}{2}\Gamma_{\omega}K_{LC}\rho\right)\lambda+\Gamma_{\omega}A_{\omega}\eta-\dfrac{3}{4}\zeta\sin^{2}{2\theta}\right)$
(28) $\displaystyle+\dfrac{3}{2}\Gamma_{\omega}K_{LC}\eta
q^{4}+\rho\lambda^{2}+\Gamma_{\omega}A_{\omega}\rho\lambda=0.$ (29)

For a solution with $\lambda>0$ to exist for some value of $q$, we require the
coefficient of $q^{2}$ to be negative. There is then a critical value of the
activity, below which an isotropic system will be stable, given by

$\displaystyle\zeta_{c}=\dfrac{4\Gamma_{\omega}A_{\omega}\eta}{3\sin^{2}{2\theta}}.$
(30)

$\zeta_{c}$ is clearly minimised when $\theta=\pi/4$. An analytical expression
may be obtained for the most unstable wave vector, $q_{m}$ which corresponds
to the local maximum of the growth rate, $\lambda_{m}$. This is given by:

$\displaystyle q_{m}^{2}=\dfrac{\sqrt{2\rho
K_{LC}\Gamma_{\omega}\zeta\sin^{2}{2\theta}J_{+}^{2}R}-4\rho
K_{LC}\Gamma_{\omega}\eta T}{2K_{LC}\Gamma_{\omega}\eta J_{-}^{2}},$ (31)
$\displaystyle J_{\pm}=2\eta\pm 3\rho\Gamma_{\omega}K_{LC},$ (32)
$\displaystyle R=3\zeta\sin^{2}{2\theta}-4\Gamma_{\omega}A_{\omega}\eta+6\rho
A_{\omega}K_{LC}\Gamma_{\omega}^{2},$ (33) $\displaystyle
T=3\zeta\sin^{2}{2\theta}-\Gamma_{\omega}A_{\omega}J_{-}.$ (34)

Despite the complicated form of the above expression, the competition between
activity versus elasticity and energy is shown to be a determining factor for
the most unstable modes.
Finally, by setting $\omega^{\prime}=\epsilon\sin{\mathbf{q}\cdot\mathbf{r}}$
and using the Fourier transform of the $\sin$ function, we may also find a
solution of the flow field, $\tilde{u}_{y}$.

$\displaystyle\tilde{u}_{y}=\dfrac{3iq\cos^{2}{\theta}\sin{\theta}\zeta}{2(\rho\lambda+\eta
q^{2})}\dfrac{1}{2i}(2\pi)^{2}\epsilon\left(\delta(\mathbf{q}+\mathbf{q^{\prime}})-\delta(\mathbf{q}-\mathbf{q^{\prime}})\right).$
(35)

Setting $\lambda=0$ for simplicity, $\theta=\pi/4$, and inverting the Fourier
transform, we may get the full flow field:

$\displaystyle u_{y}=-\dfrac{3}{4\sqrt{2}}\dfrac{\zeta\epsilon}{\eta
q}\cos\left(\dfrac{1}{\sqrt{2}}(x+y)\right),$ (36) $\displaystyle
u_{x}=\dfrac{3}{4\sqrt{2}}\dfrac{\zeta\epsilon}{\eta
q}\cos\left(\dfrac{1}{\sqrt{2}}(x+y)\right).$ (37)

## Notes and references

  * Blanch-Mercader and Casademunt 2017 C. Blanch-Mercader and J. Casademunt, _Soft Matter_ , 2017, 13, 6913–6928.
  * Blanch-Mercader _et al._ 2018 C. Blanch-Mercader, V. Yashunsky, S. Garcia, G. Duclos, L. Giomi and P. Silberzan, _Phys. Rev. Lett._ , 2018, 120, 208101.
  * Saw _et al._ 2018 T. B. Saw, W. Xi, B. Ladoux and C. T. Lim, _Adv. Mat._ , 2018, 30, 1802579\.
  * Yaman _et al._ 2019 Y. I. Yaman, E. Demir, R. Vetter and A. Kocabas, _Nat. Comm._ , 2019, 10, 1–9.
  * Pearce _et al._ 2019 P. Pearce, B. Song, D. J. Skinner, R. Mok, R. Hartmann, P. K. Singh, H. Jeckel, J. S. Oishi, K. Drescher and J. Dunkel, _Phys. Rev. Lett._ , 2019, 123, 258101.
  * Duclos _et al._ 2017 G. Duclos, C. Erlenkämper, J.-F. Joanny and P. Silberzan, _Nat. Phys._ , 2017, 13, 58–62.
  * Pearce _et al._ 2021 D. Pearce, J. Nambisan, P. Ellis, A. Fernandez-Nieves and L. Giomi, _Phys. Rev. Lett._ , 2021, 127, 197801.
  * Saw _et al._ 2017 T. B. Saw, A. Doostmohammadi, V. Nier, L. Kocgozlu, S. Thampi, Y. Toyama, P. Marcq, C. T. Lim, J. M. Yeomans and B. Ladoux, _Nature_ , 2017, 544, 212–216.
  * Copenhagen _et al._ 2021 K. Copenhagen, R. Alert, N. S. Wingreen and J. W. Shaevitz, _Nat. Phys._ , 2021, 17, 211–215.
  * Ardaševa and Doostmohammadi 2022 A. Ardaševa and A. Doostmohammadi, _Nat. Rev. Phys._ , 2022, 4, 354–356.
  * Mueller _et al._ 2019 R. Mueller, J. M. Yeomans and A. Doostmohammadi, _Phys. Rev. Lett._ , 2019, 122, 048004.
  * Thampi and Yeomans 2016 S. P. Thampi and J. M. Yeomans, _Eur. Phys. J. Special Topics_ , 2016, 225, 651–662.
  * Alert _et al._ 2022 R. Alert, J. Casademunt and J. Joanny, _Annu. Rev. Condens. Matter Phys._ , 2022, 13, 143–170.
  * Doostmohammadi _et al._ 2016 A. Doostmohammadi, S. P. Thampi and J. M. Yeomans, _Phys. Rev. Lett._ , 2016, 117, 048102.
  * De Gennes and Prost 1993 P.-G. De Gennes and J. Prost, _The Physics of Liquid Crystals_ , Oxford University Press, 1993.
  * Dell’Arciprete _et al._ 2018 D. Dell’Arciprete, M. Blow, A. Brown, F. Farrell, J. S. Lintuvuori, A. McVey, D. Marenduzzo and W. C. Poon, _Nat. Comm._ , 2018, 9, 1–9.
  * Vlassopoulos and Cloitre 2014 D. Vlassopoulos and M. Cloitre, _Curr. Opin. Colloid Interface Sci._ , 2014, 19, 561–574.
  * Lecuit and Lenne 2007 T. Lecuit and P.-F. Lenne, _Nat. Rev. Molecular Cell Bio._ , 2007, 8, 633–644.
  * Heisenberg and Bellaïche 2013 C.-P. Heisenberg and Y. Bellaïche, _Cell_ , 2013, 153, 948–962.
  * Balasubramaniam _et al._ 2021 L. Balasubramaniam, A. Doostmohammadi, T. B. Saw, G. H. N. S. Narayana, R. Mueller, T. Dang, M. Thomas, S. Gupta, S. Sonam, A. S. Yap _et al._ , _Nat. Mat._ , 2021, 20, 1156–1166.
  * Bilby and Kolbuszewski 1977 B. A. Bilby and M. Kolbuszewski, _Proc. Roy. Soc. London. A._ , 1977, 355, 335–353.
  * Gao _et al._ 2011 T. Gao, H. H. Hu and P. P. Castañeda, _J. Fluid Mech._ , 2011, 687, 209–237.
  * Gao _et al._ 2013 T. Gao, H. H. Hu and P. P. Castañeda, _J. Fluid Mech._ , 2013, 715, 573–596.
  * Beris and Edwards 1994 A. N. Beris and B. J. Edwards, _Thermodynamics of flowing systems with internal microstructure_ , Oxford University Press, New York ; Oxford, 1994\.
  * Jeffery 1922 G. B. Jeffery, _Proc. Roy. Soc. London. A._ , 1922, 102, 161–179.
  * Marenduzzo _et al._ 2007 D. Marenduzzo, E. Orlandini and J. M. Yeomans, _Phys. Rev. Lett._ , 2007, 98, 118102.
  * Simha and Ramaswamy 2002 R. A. Simha and S. Ramaswamy, _Phys. Rev. Lett._ , 2002, 89, 058101\.
  * Marenduzzo _et al._ 2007 D. Marenduzzo, E. Orlandini, M. Cates and J. M. Yeomans, _Phys. Rev. E._ , 2007, 76, 031921.
  * Santhosh _et al._ 2020 S. Santhosh, M. R. Nejad, A. Doostmohammadi, J. M. Yeomans and S. P. Thampi, _J. Stat. Phys._ , 2020, 180, 699–709.
  * Blow _et al._ 2014 M. L. Blow, S. P. Thampi and J. M. Yeomans, _Phys. Rev. Lett._ , 2014, 113, 248303.
  * Vromans and Giomi 2016 A. J. Vromans and L. Giomi, _Soft Matter_ , 2016, 12, 6490–6495.
  * Giomi 2015 L. Giomi, _Phys. Rev. X_ , 2015, 5, 031003.
  * Giomi _et al._ 2014 L. Giomi, M. J. Bowick, P. Mishra, R. Sknepnek and M. Cristina Marchetti, _Philos. Trans. Roy. Soc., A_ , 2014, 372, 20130365.
  * Bi _et al._ 2016 D. Bi, X. Yang, M. C. Marchetti and M. L. Manning, _Phys. Rev. X_ , 2016, 6, 021011.
  * Hernandez and Marchetti 2021 A. Hernandez and M. C. Marchetti, _Phys. Rev. E_ , 2021, 103, 032612\.
  * Sadhukhan and Nandi 2022 S. Sadhukhan and S. Nandi, _Elife_ , 2022, 11, e76406.
  * Czajkowski _et al._ 2018 M. Czajkowski, D. Bi, M. L. Manning and M. C. Marchetti, _Soft Matter_ , 2018, 14, 5628–5642.
  * Ladoux and Mège 2017 B. Ladoux and R.-M. Mège, _Nat. Rev. Mol. Cell Bio._ , 2017, 18, 743–757.
  * Ron _et al._ 2023.07.02 J. E. Ron, J. d’Alesandro, V. Cellerin, R. Voituriez, B. Ladoux and N. S. Gov, _bioRxiv_ , 2023.07.02, 547405.