# Clustering and ordering in cell assemblies with generic asymmetric aligning interactions

Thibault Bertrand<sup>1,*</sup> Joseph d'Alessandro<sup>2</sup> Ananyo Maitra<sup>3,4</sup> Shreyansh Jain<sup>2</sup> Barbara Mercier,<sup>2</sup> René-Marc Mège<sup>2</sup> Benoit Ladoux<sup>2,†</sup> and Raphaël Voituriez<sup>4,5,‡</sup>

<sup>1</sup>Department of Mathematics, Imperial College London, South Kensington Campus, London SW7 2AZ, United Kingdom

<sup>2</sup>Université Paris Cité, CNRS, Institut Jacques Monod, F-75013 Paris, France

<sup>3</sup>Laboratoire de Physique Théorique et Modélisation, CNRS UMR 8089, CY Cergy Paris Université, F-95032 Cergy-Pontoise Cedex, France

<sup>4</sup>Laboratoire Jean Perrin, UMR 8237 CNRS, Sorbonne Université, 75005 Paris, France

<sup>5</sup>Laboratoire de Physique Théorique de la Matière Condensée, UMR 7600 CNRS, Sorbonne Université, 75005 Paris, France

<figure><img src="image_1.png" /></figure>

(Received 15 September 2021; accepted 21 February 2024; published 4 April 2024)

Collective cell migration plays an essential role in various biological processes, such as development or cancer proliferation. While cell-cell interactions are clearly key determinants of collective cell migration, the physical mechanisms that control the emergence of cell clustering and collective cell migration are still poorly understood. In particular, observations have shown that binary cell-cell collisions generally lead to antialignment of cell polarities and separation of pairs—a process called contact inhibition of locomotion (CIL), which is expected to disfavor the formation of large-scale cell clusters with coherent motion even though the latter is often observed in tissues. To solve this puzzle, we adopt a joint experimental and theoretical approach to determine the large-scale dynamics of cell assemblies from elementary pairwise cell-cell interaction rules. We quantify experimentally binary cell-cell interactions and show that they can be captured by a minimal equilibriumlike pairwise asymmetric aligning interaction potential that reproduces the CIL phenomenology. We identify its symmetry class, build the corresponding active hydrodynamic theory, and show on general grounds that such asymmetric aligning interaction destroys large-scale clustering and ordering, leading instead to a liquidlike microphase of cell clusters of finite size and short lived polarity or to a fully dispersed isotropic phase. Finally, this shows that CIL-like asymmetric interactions in cellular systems—or general active systems—control cluster sizes and polarity, and can prevent large-scale coarsening and long-range polarity, except in the singular regime of dense confluent systems.

DOI: 10.1103/PhysRevResearch.6.023022

## I. INTRODUCTION

The emergence of collective, coordinated migration is a striking property of eukaryotic cell collectives [1–4]. It is observed in key biological processes *in vivo* such as development [5,6], cancer proliferation [7,8], and wound healing [9], and has now been reproduced in various *in vitro* setups [10–17]. This ability of cells to form large-scale cohesive, polarized, self-propelled clusters is expected to be controlled both by single cell properties (polarity and motility) and cell-cell interactions, as is confirmed experimentally [15,18–22].

The effects of cell-cell interactions on collective cell dynamics are twofold. Upon contact, cells can engage trans-membrane adhesion molecules (such as cadherins) to form junctions [15]. On the one hand, these cell-cell junctions act as an effective attractive force that opposes the separation of cell pairs and therefore favors cell clustering. On

the other hand, it was found that cell-cell junctions also impact cell polarity; indeed, binary cell-cell interaction events are reported to typically favor outward pointing, antialigned polarities and, ultimately, separation of cell pairs. In the literature, this is generically known as contact inhibition of locomotion (CIL) [23], even though quantitative experimental analysis of the phenomenon remains sparse [24–27]. Thus far, the paradigm introduced by CIL, which favors antialignment and separation of pairs, therefore seems inconsistent with the observation of large-scale cell clusters with coherent motion [17,24,25,28–30]; reconciling these observations and more generally determining, from the knowledge of basic pairwise cell-cell interaction rules, the conditions of emergence of cell clustering and collective motion remains an outstanding open question.

Cell motility generically relies on the nonequilibrium dynamics of the actin/myosin system, driven by adenosine triphosphate (ATP) hydrolysis; from a physics standpoint, this makes the cell a prototypical self-propelled particle (SPP) and cell assemblies a striking example of active matter [31,32]. Active matter based models of collective cell migration, which involve—explicit or implicit—specific choices of cell-cell interaction rules have flourished [33]; these can take various forms, from agent-based models [21,28,29,34–38] and active vertex models [18,39,40] to active hydrodynamics models [16,41,42] and

Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.

<footnote><sup>*</sup>t.bertrand@imperial.ac.uk</footnote>

<footnote><sup>†</sup>benoit.ladoux@ijm.fr</footnote>

<footnote><sup>‡</sup>voiturie@lptmc.jussieu.fr</footnote>
phase fields models [43,44]. They point, mostly through
numerical simulations, to a broad variety of possible phases
that can help interpret experimental observations. In particu-
lar, agent-based models endowed with *ad hoc* interaction rules
aiming at mimicking the CIL phenomenology have been pro-
posed, and pointed to a rich phenomenology [24,25,28,29,45].

Even the simplest interactions between SPPs can, in fact,
have striking consequences at the collective level. For in-
stance, a simple pairwise aligning interaction between SPPs,
as introduced by Vicsek et al. [46], can lead to clustering
and large-scale collective motion in settings where long-range
order and phase separation would be forbidden for systems
at equilibrium [31,32]. Another purely nonequilibrium col-
lective effect is the propensity of SPPs to cluster or undergo
phase separation in the presence of purely repulsive interac-
tions [47,48]. A systematic exploration of the phase space of
possible behaviors for more realistic models of cell assemblies
with specific interaction rules is expected to lead to an ever-
increasing complexity and therefore seems inaccessible. In
this context, hydrodynamic theories, which are insensitive to
specific microscopic choices but governed by symmetry prop-
erties and conservation rules [31,32] are promising candidates
to provide unifying principles.

Here, we adopt a joint experimental and theoretical ap-
proach that integrates quantitative multiscale *in vitro* data,
numerical simulations and active hydrodynamic theory to
determine the conditions of emergence of clustering and
collective motion in cell assemblies. We make use of micro-
fabricated *in vitro* environments to quantitatively characterize
the onset of cell clustering and collective motion from the
scale of cell pairs to the scale of large aggregates in 1D and
2D. Rather than building a computational model based on
*ad hoc* interaction rules [28,29], we experimentally analyze
pairwise cell-cell interactions and show that the observed
CIL phenomenology can be rationalized by a minimal equi-
libriumlike asymmetric aligning interaction potential whose
symmetry class we identify. Based on experimental observa-
tions, we combine such an asymmetric aligning interaction
with a classical short-range attractive potential that mimics
cell-cell junctions, build the corresponding active agent-based
model and propose a minimal active hydrodynamic theory
of this symmetry class. We show both experimentally and
theoretically that the asymmetric aligning interaction can
drastically lower the persistence of finite cell clusters and
reduce their size. We demonstrate that in the large system limit
this can lead to a transition between a dispersed (gas) isotropic
phase and a liquidlike microphase of cell clusters of finite size
and short-lived polarity, which is critically controlled by both
the strength of the asymmetric interaction and the cells’ self-
propulsion force. Our results are applicable to general active
systems of the same symmetry class and show that CIL-like
interactions can regulate cluster sizes and polarity and, in par-
ticular, prevent large-scale coarsening and long-range polarity,
except in the singular regime of dense confluent systems.

## II. CIL AS A PAIRWISE EQUILIBRIUMLIKE INTERACTION POTENTIAL IN MDCK CELLS

We first aim at characterizing quantitatively pairwise cell-
cell interactions rules. Previous studies analyzed the scattering

rules of various cell types upon contact in similar 1D ge-
ometries [24–27]. These studies were mostly restricted to
the analysis of cell behaviors immediately following con-
tact, with the aim of determining scattering rules. Here, we
extended our analysis to longer timescales, while cells dy-
namically explored all available configurations, with the aim
of determining effective energy landscapes similar to those
of equilibrium interacting spins. We used MDCK (Madine-
Darby canine kidney) cells plated on fibronectin-coated linear
strips of width *w* = 20 µm, which we obtained by micro-
contact printing on polydimethylsiloxane (PDMS) [49]. Cells
were treated with mitomycin C to prevent cell division and
maintain cell numbers. In this geometry, cell motion is re-
stricted to a single dimension and cell-cell interactions are
limited to front-rear interactions. To study the dynamics of
front-rear cell polarization, we use a fluorescent biosensor
(p21-activated kinase binding domain, PBD) of active Rac1
and Cdc42 [50], whose gradient along the cell indicates the
direction of its polarity.

At the scale of a pair of cells, CIL dictates distinctive inter-
cellular dynamical rules that have been described qualitatively
in the literature [23,51]. Upon contact, cell-cell junctions are
formed and trigger mechanotransduction signals, leading to
the repolarization of the two cells away from the contact loca-
tion. To confirm this phenomenology, we observe the possible
outcomes of binary cell collisions. Prior to collision, incom-
ing cells show clear signs of polarization, with asymmetric
internal organization and shape, and a stable lamellipodium
at the leading edge [see Figs. 1(a) and 1(b)]. Upon contact
between the lamellipodia, we typically observe inversions of
the PBD gradient, signaling an inversion of cell polarity, with
two possible outcomes: (i) both cells repolarize away from
the contact [see Fig. 1(a)] or (ii) only one of the two cells
repolarizes away from the contact, leading to the alignment of
cell polarities [see Fig. 1(b)]. Because in this cell type cell-cell
adhesion is strong, even in case (i) the cells do not separate im-
mediately after repolarizing, but instead form a stable doublet
which can stay cohesive for several hours [Fig. 1(c)]. In 1D
geometries, the sign of the polarity gradient is used to map
this continuous measure of cellular polarity to a binary ±1
spinlike variable.

Following the motion and polarities of both cells of the
doublet in time, we notice that they stochastically switch
their polarity and therefore spontaneously explore all possi-
ble accessible configurations (see Fig. 1(d) and Supplemental
Material (SM) Fig. S1 and Movie S1 [52]). From several such
time series, we measure the respective probabilities of the
four possible configurations [Fig. 1(e)]. We observe a bias in
those probabilities: tail-tail configurations (↔) are strongly
favored compared to head-head configurations (→←), while
the frequency of tail-head configurations (←↔) is not raised.
Note that, as expected, symmetric configurations (←← and
→→) are equally probable. This quantitatively confirms the
presence of CIL in those cells, although it exhibits more com-
plex consequences than simple elastic rebounds upon cell-cell
contact, even at the scale of the doublet. In analogy with
equilibrium spin dynamics, we associate an energy level $E_{ij}$
to each of the polarity states of the cell couple {i, j}, so
its probability is given by $P_{ij} = \exp(-E_{ij})/Z$, where Z is a
normalization factor [Fig. 1(f)].
<figure><img src="image_2.png" /><figcaption>FIG. 1. Scattering rules in cell doublets. Potential outcomes of a collision between two cells: (a) contact inhibition of locomotion (CIL) and (b) cell alignment shown in transmission (left) and in PBD-YFP fluorescence (right); in both panels, time progresses from top to bottom with consecutive snapshots separated by $\Delta t = 10$ min, scale bar = 50 µm. (c) Position of a cell doublet as a function of time; the shaded grey areas correspond to moments where the cells are not in contact. (d) Polarity of the cell doublet at corresponding times (see SM, Fig. S1 [52]). (e) Probabilities for the four possible polarity configurations measured from a total of 3271 data points obtained from the analysis $n = 34$ cell doublets in two independent experiments; colors correspond to the highlighted regions in (c). (f) Schematic of the model with the energy levels $E_{ij}$ of the four possible doublet configurations with values of $\alpha/\beta \approx 2$ obtained from the probabilities in (e).</figcaption></figure>

### III. MINIMAL ACTIVE PARTICLE MODEL WITH ASYMMETRIC ALIGNING INTERACTIONS

Based on these quantitative observations, we build a model of active Brownian particles (ABPs) which includes minimal asymmetric interaction rules that describe the observed CIL phenomenology. The model is introduced in arbitrary space dimension $d$, but numerical simulations will be performed mostly for $d = 1$ to reproduce the experimental setup; extensions to $d = 2$ will be discussed in the last section. Each cell is described as a particle with position $\mathbf{r}_i$ and endowed with a polarity vector $\mathbf{p}_i$ (a unit vector). We start from an equilibrium description and assume that cell-cell interactions result from a microscopic Hamiltonian given by

$$ \mathcal{H} = \sum_{i,j} U_r(\mathbf{r}_i, \mathbf{r}_j) + U_p(\mathbf{r}_i, \mathbf{r}_j, \mathbf{p}_i, \mathbf{p}_j). \quad (1) $$

The position-dependent part that models cell-cell steric repulsion and adhesion can be taken to be the classical truncated Lennard-Jones potential,

$$ U_r = \begin{cases} 4\epsilon[(\sigma/r_{ij})^{12} - (\sigma/r_{ij})^6] & r_{ij} \le r_c \\ 0 & r_{ij} > r_c, \end{cases} \quad (2) $$

where $\mathbf{r}_{ij} = \mathbf{r}_i - \mathbf{r}_j$, $r_{ij} = |\mathbf{r}_{ij}|$, and $r_c$ defines the range of interaction and $\sigma$ the particle size. Our main conclusions will be independent of this specific choice, which is used in numerical simulations for convenience. To build the polarity interaction

potential $U_p$, we note that (for example, in $d = 2$) the CIL interaction explicitly breaks the invariance under independent rotations of space and polarity vectors, which is preserved by classical aligning interactions $\propto \mathbf{p}_i \cdot \mathbf{p}_j$ characteristic of XY models or their active counterparts, the class of Vicsek-like models [31,32,46,53,54]. More explicitly, for fixed positions $\mathbf{r}_i, \mathbf{r}_j$ of an interacting pair, the CIL phenomenology dictates that the system is not invariant under the symmetry $\mathbf{p}_i, \mathbf{p}_j \to -\mathbf{p}_i, -\mathbf{p}_j$ and, therefore, cannot simply be described via a potential $\propto \mathbf{p}_i \cdot \mathbf{p}_j$. Expanding in powers of $\mathbf{r}_{ij}$ and $\mathbf{p}_i, \mathbf{p}_j$, the simplest term that breaks the invariance under independent rotations of space and polarity vectors, while being invariant under their joint rotation, has the form $(\mathbf{p}_i - \mathbf{p}_j) \cdot \mathbf{r}_{ij}$. Without loss of generality, we therefore consider an interaction potential of the form

$$ U_p = \begin{cases} -\beta \mathbf{p}_i \cdot \mathbf{p}_j - \alpha (\mathbf{p}_i - \mathbf{p}_j) \cdot \mathbf{n}_{ij} & r_{ij} \le r_c \\ 0 & r_{ij} > r_c, \end{cases} \quad (3) $$

where $\mathbf{n}_{ij} = \mathbf{r}_{ij}/r_{ij}$. This interaction potential is composed of two terms: (1) a Vicsek or XY-like alignment term with interaction strength $\beta$, which, by construction, is invariant under independent rotations of space and polarity vectors and (2) an asymmetric alignment term of amplitude $\alpha$ which explicitly breaks this symmetry and reproduces the CIL phenomenology. Of note, potentials of the same symmetry class have been considered in the study of liquid crystals [55], but their effect in active matter systems have not been systematically
examined. From the measure of the probabilities for each cell doublet configuration [see Fig. 1(e)], we can estimate the value of the relative strength of the alignment interactions and find that our experiments lead to $\alpha/\beta \approx 2$, leading to the polarity energy levels shown in Fig. 1(f) (see Appendix C).

Finally, this so-far equilibrium description is made minimally active by assuming that particles are subject to the self-propelling force $F_p$ along their polarity $\mathbf{p}_i$. Our model can thus be interpreted as a generalization of the flying XY model [56], with an additional asymmetric interaction term that couples polarity and position and a short-scale steric repulsion. More explicitly, for $d \ge 2$, the dynamics of the system in the overdamped limit is governed by the set of coupled Langevin equations (see Appendix E for $d = 1$),

$$ \zeta \dot{\mathbf{r}}_i = -\frac{\partial \mathcal{H}}{\partial \mathbf{r}_i} + F_p \mathbf{p}_i + \sqrt{2T \zeta \eta_i}, \quad (4) $$

$$ \zeta_p \dot{\mathbf{p}}_i = (\mathbf{1} - \mathbf{p}_i \mathbf{p}_i^\top) \cdot \left[ -\frac{\partial \mathcal{H}}{\partial \mathbf{p}_i} + \sqrt{2T_p \zeta_p \xi_i} \right], \quad (5) $$

where $\zeta$ is the friction coefficient, $\zeta_p$ is the rotational viscosity, $T$, $T_p$ are the translational and polarity noise strengths (which can be distinct in out-of-equilibrium systems), and $\eta_i$ and $\xi_i$ are zero mean and unit variance Gaussian white noises. The projection operator $\mathbf{1} - \mathbf{p}_i \mathbf{p}_i^\top$ ensures that the magnitude of $\mathbf{p}_i$ remains invariant under the dynamics. It is useful to introduce the Péclet number $Pe = v_0 \sigma / D$, where $v_0 = F_p / \zeta$ is the self-propulsion velocity and $D = T / \zeta$ is the self-diffusion coefficient of the ABP. We will also make use of $\gamma = \varepsilon / T$ as the ratio of the strength of the Lennard-Jones potential to the thermal fluctuations, and similarly define $\alpha, \beta$ in units of $T_p$; finally, we introduce the normalized relaxation rate $\mu = \tau / \zeta_p$ for the polarity, where $\tau = \sigma^2 / D$. Earlier agent-based models [24,25,28,29,57] that take into account CIL interactions or comparable position-dependent orientational interactions [58] can be checked retrospectively to fall within this symmetry class, even though their specific choice of dynamics cannot be re-expressed as deriving from a simple pairwise effective potential. In particular, because in our model cell-cell interactions derive from an effective Hamiltonian, interactions are reciprocal, as opposed to earlier approaches [58].

The model is thus primarily controlled by (i) the volume fraction of particles $\phi$, (ii) the competition between self-propulsion (Pe) and cohesion ($\gamma$), and (iii) the strength of the symmetric and antisymmetric alignment interaction terms $\alpha, \beta$, and (iv) the relaxation rate $\mu$. Given this relative complexity, an exhaustive exploration of the phase behavior of this model goes beyond the scope of this paper. Below, we primarily aim to discuss the effect of the asymmetric coupling $\alpha$ on the collective particle dynamics, and we restrict our analysis to regimes that are most relevant to our experimental cellular system.

## IV. DYNAMICS OF SMALL CELL TRAINS

We first focus on the effect of the asymmetric interaction (parametrized by $\alpha$) on finite-sized cell clusters (or “trains”), based on our one-dimensional (1D) setup. As shown in Figs. 2(a) and 2(b) (and Movies S2 and S3), we observe that cells at the edges of the cell trains generically have opposite

polarities, pointing away from the center of mass of the train; this is expected from the CIL phenomenology—as reported in Fig. 1—which favors $\leftarrow \rightarrow$ configurations. This is evidenced by the extension of lamellipodia (see Movie S3) and the PBD gradients [see Fig. 2(b)]. In a given cell train, we generally observe a single domain wall where the polarity changes signs ($\leftarrow \rightarrow$). This behavior can be simply accounted for by the interaction potential $U_p$ introduced in Eq. (3). For $\alpha > \beta$, as observed experimentally [see Figs. 1(d) and 1(e)], the potential $U_p$ for a 1D train of $N$ particles is minimized for all configurations with a single domain wall $\leftarrow \rightarrow$; in particular, inducing such domain wall in a fully polarized configuration leads to an energy gain $\Delta E = -2(\alpha - \beta)$. There is no energy cost incurred when the domain wall moves one step to the left or right in the bulk of the cell train (corresponding to a polarization flip of a single particle); this suggests that domain walls perform symmetric random walks and thus diffuse within the train. This directly results from the fact that the asymmetric interaction term in $U_p$ reduces to a boundary term, as is evident from summing the interaction potential over all particles of a finite 1D train of $N$ particles,

$$ \sum_{i,j} U_p = -\beta \sum_{i=1}^{N-1} p_i p_{i+1} + \alpha(p_1 - p_N), \quad (6) $$

where $p_i = \pm 1$ in 1D. This shows that edges of cell clusters induce domain walls because of the CIL interactions. We argue below that these CIL-induced domain walls have important consequences on the dynamics of cell trains.

First, Fig. 2(b) shows that domain walls can lead to the fragmentation of clusters, as expected from the outward pointing polarity and therefore propulsion force of cells on each side of the domain wall. Shortly after the separation, we observe that the cells at the newly formed edges repolarize away from the center of mass of their respective clusters, thereby inducing domain walls in the new clusters, confirming the above scenario of domain-wall nucleation at cluster edges. This mechanism suggests that even for large cohesive interactions ($\gamma \gg 1$, which would lead in equilibrium to clusters whose size diverges for $T \to 0$ in 1D), the asymmetric interaction can induce fragmentation of large clusters into smaller, finite-sized clusters in active systems. From force balance, we infer that only clusters of typical size $N \lesssim \gamma / Pe$ are insensitive to this active fragmentation mechanism. While extensive statistics of cluster sizes and controlled tuning of Pe and $\gamma$ are not accessible experimentally, we have verified that this scaling correctly predicts the average cluster size in numerical simulations of the 1D version of the model [see Fig. 2(d)]. Note that, as we discuss below, for increasing values of Pe, a competing dynamic coarsening mechanism induced by the self-propulsion of clusters occurs, and the proposed scaling is insufficient to capture the average cluster size. Finally, this shows quantitatively that the asymmetric interaction that we introduced can lead in active systems to a drastic reduction of cluster sizes, which is finite and critically controlled by Pe.

Second, we now argue that domain walls in finite clusters control their dynamics and, in particular, their self-propulsion speed and persistence. Defining the train polarity per cell as $p(t) = \frac{1}{N} \sum_i p_i(t)$, we find from force balance that a train of
<figure><img src="image_3.png" /><figcaption>FIG. 2. Dynamics of small trains of cells. (a) Cohesive cell triplet showing antialigned polarities at the edges (arrows denote polarity and are provided as a guide)—the center cell flips polarity in the middle, slowing down the train. Snapshots are separated by $\Delta t = 120$ min, scale bar = 100 µm (see Movie S2). (b) Cell train containing eight cells shown in PBD-YFP fluorescence with repolarization of edge cells after fracture of the train; cell boundaries have been highlighted for clarity $\Delta t = 85$ min, scale bar = 100 µm (see Movie S3). In (a) and (b), the stars represent the location of the domain wall. (c) Examples of ring geometries at confluence (arrows show steady rotation direction for each ring); rings have diameters $D = 100$ µm, 200 µm, and 400 µm (scale bar = 200 µm—see Movie S4). (d) Steady-state average cluster size $\langle s_c \rangle$ as a function of the Péclet number $Pe$ for various cohesiveness $\gamma = 10$ (circles), 20 (squares), and 50 (diamonds) in the regime $\alpha \gg \beta$. Inset shows that $\langle s_c \rangle \propto (Pe/\gamma)^{-1}$. (e) Position of cells in two trains for closed boundaries (ring geometry) and open boundaries (line geometry), showing high persistence in the case of closed boundaries. (f) Experimental velocity autocorrelation for open and closed boundaries. Represented are mean±STD obtained from the analysis of $n = 5$ and $n = 1$ nine-cell trains, respectively, for open and closed boundaries. (g) Velocity correlation time (measured as the time for which correlation function reaches $1/e$) as a function of train size fitted by a quadratic law (red line). We analyzed data from $n = 5, 8, 9, 1, 8, 4, 8, 9, 10, 4, 5, 5, 2, 2, 5$ trains for $N = 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14$, and $N \ge 15$-cell trains, respectively, from two independent experiments. See Appendix B for further details about the determination of error bars. (h) Cell train global polarity autocorrelation function for various train sizes (increasing from blue to yellow: $N \in [3, 100]$) for open boundaries (for cohesive trains) and closed boundaries (horizontal grey line) averaged over 50 realizations. (i) Correlation time obtained from exponential fits of the autocorrelation function as a function of cell numbers $N$; the red solid line shows the theoretical prediction given by $\tau_p \sim N^2$.</figcaption></figure>

$N$ cells is expected from the model to be self-propelled with (dimensionless) velocity $p(t)$Pe. In turn, considering a train of $N$ cells with a single domain wall ($\alpha \gg \beta \gg 1$ regime), the sign of $p(t)$ is determined by the relative position of the domain wall to the train center. Since domain walls diffuse in the model, we find that, starting from a random position in a train, a domain wall reaches the train center and thus induces a velocity sign change with a mean time $\sim N^2$. We therefore expect that the polarity or velocity autocorrelation decays with a characteristic time $\tau_p \sim N^2$; this is indeed clearly observed in both experimental data and numerical simulations [see Figs. 2(g) and 2(i)].

Together, this shows that the CIL-based asymmetric interactions have striking consequences in active systems on size selection and dynamic properties of finite trains, as suggested in Refs. [28,29]. In contrast with these previous studies, we are here able to capture this phenomenology in a minimal equilibriumlike polarity interaction potential respecting the symmetry class of CIL-induced cell doublet configurations, which allows us to derive a formal hydrodynamic description of the model (see below). We argue that size selection and dynamic properties of finite trains are direct consequences of

the nucleation of domain walls at the edges of cell trains. A very simple consequence of this analysis is that clusters with no edges, for example, for dense confluent systems in closed periodic geometries, should be completely insensitive to the CIL asymmetric interaction. Strikingly, this is what we observed both experimentally and numerically when we analyzed periodic geometries [see Figs. 2(c), 3(a)–3(d), and Movie S4]: for comparable system sizes, the persistence time was found in periodic geometries with cells at confluence to be significantly larger (larger than the observation time) than for trains with edges in open geometries [see Figs. 2(e) and 2(f)]. Indeed we observed that when the system reaches confluence, cluster edges disappear, domain walls vanish, and sustained collective motion arises as predicted [see Figs. 3(a)–3(d)]. This is consistent with the persistent collective rotational motion of confluent clusters in ring geometries, reported in Ref. [17]. In this study, focusing on closed confluent systems, however, the impact of CIL was not discussed. In the following, we describe the route to collective motion in this system and solve the apparent paradox between locally prevailing CIL and globally emerging coherent motion.
<figure><img src="image_4.png" /><figcaption>FIG. 3. Cell ordering and clustering in ring geometries. (a)–(d) Example snapshots of systems with increasing number of cells in a ring geometry with fixed diameter *D* = 400 µm; nuclei are fluorescently labeled, black arrows show the orientation of the velocity of the associated cell nucleus. As the number of cells increases (here, from 3 to 28 cells), the system reaches confluence; at confluence, cells coordinate their polarity and rotate in the same direction (scale bar = 100 µm). (e) A kymograph of the orthoradial velocity *v*θ (measured by PIV) as a function of time *t* and angular position *θ* shows a transition to persistent rotational motion. The ring diameter is *D* = 1000 µm and cell number *N* = 91. (f) Clustering index *d*<sub>g</sub> − *d*<sub>l</sub> measured at long times *t* > 48 h as a function of the number of cells in the ring geometry with diameter *D* = 1000 µm (blue circles). The data was binned and interpolated to produce the red line as a visual guide. The system transitions from a dispersed phase (*d*<sub>g</sub> − *d*<sub>l</sub> < 0) to a clustered phase (*d*<sub>g</sub> − *d*<sub>l</sub> > 0) for *N*<sup>*</sup> ≈ 20, the vertical dashed line corresponds to the number of cells needed to reach confluence (see SM, Fig. S3 [52]). Data obtained from the analysis of *n* = 94 ring patterns from three independent experiments. (g) Steady-state polar order parameter ⟨|*p*|⟩ averaged over realizations and time as a function of the number of cells in the ring geometry with diameter *D* = 400 µm (error bars are given by standard deviations); the red line shows the expected value of the order parameter for a random set of *N* polarities (see Appendix D). Data obtained from the analysis of *n* = 189 ring patterns from two independent experiments.</figcaption></figure>

V. IMPACT OF ASYMMETRIC INTERACTIONS
ON COLLECTIVE CELL BEHAVIOR
IN CLOSED GEOMETRIES

Building on the previous analysis of finite cell trains, we
now determine the impact of asymmetric interactions on the
collective dynamics of large cell assemblies. Experimentally,
the only control parameter that is adjustable quantitatively
is the cell density, which we tuned in ring geometries by
varying the number of cells and the size of the system [ring di-
ameter, see Figs. 3(a)–3(d) and S3]. Qualitatively, we observe
that an increase in cell density is associated with an increase
of both cell-clustering and coordinated polarization, as ex-
pected. As density increases, cells form larger clusters that
are more persistent, in agreement with our analysis of single
trains above. In Fig. 3(e), we clearly show the establishment
of long-range coordinated motion in a large ring geometry at
high density of cells.

To make this analysis quantitative, we chose not to use the
distribution of cluster sizes because of the limited statistics
that were accessible experimentally for each value *N* of the
cell number. Instead, to quantify the degree of clustering in the

cell assembly, we defined two phenomenological parameters:
d_l measuring the effective distance between an experimental
configuration and a perfect phase-separated state (i.e., a single
cell aggregate) and d_g measuring the distance between an
experimental configuration and an ideal Poisson distribution
(i.e., a uniform distribution of nonoverlapping cells) (see de-
tails in the SM and Fig. S2 [52]). From these quantities, we
introduce a clustering index d_g - d_l which, by construction,
differentiates between dispersed configurations (d_g - d_l < 0)
and clustered ones (d_g - d_l > 0). In Fig. 3(f), we show that
this clustering index increases with *N* and plateaus beyond
a value *N* significantly lower than *N_c* where confluence is
reached (*N_c* ≈ 65 for a ring diameter *D* = 1000 µm, see also
Fig. S3). This means that due to reduced available space, cells
tend to cluster when the density increases. In turn, as PBD-
YFP fluorescence was not accessible in those experiments, we
identified the polarity of each cell with p_i ≡ v_i/|v_i| (where
v_i is the velocity of cell i) in order to assess the average
polarity ⟨|p|⟩ within a ring. As discussed qualitatively above,
we observe that the steady-state polarity averaged over time
and realizations ⟨|p|⟩ increases with cell number *N*, which
is consistent with the reported increase of cluster size with
<figure><img src="image_5.png" /><figcaption>FIG. 4. Polar-disordered transition and microphase separation in the microscopic model. (a) Examples of steady-state structures observed in simulations of the model for various conditions $\alpha \ll \beta$ and $\alpha \gg \beta$ with $\mu = 100$ or $\mu = 1$. Cluster size statistics as a function of the control parameters of our model: the relative strength $\alpha/\beta$, the self-propulsion magnitude $Pe/\gamma$, and the relaxation rate $\mu$—(b), (e) in the ($\alpha/\beta$, $Pe/\gamma$) plane for $\mu = 100$; (c), (f) in the ($\alpha$, $Pe/\gamma$) plane for $\mu = 2$ and $\beta \ll 1$; (d), (g) in the ($\mu$, $Pe/\gamma$) plane for $\alpha/\beta \gg 1$. In (h), curves are given for logarithmically spaced values of $\mu$ from $10^{-1}$ (green) to $10^2$ (blue). In (i)–(l), we have added as visual guides for the re-entrant clustering transition a single contour at arbitrary values of $\langle s_c \rangle / N$ (solid black line).</figcaption></figure>

$N$, and the existence of typically one domain wall per cluster as reported in Fig. 2; in particular, polarity significantly exceeds the expected value for a random choice of $N$ polarity vectors $p_i$, and plateaus for confluent systems with $N > N_c$ [see Fig. 3(g)]. Hence, we can delineate a route to collective motion in this dense confluent system: as density increases, the available space decreases, which leads to larger, less numerous cell clusters; being an edge term, the CIL interaction becomes less prominent even though it dominates the dynamics of small systems, which allows collective motion to arise.

## VI. HYDRODYNAMIC THEORY AND PHASE DIAGRAM

Experiments are limited to rather small cell numbers, finite time scales and do not allow for the independent tuning of key control parameters. To explore the effect of asymmetric

interactions on the possible phases in the thermodynamic limit, we rely on both the hydrodynamic limit of the agent-based model introduced above and on a numerical analysis of the 1D version of this model [see Fig. 4(a)].

First, we place ourselves in the case of a high polarity relaxation rate $\mu = 100$; in this limit, the polarity of the particles quickly relaxes. Qualitatively, we observe in the limit where alignment dominates over antialignment the stability of macroscopic polar drops independently of the self-propulsion magnitude $Pe/\gamma$ [see Figs. 4(a)–4(e)]; all cells are moving coherently as a single drop. The observed behavior sharply transitions when the strength of the asymmetric term becomes larger than the strength of the symmetric term. For $\alpha \gg \beta$, we observe a nonmonotonic evolution of the average cluster size as a function of self-propulsion magnitude. In the limit of $Pe/\gamma \to 0$, we also observe a stable single drop but the
drop is, in this case, apolar on average. Indeed, in the limit of fast polarity relaxation $\mu \gg 1$, as discussed above the energetically favored configuration contains a single domain wall that separates groups of particles with outward pointing polarities. Over time, the domain wall diffuses in this cohesive drop. As $Pe/\gamma$ starts increasing, the average cluster size decreases sharply. In this regime, the self-propulsion magnitude competes with cohesive forces, leading to what we call cohesion-induced microphase separation. Eventually, the self-propulsion magnitude becomes large enough for the small drops to fully disassemble and a disordered dilute (gaslike) phase is reached.

Surprisingly, we observe phase reentrance; indeed, as we increase the self-propulsion magnitude again, we observe that the typical cluster size in the system increases. We observe small drops with transient single domain walls now separating regions of particles with inward pointing polarities. This second microphase separation is in fact caused by a motility induced phase separation (MIPS) mechanism. It is thus expected that the stability of the disordered dispersed (gas) phase against MIPS clustering is critically controlled by the persistence time of the cell polarities; as the polarity relaxation rate $\mu$ decreases, particles become more persistent and activity promotes microphase separation at lower self-propulsion strength [see Figs. 4(d) and 4(g)].

We now step away from this microscopic picture and construct a general hydrodynamic theory of asymmetrically interacting active polar particles in arbitrary dimensions that encompasses the agent-based model introduced above and is valid for all models with the same symmetry class [52]. Since the particles interact with a substrate, momentum is not a conserved variable; hence, the hydrodynamic fields that we must retain are the particle density $\rho(\mathbf{r}, t)$, which is a conserved quantity, and the polarity $\mathbf{p}$. Note that even though the polarity $\mathbf{p}_i$ of a single particle is a unit vector, the corresponding coarse-grained variable $\mathbf{p}$ is not. The microscopic interaction potential $H$ defined in Eq. (1) can be coarse-grained to yield an effective free-energy in terms of the variables $\rho$ and $\mathbf{p}$ [59]. We do not perform this procedure explicitly here and instead provide a phenomenological expression based on symmetry arguments. The classical alignment term $-\beta\mathbf{p}_i \cdot \mathbf{p}_j$ in the microscopic polarity potential $U_p$ defined in Eq. (3) generically yields the usual Landau free energy [60]

$$F_p = \int d\mathbf{r} \left[ \frac{A}{2} p^2 + \frac{B}{4} p^4 + \frac{K}{2} (\nabla \mathbf{p})^2 \right], \quad (7)$$

where all terms, in particular, the parameter $A$ that controls the isotropic-polar transition, can be functions of the density. Note that the $(\nabla \mathbf{p})^2$ term stands for the classical Frank elasticity in the one constant approximation, which comprises splay and bend contributions [60].

As stated above, the asymmetric alignment term $-\alpha(\mathbf{p}_i - \mathbf{p}_j) \cdot \mathbf{n}_{ij}$ in Eq. (3) breaks the invariance under independent rotations of space and polarity, but preserves the required invariance under joint rotations of space and polarity that characterizes equilibrium polar liquid crystals. Upon coarse graining, any asymmetric alignment term with such symmetry generically yields a coupling between density $\rho$ and space derivatives of $\mathbf{p}$; to lowest order in these hydrodynamic fields,

we therefore write without loss of generality

$$F_{\rho p} = - \int d\mathbf{r} \bar{\alpha} \delta\rho \nabla \cdot \mathbf{p}, \quad (8)$$

where the coefficient $\bar{\alpha}$ is proportional to $\alpha$ and $\delta\rho$ is the deviation of the density from a steady state, spatially homogeneous density $\rho_0$ (a free-energy $\propto \rho_0 \nabla \cdot \mathbf{p}$ with a spatially homogeneous density $\rho_0$ would yield only a boundary term and is omitted). This is the spontaneous splay term, well-known in the context of equilibrium polar liquid crystals [61]; in particular, it is apparent from Eq. (8) that for finite particle clusters of uniform density, $F_{\rho p}$ is merely a boundary term, as we argued on the basis of the agent-based model above. As we show below, this term has important consequences in active systems that have remained unexplored so far.

Finally, the coarse-graining of the position-dependent interaction potential $U_r$ yields

$$F_{\rho} = \int d\mathbf{r} \left[ U(\rho) + \frac{\kappa}{2} (\nabla \rho)^2 \right] \quad (9)$$

as in standard equilibrium theories of fluids, where $U(\rho)$ is the internal energy density and $\kappa$ is the interface energy constant. Since we are interested in model-independent properties, we choose a simple phenomenological form for $U(\rho)$ that accommodates a liquid-gas transition near a steady-state density $\rho_0$ controlled by a phenomenological parameter $A_c$, i.e.,

$$U(\rho) = \frac{A_c}{2} (\delta\rho)^2 + \frac{B_c}{4} (\delta\rho)^4. \quad (10)$$

The effect of the asymmetric contribution $F_{\rho p}$ to the total free energy $F = F_p + F_{\rho p} + F_{\rho}$ at equilibrium is summarized in the SM for completeness, and shown to be unimportant for $d=1$, and in the homogeneous disordered phase for $d \ge 2$ case that we consider below [52]. We now construct the phenomenological active dynamics of $\rho$ and $\mathbf{p}$. The activity in the model enters primarily as a coarse-grained self-propulsion velocity of the polar particles that we assume enslaved to the polarity $\mathbf{v} = v_0 \mathbf{p}$. The coupled equations for the conserved particle density and the polarity fields are

$$\partial_t \rho = - \nabla \cdot (v_0 \rho \mathbf{p}) + D \nabla^2 \frac{\delta F}{\delta \rho} + \Lambda_1 \nabla \cdot \frac{\delta F}{\delta \mathbf{p}}, \quad (11)$$

$$\partial_t \mathbf{p} = -\bar{\mu} \frac{\delta F}{\delta \mathbf{p}} - \Lambda_2 \nabla \frac{\delta F}{\delta \rho}. \quad (12)$$

Here $D$ is an effective diffusivity and $\bar{\mu}$ is the polarity relaxation rate which is the coarse-grained counter part of the parameter $\mu$ in the agent-based model. We have included phenomenological couplings between $\mathbf{p}$ and $\rho$ with coefficients $\Lambda_1$ and $\Lambda_2$. The detailed-balance obeying limit is defined by $v_0 = 0$ and $\Lambda_1 = \Lambda_2$, imposed by Onsager symmetry, and corresponds to the linear response of a passive system of interacting polar particles. However, out of equilibrium, there is no symmetry to enforce this equality. For instance, these coefficients could be distinct due to microscopic nonreciprocal interactions between polarities [62]. Note that in principle, Eqs. (11) and (12) also contain nonlinear terms that cannot be derived from a potential (such as advective and self-advective
<figure><img src="image_6.png" /><figcaption>FIG. 5. Phase diagram represented, respectively, in (a) the (ā, v₀) plane with µ̄ ≫ 1 (sketch), (b) the (ā, v₀) plane with µ̄ ~ 1, and (c) the (µ̄, v₀) plane with ā ≫ 1. The black solid line in (b) and (c) is the spinodal line associated to Eq. (13), which defines the domain of linear stability of the homogeneous disordered phase (red region). Here, we assumed that Λ₂(v₀) is a decreasing quadratic function of v₀ with Λ₂(v₀ = 0) = 0. All three phase diagrams clearly display rerentrance to a clustered phase as the self-propelling velocity v₀ is increased.</figcaption></figure>

nonlinear terms). These, however, do not affect the linear analysis below and are omitted for simplicity. Note, finally, that while the microscopic model has no nonreciprocal interaction, unlike the model introduced in Ref. [58], the coarse-grained equations we obtain are equivalent. In other words, in this system, microscopic nonreciprocity does not yield any effect not already accounted for by the combination of generic asymmetric interaction and motility that we introduce in our model.

We now linearize Eqs. (11) and (12) about a disordered but homogeneous phase with density ρ₀ (hence A > 0). Upon elimination of the polarity field which relaxes fast, one obtains [52]

$$ \partial_t \rho = \left[ v_0 \rho_0 \left( \frac{\bar{\alpha}}{A} + \frac{\Lambda_2 A_c}{A \bar{\mu}} \right) - \frac{\Lambda_1 \Lambda_2 A_c}{\bar{\mu}} + A_c D \right] \nabla^2 \rho. \quad (13) $$

In equilibrium, i.e., for v₀ = 0 and Λ₁ = Λ₂, the diffusivity is simply renormalized to D - Λ₁²/μ̄ > 0 and the homogeneous disordered state is linearly stable, as expected, for A_c > 0; the usual spinodal line of the liquid-gas transition is then simply given by A_c = 0.

At this stage, it is useful to compare this analysis with the equation of motion for the density and polarity fields in Ref. [63]. These latter equations can be recovered in our formalism by taking Λ₁ = ᾱ = 0 and Λ₂ ≠ 0. The condition A_c(D + v₀ρ₀Λ₂/Aμ̄) < 0 is then equivalent to the MIPS [47,48] spinodal obtained in Ref. [63], which arises in our notation when Λ₂A_c < 0, and which is indeed expected to be present in our system, which has both self-propulsion and hard-core repulsion. To be consistent with Ref. [63], we therefore assume from now on that Λ₁ = 0, and that A_cΛ₂(v₀) is a decreasing function of v₀ with Λ₂(v₀ = 0) = 0.

With these prescriptions, the spinodal line is explicitly determined by Eq. (13) and defines the domain of linear stability of the homogeneous disordered phase, which we plot in the (ā, v₀) and (μ̄, v₀) planes in Fig. 5. This analysis, however, does not allow us to fully describe the expected phases beyond the spinodal lines. Therefore, we make use of the results of numerical simulations of the 1D version of the agent-based model to characterize the inhomogeneous phases (see Fig. 4 and Appendix F for details). Confronting Figs. 4 and 5, we observe that our numerical results qualitatively confirm the analytical predictions. Both analytical and numerical analysis show that the asymmetric term ᾱ, when coupled to activity (v₀ ≠ 0) stabilizes the disordered homogeneous phase, both

with respect to the usual, passive, mean-field transition to liquid induce by cohesive cellular forces (the homogeneous gas phase can be stabilized even for A_c < 0 when ᾱ is increased) and with respect to MIPS (whose spinodal can be suppressed by increasing ᾱ).

Strikingly, starting at equilibrium (v₀ = 0) deep in the liquid phase (A_c < 0) induced by cohesive interactions (CI), we find that for ᾱ > 0, increasing self-propulsion first destabilizes clusters and induces a microphase of finite self-propelled clusters, which can ultimately lead to a fully dispersed phase. This is in agreement with the mechanism of cluster fragmentation analyzed in Fig. 2(b), and proves its impact at larger scales. For this mechanism to occur, the asymmetric interactions have to dominate over the aligning ones (α > β) for β > 1 or over the thermal fluctuations (α > 1) for β < 1 (where α and β are expressed in units of T_p). In addition to this fragmentation mechanism, we show in Fig. 3 that activity increases self-propulsion and persistence of clusters, which thus behave as mesoscale SPPs with hard-core repulsion. Upon increasing activity, a competing MIPS mechanism of either single particles or clusters favoring aggregation is therefore expected. Such reentrance into a clustered phase is indeed predicted by our stability analysis and numerical simulations (Figs. 4 and 5). This MIPS-induced (MI) clustering mechanism is expected to be critically controlled by the persistence time of the clusters, which is in turn controlled by the timescale μ⁻¹ of polarity dynamics, as confirmed by the linear stability analysis—see the term ρ₀v₀A_cΛ₂/(Aμ̄) in Eq. (13). Indeed, we find that for large values of μ, the MIPS phase is pushed to larger values of v₀, thus increasing the stability domain of the disordered dispersed phase; accordingly, for low values of μ, the MIPS phase (MI) appears even at low v₀, thereby completely masking the fragmentation mechanism and preventing the emergence of a microphase (Figs. 4 and 5).

## VII. EXTENSION TO TWO-DIMENSIONAL SYSTEMS

Up to now, we have used simple 1D geometries to make a detailed analysis of the system tractable, both experimentally and computationally; in our integrative bottom-up approach, we have done this systematically from the scale of pairs of cells up to large confluent systems. Yet, due to the generic nature of the d-dimensional hydrodynamic theory presented in the previous section, we expect that several of our key
<figure><img src="image_7.png" /><figcaption>FIG. 6. Two-dimensional cell clusters.—We analyze examples of experimental and numerical thick ring geometries and finite clusters. (a) Example of a steady-state configuration in a thick ring geometry (where the thickness of the ring is $R = 246$ µm); red arrows show the local velocity in the ring as measured by PIV (scale bar = 100 µm). (b) Steady-state polar order parameter $\langle|p|\rangle$ averaged over realizations and time as a function of the number of cells in thick ring geometries with ring thicknesses, $R = 26$ µm (blue circles, 82 patterns per experiment), 49 µm (red squares, 76 patterns per experiment), 124 µm (orange diamonds, 85 patterns per experiment), 246 µm (green lower triangles, 97 patterns per experiment), and 395 µm (grey upper triangles, 92 patterns per experiment). For each geometry, data was binned by the number of cells. Vertical error bars are given by standard deviations over the polarization, while horizontal error bars are standard deviations over the number of cells. Data obtained from the analysis of three independent experiments. (c) Example realizations of small two-dimensional cohesive cell clusters. Red arrows represent the cell polarity and are estimated by the nucleus-Golgi axis. (d) Polarity of the cells in experimental two-dimensional cohesive clusters relative to the edge of the cluster as a function of the distance from the edge. We show here the median $\mu(p \cdot \hat{n})$ as blue symbols, where positive values of $p \cdot \hat{n}$ denote cells pointing towards the outward normal vector to the edges of the clusters. Error bars are given by the first and third quartiles. (e) Example realizations of two-dimensional cohesive clusters in the model with $\alpha > \beta$ showing polarities at the edge of the cluster pointing along the local outward normal direction. Given boundary conditions, the polarities of the cells inside the clusters are organized around a single +1-charge nematic defect (aster).</figcaption></figure>

predictions will still be valid in higher physically relevant dimensions, in particular, for $d = 2, 3$. To prove this point, we focus on demonstrating in two-dimensional (2D) cell monolayers two of our main predictions characterizing CIL as an edge term: First, in finite-size clusters, we expect cells close to the edge to point preferentially along the normal vector to the edges of the cluster and, second, in a large system with periodic boundary conditions (i.e., in closed-ring geometries), large-scale coordinated flow (global polarity) will emerge as the cell surface coverage approaches full confluency, because in this regime cell cluster edges are absent and the effect of CIL thus vanishes.

First, we extend to two spatial dimensions the analysis of the dynamics of cells in our ring geometries. To do so, we study annuli whose thicknesses are large compared to a single cell typical size; an example of such a geometry is given in Fig. 6(a). This experimental setup allows us to explore cell dynamics in two dimensions while imposing periodic boundary conditions, albeit only in the orthoradial direction whereas the other boundaries impose a tangential flow con-

dition. Note that imposing periodic boundary conditions in all spatial directions would require us to study the motion of cells on the surface of a three-dimensional torus which is outside the realm of this study. In our experiments, we use annuli with an outer diameter of 1000 µm and five different thicknesses ranging from 26 µm to ~395 µm. We observe that at low density, cells exhibit 2D random walks, while in confluent rings, a global polar flow emerges in the orthoradial direction. To quantify the polar order in these systems, we proceed as done previously and measure the instantaneous ring polarity $p(t)$ following Eq. (B4). In Fig. 6(b), we show that, in the steady state, $\langle|p|\rangle$ increases as the number of cells in the annuli increases (i.e., with increasing cell volume fraction), and that for all ring thicknesses. This confirms the emergence of large-scale coordinated motion, i.e., the emergence of coordinated polarity when the system reaches confluency (see SM Movie S5).

Second, we generate long-lived, cohesive two-dimensional (2D) clusters. Starting from sparse single cells, we let them grow until small colonies of up to few hundred cells are
generated. For $d = 2$, we found the PBD-YFP signal to be
quite unreliable at accurately predicting the cell polarities.
Instead, we fix cells and fluorescently label the nucleus and
Golgi apparatus [see Fig. 6(c)], as the nucleus-Golgi axis has
been previously shown to quite reliably reflect cell polarity
[64–66]. With this definition in mind, we show in Fig. 6(d)
that the outermost cells significantly polarize outwards at the
edge of the cluster by quantifying the value of $\mathbf{p} \cdot \hat{\mathbf{n}}$, where $\mathbf{p}$
denotes the cell polarization vector and $\hat{\mathbf{n}}$ the vector joining
the center of mass of the cell to the closest point on the outer
edge of the cell cluster; positive values of $\mathbf{p} \cdot \hat{\mathbf{n}}$ denote cells
pointing towards the outward normal vector to the edges of
the clusters. We show in Fig. 6(e) sample realizations of such
long-lived cohesive clusters in our microscopic numerical
model; we confirm that as long as $\alpha > \beta$, the clusters will
show cell polarities at the edge of the cluster pointing along
the local outward normal direction. Because of this boundary
condition imposed by the CIL asymmetric interaction, the
polarities of the cells inside the clusters self-organize around a
single +1-charge nematic defect (aster), which is the expected
ground-state solution of the polarity dynamics when $\alpha > \beta$
in 2D. The dynamics of this topological defect thus controls
the dynamics of the overall cluster. Of note, for $d = 2$, the
motion of this topological defect inside the cohesive cluster is
not free diffusion because of elastic interactions and thus we
expect the dynamics to be different from that of clusters for
$d = 1$.

Overall, we conclude that key features of the model, observed for $d = 1$ and predicted by our hydrodynamic theory, are conserved for 2D cell assemblies: CIL interactions act as boundary terms only and lead to edge cells pointing outwards in finite clusters, and to the emergence of coordinated motion at confluency where cluster edges disappear and CIL effects become negligible.

VIII. CONCLUSION

Our results illustrate how cell-cell interactions regulate the
collective behavior of cellular systems and their organization.
Based on a joint experimental and theoretical approach, we
analyzed the impact of generic asymmetric interactions remi-
niscent of the CIL interactions reported for various cell types
on the collective dynamics of cell assemblies, and more gener-
ally of dry active systems. We made use of microfabricated 1D
*in vitro* environments to characterize quantitatively pairwise
cell-cell interactions and showed that the observed CIL-type
phenomenology can be captured by a generic equilibriumlike
asymmetric aligning interaction potential that breaks the usual
invariance under *independent* rotations of space and polarities.
Based on experimental observations, numerical simulations
and analysis of the relevant active hydrodynamic theory, we
demonstrated that such an asymmetric aligning interaction
can drastically lower the size of cell clusters and control their
self-propulsion speed and persistence. By carefully inspecting
the role of cell density in closed geometry, we added a piece
to the puzzle of persistent coherent motion of cells confined to
a ring: while CIL-like interactions should tend to destroy any
large-scale order, their effect is mitigated by the disappearance
of cluster edges. We have also confirmed our predictions by

studying cell polarity in simple 2D systems, including thick
ring geometries and isolated cellular clusters.

In the large system limit, we found that this can lead
to the emergence of a liquidlike microphase of cell clusters
of finite size and short-lived polarity, and ultimately stabi-
lish a fully dispersed apolar phase. Altogether, this analysis
suggests that CIL-like asymmetric interactions in generic ac-
tive systems—cellular or artificial—can control cluster sizes
and polarity, and can thus prevent large-scale coarsening and
long-ranged polarity, except in the singular regime of dense
confluent systems. While our experimental and numerical
analysis was focused on 1D and 2D geometries, we expect
that several key features, predicted by the generic analysis
of the *d*-dimensional hydrodynamic theory, are still valid for
$d = 3$, and therefore relevant to *in vivo* biological systems. We
anticipate that the mechanism of active cluster fragmentation
induced by the CIL interaction, and the critical dependence of
size, speed, and persistence of cell clusters on the CIL inter-
action may provide a unique mechanism to interpret directed
cell migration during development, epithelial-mesenchymal
transition, and collective cancer cell invasion.

ACKNOWLEDGMENTS

We acknowledge the ImagoSeine core facility of the
Institut Jacques Monod and members of IBiSA and France-
BioImaging (ANR-10-INBS-04) infrastructures. The authors
are grateful to F. Martin-Belmonte for their generous gift of
MDCK cell lines. We acknowledge S. Grigolon for useful
discussions. We also acknowledge financial support from the
Agence Nationale de la Recherche (ANR) POLCAM (ANR-
17-CE13-0013) (R.M.M., B.L., R.V.), the LABEX “Who
am I?” (ANR-11-LABX-0071), Université de Paris IdEx
ANR-18-IDEX-0001, the Ligue Contre le Cancer (Equipe la-
bellisée) (B.L., R.M.M.), Défi Mécanobio CNRS MECAPOL
2016 (R.M.M., R.V.) and Institut National du Cancer
(INCA_16712, B.L.) and CNRS Program 80 Prime (R.V.,
B.L.). R.V. also acknowledges support from Fondation pour
la Recherche Médicale (FRM) and Institut National du Cancer
(INCA).

T.B. and J.d'A. contributed equally to this work. T.B.,
J.d'A., B.L., and R.V. designed the research. B.L. and R.V.
supervised the project. T.B., A.M., and R.V. developed the
theoretical model and numerical simulations. J.d'A., B.M.,
and S.J. performed the experiments. T.B., J.d'A., R.M.M.,
B.L., and R.V. analyzed the data.

APPENDIX A: EXPERIMENTAL METHODS

1. Cell culture

We used MDCK wild-type, MDCK histon-GFP, and
MDCK PBD-YFP (a gift from F. Martin-Belmonte lab-
oratory). The cells were cultured in DMEM GlutaMAX
high-glucose (Gibco, Waltham, MA) supplemented with 10%
fœtal bovine serum (BioWest, Nuaillé, France). Prior to ex-
periments, the cells were treated with mitomycin C at final
concentration of 10 µg.mL added in the medium for 1 h,
then rinsed before subsequent detachment and seeding on the
experimental samples.
## 2. Sample preparation

All micropatterns were prepared using standard microcontact printing on PDMS, as described in Ref. [44]. The substrates used were (i) nonculture treated plastic dishes (Greiner Bio-One, Kremsmünster, Austria) for wide-field microscopy or (ii) glass coverslips (Menzel-Gläser) for spinning-disk microscopy. The substrates were first covered with a thin layer of poly-dimethyl-syloxane (PDMS, Sylgard, Dow Corning, Midland, MI) using a spin-coater and crosslinked at 80 °C for 2 h. PDMS stamps were made by pouring PDMS on a mold featuring the patterns to be printed and crosslinked as described. After cooling down, a fibronectin solution was prepared by adding 5 mg.mL⁻¹ of fibronectin and 2.5 mg.mL⁻¹ of Cy3- or Cy5-labelled fibronectin into sterile milliQ water. The solution was then incubated on the stamps for 40 min at room temperature. Before stamping, the substrates were activated using UV-ozone for 10 min; the stamps were rinsed to remove any excess fibronectin and dried using an air-gun. The stamps were briefly put in contact with the surface of the substrate, then removed, and the substrates immerged in a 2% pluronics F127 (Sigma-Aldrich, Saint Louis, MO) solution in PBS for 2 h. Finally, the substrates were rinsed in PBS and sterilized under the UV lamp of a culture hood before use.

## 3. Stencils

To be able to track cell trains of definite length over a long period of time, we prepared isolated trains as follows. We fabricated PDMS microstencils by cutting a trapezoidal shape through a thin (approximately 100 µm) layer of PDMS using a cutting plotter (Graphtec CE6000-40, Graphtec Corp., Yokohama, Japan). The stencil was then placed on top of the linear patterns at a 90° angle to leave from a few dozens of microns up to 400 µm of uncovered space at the middle of each line. The stencil was removed after the standard cell seeding, attachment, and rinsing steps, so as to let the trains move freely without future encounters.

## 4. Cell seeding

The cells were enzymatically detached, then concentrated using a centrifuge and seeded on the substrates in culture medium, at a controlled density: medium-low density for the doublet and random small train experiments, very high density for the stencil experiments, a wide range of densities for the ring experiments. The cells were let to adhere in the incubator for approximately 45 min, then rinsed thoroughly (but carefully) to remove excess floating cells without affecting adhered cells.

## 5. Time-lapse microscopy

All experiments were run at 37 °C in 5% CO₂. The experiments in Figs. 1 and 2(b) were done using an inverted microscope (Leica, Wetzlar, Germany) with a CSU-W1 confocal spinning-disk module (Nikon, Tokyo, Japan) and a 40X oil-immersion objective. The acquisition was done using Metamorph (Molecular Devices, San Jose, CA), at a 6- to 10-min. acquisition rate. The focus was done on the basal

plane of the cells and the microscope's hardware autofocus was used to ensure the absence of defocusing.

All the other experiments were performed with a wide-field inverted microscope (Olympus, Tokyo, Japan) using a 10X air objective. Phase contrast and green fluorescent protein (GFP) fluorescence images were acquired using Metamorph (Molecular Devices, San Jose, CA), at a 6- to 12-min. acquisition rate. An image of the labeled patterns was done at least at the beginning of the experiment to allow further alignment.

## 6. Cell colony preparation and immunostaining

To prepare 2D cell colonies, cells were seeded at very low density on glass coverslips and let to divide for four days to yield on average 100 cells per colony. At this point, the cells were fixed with 4% paraformaldehyde for 20 min. The cells were then immuno-stained using classic procedure: after rinsing in PBS, the cells were treated 10 min in 0.5% Triton X100, then rinsed and incubated in 10% FBS, 1% BSA for 1 h. After that, they were incubated with primary antibodies overnight with 1% FBS and 1% BSA at 4 °C. The next day, primary antibodies were rinsed; they were further incubated with secondary antibodies, together with Hoechst 33342 and Phalloidin for 2 h at room temperature. Finally, the samples were rinsed and mounted in Prolong before imaging. The antibodies used were anti-GM130 to mark the Golgi. The colonies were then imaged on a wide-field fluorescence microscope at 20x magnification with a 2 x 2 binning of the pictures.

# APPENDIX B: IMAGE ANALYSIS

## 1. Analysis of isolated trains

The line patterns were first detected using an in-house macro allowing either automatic thresholding or semiautomatic drawing. After rotation and stitching of the images using in-house macros, the phase contrast images were binarized using a simple thresholding approach with adapted smoothing, dilatation-erosion steps and filters on both the size and circularity of the detected objects. This way, the pixels occupied by cohesive trains could be detected. The GFP pictures were then treated using Imaris (Belfast, UK) to get the trajectories of single nuclei, with a manual correction step allowing to get experiment-long clean trajectories. Further analysis used a combination of those two data sets: tracked position of the single nuclei and position and extension of cohesive trains.

## 2. Analysis of ring geometries

The analysis of ring experiments followed the same spirit but with slightly different procedures. The detection of the rings was done using in-house macros performing template matching; the images were then cropped to get single rings of all given diameters in different movies. The train detection was done as previously, but the nuclei detection relied on an in-house macro based on the Find Maxima function of ImageJ [38], and subsequent tracking was done using the track.m function in MATLAB [67]. This procedure induced—low rate—random errors in the trajectories but allowed much
larger outputs in terms of number of cells analyzed. All quantities, in particular, the velocities, were projected on the orthoradial direction to get a 1D data set. In parallel, we performed particle-image velocimetry (PIV) using MatPIV [68] on the phase-contrast images. The velocity field, we obtained, was similarly projected onto the orthoradial direction, and the signal was filtered using the trains location data to remove spurious velocities outside of the areas covered by cells [see Fig. 3(e)]. We thus obtained the cell velocities from two independent protocols: a discrete Langrangian method and a continuous Eulerian method.

### 3. Velocity autocorrelation functions

In both isolated trains and rings, the velocity autocorrelation functions were computed in direct space for each cell i as follows:

$$C_i(\Delta t) = \langle \mathbf{v}_i(t) \cdot \mathbf{v}_i(t + \Delta t) \rangle_i. \quad (B1)$$

The autocorrelation function $C_i$ was then normalized to obtain $c_i(\Delta t) = C_i(\Delta t)/C_i(0)$ and for each cell number $N$, the average $c_{avg}(\Delta t)$ was taken over all cells belonging to a train of (initially) $N$ cells. The persistence time $\tau_v$ was defined as the first $\Delta t$ this average function decayed below $e^{-1}$. The error in $\tau_v$ was obtained by applying the same threshold detection to the functions $c_{avg} \pm c_{sem}$ where $c_{avg}$ and $c_{sem}$ are, respectively, the average and standard error of the mean of $c_i$ for a given $N$.

### 4. Polar order parameter

The polar order parameter was computed on single rings with two different definitions based on the two available velocity data sets. When using individual cell tracking data, we computed the polar order parameter as follows:

$$\mathbf{p}_{\text{tracking}}(t) = \frac{\sum_{i=1}^{N} \mathbf{v}_i(t)}{\sum_{i=1}^{N} |\mathbf{v}_i(t)|}, \quad (B2)$$

where $\mathbf{v}_i(t)$ denotes the velocity of cell i at time $t$ in a ring of $N$ cells. When using the velocity field obtained via PIV, we defined the polar order parameter as

$$\mathbf{p}_{\text{PIV}}(t) = \frac{\sum_{i \in \mathbb{T}(t)} \tilde{\mathbf{v}}_i(t)}{\sum_{i \in \mathbb{T}(t)} |\tilde{\mathbf{v}}_i(t)|}, \quad (B3)$$

where $\tilde{\mathbf{v}}_i(t)$ denotes the velocity measured on a coarse-grained “pixel” i at time $t$ and $\mathbb{T}(t)$ is the ensemble of pixels covered by cells at time $t$. We checked that those two definitions were consistent. In what follows, we mostly used $\langle |\mathbf{p}| \rangle_{\text{PIV}}$— while checking again that the subsequent results were not affected by the definition of $\langle |\mathbf{p}| \rangle$, where $\langle \cdots \rangle$ denotes an average over configurations. To plot the steady-state values of $\langle |\mathbf{p}| \rangle$ against $N$, we waited for $t \simeq 30$ h for the system to stabilize and reach stationarity. We then averaged the $\langle |\mathbf{p}| \rangle$ value for all $t \ge 30$ h in individual rings, then binned the data obtained over number of cells $N$ (using six bins). We then computed both average and standard deviation over these bins, as shown in Fig. 3(g) for rings of diameter $D = 400$ µm.

For thick annuli experiments, we used a similar definition for the polar order parameter, namely,

$$p(t) = \frac{\sum_{i \in \mathbb{T}(t)} \tilde{v}_i^{\theta}(t)}{\sum_{i \in \mathbb{T}(t)} |\tilde{v}_i(t)|}, \quad (B4)$$

where $\tilde{v}_i(t)$ and $\tilde{v}_i^{\theta}(t)$ are, respectively, the velocity vector and its orthoradial component measured on coarse-grained pixel i at time $t$ and $\mathbb{T}(t)$ is the ensemble of pixels covered by the cells at time $t$.

### 5. Cell polarity in 2D clusters

The nucleus-Golgi axis was used to estimate the direction of cell polarity in immunolabelled 2D clusters. First, the positions of the nuclei and the Golgi bodies were detected using an in-house ImageJ macro based on the Find Maxima function that locates local maxima in the image. Then, nuclei and Golgi bodies were paired using a squared-distance minimization algorithm [67] and the nucleus-Golgi vector was used to determine the orientation of cell polarity. The colonies were segmented based on the phalloidin (actin) signal, so as to define cell colony edge and center. Then, the polarization parameter is simply $\mathbf{p} \cdot \hat{\mathbf{n}} = \cos(\theta - \phi)$, where $\phi$ is the orientation of cell polarity $\mathbf{p}$ while $\theta$ is the orientation cell-to-nearest-colony-edge vector $\hat{\mathbf{n}}$. The midpoint between nucleus and Golgi was used to define the cell position to remove any potential bias in defining $\theta$. Hence, $\mathbf{p} \cdot \hat{\mathbf{n}} = 1$ if the cell points to the colony edge (or away from the colony center), while $\mathbf{p} \cdot \hat{\mathbf{n}} = -1$ if it points to the colony center (or away from the edge).

## APPENDIX C: RELATIVE STRENGTH OF THE ALIGNMENT INTERACTIONS

The evolution of the polarities is governed by an equilibrium process. From the measure of the probabilities for each polarity doublet configuration [see Fig. 1(e)], we can estimate the value of the relative strength of the alignment interactions as follows. We denote the four possible configurations of polarity doublets: $A(\to \to)$, $B(\leftarrow \leftarrow)$, $C(\leftarrow \to)$, and $D(\to \leftarrow)$. The respective probabilities of these configurations are given by

$$p_A = \frac{e^\beta}{Z}, \quad p_B = \frac{e^\beta}{Z}, \quad p_C = \frac{e^{2\alpha-\beta}}{Z}, \quad p_D = \frac{e^{-2\alpha-\beta}}{Z}, \quad (C1)$$

where $\alpha$ and $\beta$ are expressed in units of the polarity temperature and $Z$ is the partition function. As we do not have access to the partition function, we take ratios of these probabilities to obtain the following final expressions:

$$\alpha = \frac{1}{4} \log(p_C/p_D), \quad (C2)$$

$$\beta = \alpha + \frac{1}{2} \log(p_A/p_C). \quad (C3)$$

We obtain $\alpha \approx 0.44$ and $\beta \approx 0.22$, i.e., $\alpha/\beta \approx 2$.

## APPENDIX D: POLAR ORDER PARAMETER FOR RANDOM CONFIGURATIONS

To compute the value of the polar order parameter expected for random configurations [see Fig. 3(g)], we consider a set
of N cells of which $n_+$ have a positive polarity and $n_-$ have a negative polarity. The global polarization of the system is given by

$$p(n_+, n_-) = \frac{n_+ - n_-}{N}. \qquad (D1)$$

As $N = n_+ + n_-$, we can characterize the polarization using solely the number of negative polarities, which we denote $n$ and write

$$p_n = \frac{N - 2n}{N}. \qquad (D2)$$

If we consider that each cell polarity is assigned randomly, we can write that

$$\langle |p_N| \rangle = \sum_{n=0}^{N} a_n |p_n|, \qquad (D3)$$

where $a_n$ is the probability to randomly pick $n$ cells (those with negative polarities) out of the $N$ cells given by the binomial coefficient:

$$a_n = \frac{1}{2^N} \binom{N}{n}. \qquad (D4)$$

This leads to

$$\langle |p_N| \rangle = \frac{2}{2^N N} \binom{N}{1 + \lfloor N/2 \rfloor} \left[ 1 + \lfloor \frac{N}{2} \rfloor \right]. \qquad (D5)$$

## APPENDIX E: NONDIMENSIONALIZATION OF THE MICROSCOPIC MODEL

In 1D, we nondimensionalize the Langevin equations describing our microscopic model [see Eqs. (4) and (5)] using $\sigma$ (the particle size) and $\epsilon$ (the Lennard-Jones energy scale) as basic units of length and energy and as unit of time $\tau = \sigma^2/D$, where $D = T/\zeta$ is the self-diffusion coefficient of the ABP. To do so, we introduce the self-propulsion velocity $v_0 = DF_p/T = F_p/\zeta$. We can thus define the nondimensional Péclet number $Pe = v_0\tau/\sigma = v_0\sigma/D$ which measures the ratio between the strength of the self-propulsion and thermal fluctuations,

$$\dot{r}_i = \gamma F_i + Pe p_i + \sqrt{2\eta_i}, \qquad (E1)$$

where we defined $\gamma = \epsilon/T$ as the ratio of the strength of the Lennard-Jones potential to the thermal fluctuations and $F_i$ as the total nondimensionalized truncated Lennard-Jones force on particle $i$.

We consider the polarities to be in contact with a heat reservoir temperature $T_p$; the dynamics of the polarities is thus governed by flips between the values $p = +1$ and $p = -1$ with a given rate per unit of time $\mu$ [69]. Here, we use the single spin-flip kinetic Ising model (also known as Glauber dynamics). This model is defined in terms of a Markovian master equation for the probability distribution

$P(p_1, \dots, p_N, t)$, where $p_i \in \{-1, +1\}$,

$$\begin{aligned} \frac{d}{dt} P(p_1, \dots, p_N, t) \\ = & - \sum_i w(p_i \rightarrow -p_i) P(p_1, \dots, p_i, \dots, p_N, t) \\ & + \sum_i w(-p_i \rightarrow p_i) P(p_1, \dots, -p_i, \dots, p_N, t), \end{aligned} \qquad (E2)$$

where the transition rates $w$ are proportional to $\mu$. The transitions rates satisfy the detailed balance condition

$$\frac{P_0(p_1, \dots, -p_i, \dots, p_N)}{P_0(p_1, \dots, p_i, \dots, p_N)} = \frac{w(p_i \rightarrow -p_i)}{w(-p_i \rightarrow p_i)}, \qquad (E3)$$

where the equilibrium distribution $P_0(p_1, \dots, p_N) = (1/Z) \exp(-U_p/T_p)$ with $T_p$ the polarity temperature, $Z$ is the partition function, and $U_p$ the polarity-polarity interactions which are governed by the following Hamiltonian:

$$U_p = \begin{cases} -\beta p_i \cdot p_j - \alpha(p_i - p_j) \cdot n_{ij} & r_{ij} \le r_c \\ 0 & r_{ij} > r_c, \end{cases} \qquad (E4)$$

where $n_{ij} = r_{ij}/|r_{ij}|$. We can nondimensionalize this Hamiltonian by expressing the symmetric alignment $\beta$ and asymmetric alignment $\alpha$ interaction strengths in units of the polarity temperature $T_p$ (where, in general, $T_p \neq T$). To nondimensionalize the equations governing the polarity dynamics in a way that is consistent with the basic units we detailed above, we express the polarity flipping rate $\mu$ in units of the diffusion timescale $\tau$. The parameter $\mu$ is thus a nondimensional number representing the ratio of the spatial diffusion timescale to the average polarity lifetime. In the limit where $\mu \gg 1$, particles will attempt to flip their polarities a large number of times in the time they require to diffuse by a distance corresponding to their size. Conversely, the limit where $\mu \ll 1$ corresponds to the limit of very persistent polarities.

## APPENDIX F: NUMERICAL METHODS

In 1D, our simulations employed simultaneously: (i) the stochastic Runge-Kutta method [70] to solve Eq. (E1) and (ii) the Glauber algorithm to solve the dynamics of the polarities [69]. In all simulations, we use periodic boundary conditions and initial conditions are chosen so the particles form a single drop with randomly chosen polarities, $p_i = \pm 1$. The maximum time step used is $\delta t = 0.5 \times 10^{-5}$ and we run simulations for 200 $\tau$.

In our implementation of spin-flip dynamics, each particle is endowed with an internal clock. The lifetime of a given polarity is drawn from an exponential distribution with parameter $\mu$. At each time step of the simulation, we start by updating the polarities before moving the particles. To do so, we decrease the clocks of all particles by $\delta t$, list (in order) the particles whose clocks have timed out, and update them as follows: we first calculate the energy difference $\Delta E$ (in units of the temperature) resulting from the flipping of the polarity of the particle and the transition probability $P_t = 1/(1 + e^{\Delta E})$; we accept the flip if $X < P_t$ where $X$ is a uniformly distributed random number between 0 and 1. We then draw from the exponential distribution a new lifetime for the polarity of the particle.
[1] P. Friedl, Y. Hegerfeldt, and M. Tusch, Collective cell migration in morphogenesis and cancer, *Int. J. Dev. Biol.* **48**, 441 (2004).

[2] P. Friedl and D. Gilmour, Collective cell migration in morphogenesis, regeneration and cancer, *Nat. Rev. Mol. Cell Biol.* **10**, 445 (2009).

[3] R. Mayor and S. Etienne-Manneville, The front and rear of collective cell migration, *Nat. Rev. Mol. Cell Biol.* **17**, 97 (2016).

[4] V. Hakim and P. Silberzan, Collective cell migration: A physics perspective, *Rep. Prog. Phys.* **80**, 076601 (2017).

[5] R. McLennan, L. Dyson, K. W. Prather, J. A. Morrison, R. E. Baker, P. K. Maini, and P. M. Kulesa, Multiscale mechanisms of cell migration during development: theory and experiment, *Development* **139**, 2935 (2012).

[6] A. Shellard, A. Szabó, X. Trepat, and R. Mayor, Supracellular contraction at the rear of neural crest cell groups drives collective chemotaxis, *Science* **362**, 339 (2018).

[7] A. G. Clark and D. M. Vignjevic, Modes of cancer cell invasion and the role of the microenvironment, *Curr. Opin. Cell Biol.* **36**, 13 (2015).

[8] A. Palamidessi, C. Malinverno, E. Frittoli, S. Corallino, E. Barbieri, S. Sigismund, G. V. Beznoussenko, E. Martini, M. Garre, I. Ferrara, C. Tripodo, F. Ascione, E. A. Cavalcanti-Adam, Q. Li, P. P. Di Fiore, D. Parazzoli, F. Giavazzi, R. Cerbino, and G. Scita, Unjamming overcomes kinetic and proliferation arrest in terminally differentiated cells and promotes collective motility of carcinoma, *Nat. Mater.* **18**, 1252 (2019).

[9] D. Cai, W. Dai, M. Prasad, J. Luo, N. S. Gov, and D. J. Montell, Modeling and analysis of collective cell migration in an *in vivo* three-dimensional environment, *Proc. Natl. Acad. Sci. USA* **113**, E2134 (2016).

[10] T. Omelchenko, J. M. Vasiliev, I. M. Gelfand, H. H. Feder, and E. M. Bonder, Rho-dependent formation of epithelial "leader"cells during wound healing, *Proc. Natl. Acad. Sci. USA* **100**, 10788 (2003).

[11] R. Farooqui and G. Fenteany, Multiple rows of cells behind an epithelial wound edge extend cryptic lamellipodia to collectively drive cell-sheet movement, *J. Cell Sci.* **118**, 51 (2005).

[12] M. Poujade, E. Grasland-Mongrain, A. Hertzog, J. Jouanneau, P. Chavrier, B. Ladoux, A. Buguin, and P. Silberzan, Collective migration of an epithelial monolayer in response to a model wound, *Proc. Natl. Acad. Sci. USA* **104**, 15988 (2007).

[13] X. Trepat, M. R. Wasserman, T. E. Angelini, E. Millet, D. A. Weitz, J. P. Butler, and J. J. Fredberg, Physical forces during collective cell migration, *Nat. Phys.* **5**, 426 (2009).

[14] S. R. K. Vedula, M. C. Leong, T. L. Lai, P. Hersen, A. J. Kabla, C. T. Lim, and B. Ladoux, Emerging modes of collective cell migration induced by geometrical constraints, *Proc. Natl. Acad. Sci. USA* **109**, 12974 (2012).

[15] B. Ladoux and R.-M. Mège, Mechanobiology of collective cell behaviors, *Nat. Rev. Mol. Cell Biol.* **18**, 743 (2017).

[16] G. Duclos, C. Blanch-Mercader, V. Yashunsky, G. Salbreux, J. F. Joanny, J. Prost, and P. Silberzan, Spontaneous shear flow in confined cellular nematics, *Nat. Phys.* **14**, 728 (2018).

[17] S. Jain, V. M. L. Cachoux, G. H. N. S. Narayana, S. de Beco, J. D'Alessandro, V. Cellerin, T. Chen, M. L. Heuzé, P. Marcq, R.-M. Mège, A. J. Kabla, C. T. Lim, and B. Ladoux, The role of single-cell mechanical behavior and polarity in driving collective cell migration, *Nat. Phys.* **16**, 802 (2020).

[18] R. Farhadifar, J.-C. Röper, B. Aigouy, S. Eaton, and F. Jülicher, The influence of cell mechanics, cell-cell interactions,

and proliferation on epithelial packing, *Curr. Biol.* **17**, 2095 (2007).

[19] T. E. Angelini, E. Hannezo, X. Trepat, J. J. Fredberg, and D. A. Weitz, Cell migration driven by cooperative substrate deformation patterns, *Phys. Rev. Lett.* **104**, 168104 (2010).

[20] T. E. Angelini, E. Hannezo, X. Trepat, M. Marquez, J. J. Fredberg, and D. A. Weitz, Glass-like dynamics of collective cell migration, *Proc. Natl. Acad. Sci. USA* **108**, 4714 (2011).

[21] M. Basan, J. Elgeti, E. Hannezo, W.-J. Rappel, and H. Levine, Alignment of cellular motility forces with tissue flow as a mechanism for efficient wound healing, *Proc. Natl. Acad. Sci. USA* **110**, 2452 (2013).

[22] R. Sunyer, V. Conte, J. Escribano, A. Elosegui-Artola, A. Labernadie, L. Valon, D. Navajas, J. M. García-Aznar, J. J. Muñoz, P. Roca-Cusachs, and X. Trepat, Collective cell dura-taxis emerges from long-range intercellular force transmission, *Science* **353**, 1157 (2016).

[23] B. Stramer and R. Mayor, Mechanisms and *in vivo* functions of contact inhibition of locomotion, *Nat. Rev. Mol. Cell Biol.* **18**, 43 (2017).

[24] R. A. Desai, S. B. Gopal, S. Chen, and C. S. Chen, Contact inhibition of locomotion probabilities drive solitary versus collective cell migration, *J. R. Soc. Interface* **10**, 20130717 (2013).

[25] E. Scarpa, A. Roycroft, E. Theveneau, E. Terriac, M. Piel, and R. Mayor, A novel method to study contact inhibition of locomotion using micropatterned substrates, *Biol. Open* **2**, 901 (2013).

[26] B. Lin, T. Yin, Y. I. Wu, T. Inoue, and A. Levchenko, Interplay between chemotaxis and contact inhibition of locomotion determines exploratory cell migration, *Nat. Commun.* **6**, 6619 (2015).

[27] David Li and Yu-Li Wang, Coordination of cell migration mediated by site-dependent cell-cell contact, *Proc. Natl. Acad. Sci. USA* **115**, 10678 (2018).

[28] B. Smeets, R. Alert, J. Pešek, I. Pagonabarraga, H. Ramon, and R. Vincent, Emergent structures and dynamics of cell colonies by contact inhibition of locomotion, *Proc. Natl. Acad. Sci. USA* **113**, 14621 (2016).

[29] J. Zimmermann, B. A. Camley, W.-J. Rappel, and H. Levine, Contact inhibition of locomotion determines cell–cell and cell–substrate forces in tissues, *Proc. Natl. Acad. Sci. USA* **113**, 2660 (2016).

[30] V. E. Debets, L. M. C. Janssen, and C. Storm, Enhanced persistence and collective migration in cooperatively aligning cell clusters, *Biophys. J.* **120**, 1483 (2021).

[31] J. Toner, Y. Tu, and S. Ramaswamy, Hydrodynamics and phases of flocks, *Ann. Phys.* **318**, 170 (2005).

[32] M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao, and R. A. Simha, Hydrodynamics of soft active matter, *Rev. Mod. Phys.* **85**, 1143 (2013).

[33] B. A. Camley and W.-J. Rappel, Physical models of collective cell motility: From cell to tissue, *J. Phys. D* **50**, 113002 (2017).

[34] B. Szabó, G. J. Szöllösi, B. Gönci, Zs. Jurányi, D. Selmeczi, and T. Vicsek, Phase transition in the collective migration of tissue cells: Experiment and model, *Phys. Rev. E* **74**, 061908 (2006).

[35] J. M. Belmonte, G. L. Thomas, L. G. Brunnet, R. M. C. de Almeida, and H. Chaté, Self-propelled particle model for cell-sorting phenomena, *Phys. Rev. Lett.* **100**, 248702 (2008).
[36] S. Henkes, Y. Fily, and M. C. Marchetti, Active jamming: Self-propelled soft particles at high density, Phys. Rev. E **84**, 040301(R) (2011).

[37] S. Garcia, E. Hannezo, J. Elgeti, J.-F. Joanny, P. Silberzan, and N. S. Gov, Physics of active jamming during collective cellular motion in a monolayer, Proc. Natl. Acad. Sci. USA **112**, 15314 (2015).

[38] J. d'Alessandro, A. P. Solon, Y. Hayakawa, C. Anjard, F. Detcheverry, J.-P. Rieu, and C. Rivière, Contact enhancement of locomotion in spreading cell colonies, Nat. Phys. **13**, 999 (2017).

[39] D. Bi, J. H. Lopez, J. M. Schwarz, and M. L. Manning, A density-independent rigidity transition in biological tissues, Nat. Phys. **11**, 1074 (2015).

[40] D. Bi, X. Yang, M. C. Marchetti, and M. L. Manning, Motility-driven glass and jamming transitions in biological tissues, Phys. Rev. X **6**, 021011 (2016).

[41] T. B. Saw, A. Doostmohammadi, V. Nier, L. Kocgozlu, S. Thampi, Y. Toyama, P. Marcq, C. T. Lim, J. M. Yeomans, and B. Ladoux, Topological defects in epithelia govern cell death and extrusion, Nature (London) **544**, 212 (2017).

[42] C. Pérez-González, R. Alert, C. Blanch-Mercader, M. Gómez-González, T. Kolodziej, E. Bazellieres, J. Casademunt, and X. Trepat, Active wetting of epithelial tissues, Nat. Phys. **15**, 79 (2019).

[43] B. A. Camley, Y. Zhang, Y. Zhao, B. Li, E. Ben-Jacob, H Levine, and W.-J. Rappel, Polarity mechanisms such as contact inhibition of locomotion regulate persistent rotational motion of mammalian cells on micropatterns, Proc. Natl. Acad. Sci. USA **111**, 14770 (2014).

[44] G. Peyret, R. Mueller, J. d'Alessandro, S. Begnaud, P. Marcq, R.-M. Mêge, J. M. Yeomans, A. Doostmohammadi, and B. Ladoux, Sustained oscillations of epithelial cell sheets, Biophys. J. **117**, 464 (2019).

[45] Youyuan Deng, Herbert Levine, Xiaoming Mao, and Leonard M. Sander, Collective motility and mechanical waves in cell clusters, Eur. Phys. J. E **44**, 137 (2021).

[46] T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen, and O. Shochet, Novel type of phase transition in a system of self-driven particles, Phys. Rev. Lett. **75**, 1226 (1995).

[47] M. E. Cates and J. Tailleur, Motility-induced phase separation, Annu. Rev. Condens. Matter Phys. **6**, 219 (2015).

[48] E. Tjhung, C. Nardini, and M. E. Cates, Cluster phases and bubbly phase separation in active fluids: Reversal of the Ostwald process, Phys. Rev. X **8**, 031080 (2018).

[49] S. R. K. Vedula, A. Ravasio, E. Anon, T. Chen, G. Peyret, M. Ashraf, and B. Ladoux, Microfabricated environments to study collective cell behaviors, in Micropatterning in Cell Biology Part B, Methods in Cell Biology, Vol. 120, edited by M. Piel and M. Théry (Academic Press, Cambridge, MA, 2014), Chap. 16, pp. 235–252.

[50] M. Machacek, L. Hodgson, C. Welch, H. Elliott, O. Pertz, P. Nalbant, A. Abell, G. L. Johnson, K. M. Hahn, and G. Danuser, Coordination of Rho GTPase activities during cell protrusion, Nature (London) **461**, 99 (2009).

[51] M. Abercrombie, Contact inhibition and malignancy, Nature (London) **281**, 259 (1979).

[52] See Supplemental Material at http://link.aps.org/supplemental/10.1103/PhysRevResearch.6.023022 for movies S1–S4, and additional details about the experimental methods and the numerical simulations as well as the active hydrodynamic theory.

[53] L. P. Dadhichi, A. Maitra, and S. Ramaswamy, Origins and diagnostics of the nonequilibrium character of active systems, J. Stat. Mech.: Theory Exp. (2018) 123201.

[54] Vicsek-like models break the invariance under independent rotations of space and polarity vectors, but only dynamically.

[55] S. Dhakal and J. V. Selinger, Statistical mechanics of splay flexoelectricity in nematic liquid crystals, Phys. Rev. E **81**, 031704 (2010).

[56] F. Peruani, A. Deutsch, and M. Bär, A mean-field theory for self-propelled particles interacting by velocity alignment mechanisms, Eur. Phys. J.: Spec. Top. **157**, 111 (2008).

[57] M. George, F. Bullo, and O. Campàs, Connecting individual to collective cell migration, Sci. Rep. **7**, 9720 (2017).

[58] J. Zhang, R. Alert, J. Yan, N. S. Wingreen, and S. Granick, Active phase separation by turning towards regions of higher density, Nat. Phys. **17**, 961 (2021).

[59] E. Bertin, M. Droz, and G. Grégoire, Boltzmann and hydrodynamic description for self-propelled particles, Phys. Rev. E **74**, 022101 (2006).

[60] P. G. de Gennes and J. Prost, The Physics of Liquid Crystals (Oxford University Press, Oxford, UK, 1993).

[61] W. Kung, M. C. Marchetti, and K. Saunders, Hydrodynamics of polar liquid crystals, Phys. Rev. E **73**, 031708 (2006).

[62] L. P. Dadhichi, J. Kethapelli, R. Chajwa, S. Ramaswamy, and A. Maitra, Nonmutual torques and the unimportance of motility for long-range order in two-dimensional flocks, Phys. Rev. E **101**, 052601 (2020).

[63] D. Geyer, D. Martin, J. Tailleur, and D. Bartolo, Freezing a flock: Motility-induced phase separation in polar active liquids, Phys. Rev. X **9**, 031043 (2019).

[64] Miguel Vicente-Manzanares, Donna J. Webb, and A. Rick Horwitz, Cell migration at a glance, J. Cell Sci. **118**, 4917 (2005).

[65] M. E. Brasch, G. Passucci, A. C. Gulvady, C. E. Turner, M. L. Manning, and J. H. Henderson, Nuclear position relative to the Golgi body and nuclear orientation are differentially responsive indicators of cell polarized motility, PLoS One **14**, e0211408 (2019).

[66] K. Vaidžiulytė, A.-S. Macé, A. Battistella, W. Beng, K. Schauer, M. Coppey, Persistent cell migration emerges from a coupling between protrusion dynamics and polarized trafficking, eLife **11**, e69229 (2022).

[67] http://site.physics.georgetown.edu/matlab/.

[68] https://www.mn.uio.no/math/english/people/aca/jks/matpiv/.

[69] R. J. Glauber, Time-dependent statistics of the Ising model, J. Math. Phys. **4**, 294 (1963).

[70] A. C. Brańka and D. M. Heyes, Algorithms for Brownian dynamics computer simulations: Multivariable case, Phys. Rev. E **60**, 2381 (1999).