Collective epithelial migration mediated by the unbinding of hexatic defects

eLife Assessment

This important theoretical study shows that active hexatic topological defects in epithelia enable collective cell flows. Within the general limitations of coarse-grained hydrodynamic models in fully capturing cell-scale behavior, the study provides compelling evidence supporting its conclusions. These findings will be of interest to both biophysicists studying collective cell behaviors and biologists investigating epithelial flows during development.

https://doi.org/10.7554/eLife.105397.3.sa0

Significance of the findings:

Important: Findings that have theoretical or practical implications beyond a single subfield

Landmark
Fundamental
Important
Valuable
Useful

Strength of evidence:

Compelling: Evidence that features methods, data and analyses more rigorous than the current state-of-the-art

Exceptional
Compelling
Convincing
Solid
Incomplete
Inadequate

During the peer-review process the editor and reviewers write an eLife Assessment that summarises the significance of the findings reported in the article (on a scale ranging from landmark to useful) and the strength of the evidence (on a scale ranging from exceptional to inadequate). Learn more about eLife Assessments

Abstract
Introduction
Results
Discussion
Methods
Data availability
References
Article and author information
Metrics

Abstract

Collective cell migration in epithelia relies on cell intercalation: a local remodeling of the cellular network that allows neighboring cells to swap their positions. Unlike foams and passive cellular fluid, in epithelial intercalation, these rearrangements crucially depend on activity. During these processes, the local geometry of the network and the contractile forces generated therein conspire to produce a burst of remodeling events, which collectively give rise to a vortical flow at the mesoscopic length scale. In this article, we formulate a continuum theory of the mechanism driving this process, built upon recent advances toward understanding the hexatic (i.e., sixfold ordered) structure of epithelial layers. Using a combination of active hydrodynamics and cell-resolved numerical simulations, we demonstrate that cell intercalation takes place via the unbinding of topological defects, naturally initiated by fluctuations and whose late-times dynamics is governed by the interplay between passive attractive forces and active self-propulsion. Our approach sheds light on the structure of the cellular forces driving collective migration in epithelia and provides an explanation of the observed extensile activity of in vitro epithelial layers.

Introduction

From humble soap froths (Weaire and Hutzler, 1999; Graner et al., 2008; Marmottant et al., 2008) down to epithelial layers, confluent cellular fluids use intercalation to flow, even in the absence of gaps and interstitial structures. At the heart of the process is a mechanism known as topological rearrangement process of the first kind, or T1 for brevity, through which the vertices of a honeycomb network merge and then split, thereby leading to a remodeling of the network’s topology. In epithelia, this strategy allows cells to migrate, while preserving the structural integrity of the tissue, hence its major biological functionalities. These include the barrier function, which allows epithelial layers to maintain homeostasis, ensure nutrient transport, and filter out harmful pathogens, as well as extrusion and replacement of apoptotic cells (Irvine and Wieschaus, 1994; Keller et al., 2000; Walck-Shannon and Hardin, 2014; Tetley et al., 2016; Tetley and Mao, 2018; Paré and Zallen, 2020; Rauzi, 2020).

While relying exclusively on local rearrangements of the cellular network, intercalation gives rise to surprisingly organized patterns, where cells are able to migrate collectively over distances orders of magnitude larger than the average cell size. In the morphogenesis of Drosophila, for instance, cell intercalation drives a major cellular rearrangement known as germ-band extension, in which a layer of cells initially localized in the ventral region of the developing embryo folds over the dorsal region upon extending by approximately two and a half times its initial length. In this process, individual cells persistently move across the anterior-posterior axis more than 10 times their body length. In Drosophila, intercalation is also the main process driving the salivary gland tube formation, in which cells radially converge toward a central pit and eventually escape from the tangent plane of the embryo, thereby giving rise to the tubular structure (Sanchez-Corrales et al., 2018).

In the context of cancer progression, the role of cell intercalation has been recently debated in relation to various stages of the so-called metastatic cascade: that is, biomechanical pathway leading to the formation of a secondary tumor (Cheung and Ewald, 2016). The latter is schematically divided into three main phases: (1) detachment of cell clusters from a primary tumor and invasion of the extracellular matrix (ECM); (2) penetration (i.e., intravasation), circulation, and expulsion (i.e., extravasation) of the clusters in and from the blood stream; and (3) colonization of a healthy tissue and the proliferation of a secondary tumor. Along the cascade, metastatic cells undergo multiple phenotypic switches, aimed at maximizing their chances of success and survival within the surrounding microenvironment. These, in turn, determine the cells’ motility mode, which can vary from individual (e.g., ameboid or mesenchymal) to collective (e.g., intercalation-based or flocking guided by a small number leader cells localized at the front of the cluster) (Friedl and Gilmour, 2009; Haeger et al., 2020; Serra-Picamal et al., 2012; Murugan et al., 2024).

Yet, while being highly regulated at the biochemical level (Cavey et al., 2008; Yamada et al., 2005; Yonemura et al., 2010; Buckley et al., 2014; Engl et al., 2014; Zallen and Wieschaus, 2004; Bertet et al., 2004), cell intercalation cannot be separated from the mechanical forces originating it and whose nature, spatial organization, and dynamics are still largely unknown. In this article, we provide a topological insight into the mechanics of cell intercalation, by leveraging recent advances toward deciphering orientational order in epithelial layers (Li and Ciamarra, 2018; Pasupalak et al., 2020; Durand and Heu, 2019; Armengol-Collado et al., 2023a; Armengol-Collado et al., 2023b; Eckert et al., 2023; Cislo et al., 2023). The latter originates from the cells’ anisotropic shape and results in the emergence of liquid crystal phases collectively known as p-atics, with $p$ an integer reflecting the symmetry of the system under rotation by $2 π / p$ . The honeycomb structure of epithelial layers, in particular, has been shown to give rise to hexatic order (i.e., $p = 6$ ) at length scales ranging from one to dozens of cells, depending on the cells’ density and molecular repertoire, as well as the mechanical properties of the substrate (Armengol-Collado et al., 2023a; Armengol-Collado et al., 2023b; Eckert et al., 2023). In the language of liquid crystals, the cell-wide morphological transformations underlying T1 processes can be described in terms of topological defects known as disclinations: that is, point-like singularities in the otherwise regular configuration of a continuous p-atic order parameter – that is, $Ψ_{p} = ⟨ e^{i p ϑ} ⟩$ , with $ϑ$ the orientation of the individual building blocks and $⟨ \dots ⟩$ the ensemble average (Giomi et al., 2022a; Giomi et al., 2022b) – around which the average cellular orientation rotates by $2 π s$ , with $s = \pm 1 / p, \pm 2 / p \dots$ the winding number or ‘strength’ of the defect. In passive two-dimensional matter, disclinations mediate the transition from solid to liquid via a process known as Kosterlitz–Thouless–Halperin–Nelson–Young (KTHNY) melting scenario (Kosterlitz and Thouless, 1972; Kosterlitz and Thouless, 1973; Nelson and Halperin, 1979; Young, 1979; Kosterlitz, 2016). According to this, the hierarchical unbinding of neutral defect complexes – that is, for which $\sum_{i} s_{i} = 0$ – renders the system progressively more disordered. Our working hypothesis is that the competition between active and passive forces drives a similar unbinding mechanism in epithelial layers. In the following, we clarify the various steps and possible outcomes of this process and test this hypothesis against both hydrodynamic and cell-resolved numerical simulations.

Results

Cell intercalation and T1 cycle

A typical cell intercalation is illustrated in Figure 1a for a cluster of four cells, hereafter referred to as primary cell cluster. A T1 occurs when the junction that connects the internal threefold coordinated vertices shrinks until they merge into a fourfold vertex and then split once more along the orthogonal direction. Such an internal T1, however, leaves the number and, importantly, the positions of external vertices of each cell unchanged. Thus, despite this concept having received little attention in the literature (see, e.g., Staple et al., 2010; Fletcher et al., 2014; Rauzi, 2020; Duclut et al., 2022; Sknepnek et al., 2023; Jain et al., 2023), it is impossible to achieve collective migration by means of isolated T1 processes.

Figure 1

Download asset Open asset

Cell intercalation and T1 cycle.

(a) A full cell intercalation consists of an *internal* and four *external* T1 processes. The latter reconfigures the peripheral vertices of the primary cluster, thereby triggering new T1 processes across the neighboring cells. In the language of topological defects, the T1 translates to the (i) unbinding of a $\pm 1 / 6$ defect quadrupole and (ii) a further unbinding of the quadrupole into a pair of dipoles. These two processes are schematically presented in a specific temporal order, but, in practice, they occur simultaneously or nearly so. (b) In a T1 cycle, the primary cell cluster undergoes a T1, followed by an inverse T1, which restores its initial configuration. The process corresponds to (i) the unbinding of a defect quadrupole and (ii) its annihilation.

A full cell intercalation then consists of a burst of internal and a peripheral T1 processes, as is schematically summarized in Figure 1a. These steps do not necessarily occur in a sequential order but are most often simultaneous. Furthermore, since collective migration is a phenomenon that occurs roughly homogeneously across the entire cell layer, there is no single initial T1, but a uniform distribution of seeds. Crucially, internal T1 processes do not always trigger a full intercalation. After the first T1, where the internal threefold vertices shrink to a fourfold vertex, the cluster may reverse its dynamics and return to its original configuration (see Figure 1b). In the following, we refer to this scenario as the T1 cycle: that is, a direct T1 followed by an inverse T1, which does not permanently alter the configuration of the honeycomb network. Intuitively, and as we show next, whether the initial T1 triggers a full intercalation, hence collective migration, or a cycle, is determined by the configuration of contractile stresses exerted by the cells, which, in turn, are modulated by the local geometry of the primary cluster.

As a starting point, we focus on the intermediate configuration of primary cell cluster comprising two orthogonal pairs of pentagonal and heptagonal cells, as shown in the central column of Figure 1. This configuration, which corresponds to the most elementary short-ranged excitation of a honeycomb network, can be described in the language of topological defects as a quadrupole of $\pm 1 / 6$ disclinations. As in two-dimensional melting, such a defect structure can arise spontaneously as a consequence of spatiotemporal fluctuations of physical or biological nature. The two complementary remodeling events following the intermediate configuration, that is, cell intercalation (Figure 1a) and the T1 cycle (Figure 1b), correspond instead to the two possible ‘fates’ of this initial excitation. Upon unbinding (Figure 1ai and bi), the four defects comprising the quadrupole can further unbind into two $\pm 1 / 6$ disclination dipoles (Figure 1aii), or annihilate (Figure 1bii), thereby restoring the initial defect-free configuration. As we demonstrate next, the former scenario corresponds to cell intercalation and the latter to the T1 cycle. To this end, we identify three geometric requirements that a model of cell intercalation must fulfill, regardless of the desired level of biophysical accuracy. (1) The average orientation of the cells must rotate by $π / 6$ with respect to initial configuration. (2) Both the primary cluster and its surrounding cells must perform a local convergent extension: that is, move inward in one direction, and outward along the orthogonal one (Keller et al., 2008; Blanchard, 2017; Wang et al., 2020; Ioratim-Uba et al., 2023). We stress that the adjective local is used here to distinguish this process from convergent extension as intended in developmental biology, where the same rearrangement occurs at the scale of the entire organism. (3) In order to remodel the external vertices and initiate cell intercalation, the primary cluster must undergo a spontaneous shear deformation.

In the following, we show that our construction not only fulfills these requirements, but, harnessing the predictive power of active hydrodynamics, provides readily testable experimental predictions. To this end, we numerically integrate the hydrodynamic equations of active hexatic liquid crystals, introduced by Armengol-Collado et al. in Armengol-Collado et al., 2023b (see ‘Methods’). To follow the fate of the primary cell cluster after an internal T1, we assume the cells to be initially horizontally oriented, so that the phase θ of the hexatic order parameter, $Ψ_{6} = | Ψ_{6} | e^{6 i θ}$ , is $θ = 0$ , and construct a configuration featuring a quadrupole of $\pm 1 / 6$ disclinations (see Figure 2ai and bi). Along one full loop encircling each of these defects, θ changes by $\pm π / 3$ , with the sign reflecting that of the defect’s winding number. Since the external boundary of the primary cluster consists of four half loops, this implies that θ varies in the range $- π / 6 \leq θ \leq π / 6$ around the quadruple, as indicated by the alternating blue and red tones in Figure 2aii and bii. Once the defect quadrupole breaks into two dipoles, this new orientation propagates from the boundary of the primary cluster into the space between the dipoles (see Figure 2aiii), while leaving the orientation of the cells in the exterior essentially undistorted. Thus, the unbinding of a $\pm 1 / 6$ defect quadrupole from a defect-free configuration and its breakup into two dipoles drives a $π / 6$ rotation of the cells between the dipoles (see Figure 2aiv), consistently with our first requirement. Conversely, if defects annihilate (see Figure 2biii), the cells’ initial orientation is restored after a transient orientational fluctuation (see Figure 2biv).

Figure 2

Download asset Open asset

Cell intercalation and T1 cycle as defect unbinding and annihilation.

(a) Cell intercalation. (i) Backflow velocity field generated during the unbinding of an *active*, hexatic defect quadrupole (Video 1). The three panels below show the orientation field associated with (ii) the quadruple in the initial configuration, (iii) as it unbinds in a pair of $\pm 1 / 6$ dipoles, and (iv) after the dipoles have moved outside of the region of interest, together with the corresponding configuration of the primary cluster (Video 2). As the dipoles move away from each other, the cells surrounding the primary cluster rotate clockwise (blue) and counterclockwise (red). (b) T1 cycle. (i–iv) Analogous sequence as in panel (a), but associated with the annihilation of the defect quadrupole. Notice that, in panel (iii), the direction of the flow is reversed. The details of the finite difference simulations can be found in ‘Methods’.

Figure 2—source data 1 Values of the hexatic orientation field and velocity field used to produce the panels in this figure, obtained from finite-difference numerical solutions of Equation 3a, Equation 3b, Equation 3c.: https://cdn.elifesciences.org/articles/105397/elife-105397-fig2-data1-v1.zip
Download elife-105397-fig2-data1-v1.zip

Video 1

Download asset

posterframe for video — Backflow velocity field (vector plot) generated during the unbinding of an active, hexatic defect quadrupole (color plot, magnitude of order parameter).

Video 2

Download asset

In order to address the second and third requirements, we look at the configuration of the velocity field $v$ , corresponding to the average velocity of the cells in the surroundings of the primary cluster. As well documented in the theoretical (Giomi et al., 2014; Giomi, 2015; Hoffmann et al., 2020) and experimental (Saw et al., 2017; Kawaguchi et al., 2017; Blanch Mercader et al., 2018; Balasubramaniam et al., 2021; Yashunsky et al., 2022) literature of active nematic liquid crystals, the distortion induced by topological defects drives a flow, whose structure and direction is determined by the defect’s strength and the magnitude of the active stresses collectively exerted by the cells. Immediately after unbinding, an approximated expression for the velocity of the flow caused by the defect quadrupole can be analytically calculated (see ‘Methods’). Calling $r = | r | (\cos ϕ e_{x} + \sin ϕ e_{y})$ the distance from the center of the primary cluster and $ℓ$ the cluster’s size, that is, the distance between the center of the central junction and the center of any cell in the cluster, this velocity is given by

v \approx \frac{120 α_{6} ℓ^{4}}{η} [(4 - 6 \log \frac{| r |}{ℓ} - 3 \frac{{| r |}^{2}}{ℓ^{2}}) \cos 6 ϕ] \frac{r}{| r |^{8}} .

The constant $α_{6}$ has dimensions of force over volume, embodying the active stresses exerted by the cells and modulated by the local hexatic order, and η the shear viscosity (see ‘Methods’ for details). Thus, in close proximity to the primary cluster, where $| r | / ℓ ⪆ 1$ , Equation 1 gives $v_{x} (x, 0) \sim α_{6} / (η ℓ^{4}) x$ and $v_{y} (y, 0) \sim - α_{6} / (η ℓ^{4}) y$ . In agreement with our second requirement, this is a typical structure of a stagnation flow, whose realization at the cellular scale is a local convergent extension. We stress that, consistently with the cooperative and mesoscale nature of cell intercalation, Equation 1 cannot be obtained from the mere superposition of the flows individually sourced by the defects, as a consequence of the non-linear dependence of the order parameter on the average cellular orientation. Finally, we notice that the specific direction of motion – that is, the sign of $v$ – depends solely on $α_{6}$ , which, in turn, is either positive or negative depending on whether the active stresses exerted within the cell layer are respectively contractile or extensile.

As the cellular layer starts remodeling and the defects comprising the quadrupole migrate away from their original position, a solution of the hydrodynamic equations becomes analytically inaccessible, but can be obtained from a numerical integration of the hydrodynamic equations and is displayed in Figure 2a and b, for two different $α_{6}$ values. When $α_{6}$ is large and negative, the quadrupole splits into two $\pm 1 / 6$ dipoles moving away from each other at an angle of approximately $5 π / 6$ (see ‘Methods’ for details). Such a shear deformation is further enhanced by the coupling between hexatic order and flow, which, in a way not dissimilar to flow alignment effects in nematic liquid crystals, drives a rotation of the local orientation (Gennes and Prost, 1993). This biases the unbinding dynamics of the defect dipoles, thereby setting, in concert with the passive Coulomb-like forces at play, the direction along which the dipoles move away from each other (Krommydas et al., 2023). Finally, switching off active stresses – that is, $α_{6} = 0$ , shown in Figure 2b – suppresses both defect unbinding and shear. Consistently, inverting the direction of active forces from extensile to contractile – that is, $α_{6} > 0$ , see ‘Methods’ – results in a speed up of defects annihilation, hence of the T1 cycle, but never leads to the unbinding of $\pm 1 / 6$ pairs, thus to the onset of cell intercalation.

We conclude this section by stressing that the scenario emerging from our hydrodynamic analysis, of which Equation 1 represents a central outcome, is further corroborated by recent numerical work by Erdemci-Tandogan and Manning, 2021 and Das et al., 2021, who, using a different cell-resolved model of epithelia – that is, the Vertex model (Honda and Eguchi, 1980) – showed that the rate of T1 processes enhances the fluidity of tissues. Consistently, the speed of the stagnation flow sourced by cell intercalation increases like $η^{- 1}$ , indicating that the faster cells intercalate, the smaller η, the more fluid is the epithelial layer.

Collective cell migration

Having demonstrated the viability of our hydrodynamic approach, we next outline a number of general predictions, the most striking of which is that collective epithelial migration, as it results from cell intercalation, is a process of activity-guided defect unbinding. To this end, we perform numerical simulations of the multiphase field model (MPF) (Loewe et al., 2020), which have been proved to capture various aspects of epithelial organization, including the recently found hexanematic multiscale order (Armengol-Collado et al., 2023a; Armengol-Collado et al., 2023b; Eckert et al., 2023). In this approach, cells are modeled as droplets of immiscible fluid phases, undergoing a persistent random walk. Each cell is characterized by a nominal propulsion speed $v_{0}$ , equal for all cells, and a fluctuating direction of motion, with rotational diffusion coefficient $D_{r}$ . A detailed description of the MPF model and numerical details is provided in ‘Methods’. In order to test the correlation between topological defects, T1 processes, and collective migration, we fix the speed $v_{0}$ and the total number of cells and vary the persistence length of the cells' trajectories by tuning the rotational diffusion coefficient $D_{r}$ . Figure 3a and b show a typical configuration of MPF simulations, colored according to the normalized hexatic and nematic longitudinal stress, respectively, experienced by cells in the tissue. Note that the negative value of the longitudinal stress corresponds to extensile stresses (see the ‘Methods’ section for details).

Figure 3

Download asset Open asset

Collective cell migration as defect unbinding in the multiphase field model (MPF).

(**a, b**) Color plots illustrating the longitudinal hexatic (a) and nematic (b) stresses in MPF simulations (refer to ‘Methods’). The color bar is normalized to the largest stress magnitude observed in the configuration. Notably, the stress is uniformly negative, reflecting the extensile characteristics of both hexatic and nematic stresses. (c) Example of a four-cell cluster as it undergoes a T1 process, together with (d) the reconstructed sixfold orientation field. The six-legged stars mark the local sixfold orientation of the cells (see ‘Methods’), while the red and blue dots denote the $+ 1 / 6$ and $- 1 / 6$ defects. For such a four-cell cluster in real epithelial cell monolayer, please see Armengol-Collado et al., 2023a. (e) Probability distribution of finding a T1 (red tones) and a random cell (yellow tones) at a given distance from a defect, for four different values of the rotational noise $D_{r}$ . The data indicate a prominent correlation between T1 process and topological defects. (f) The mean square displacement (m.s.d) of cells versus defect density computed over a time window of $Δ t = 25 \times 10^{3}$ iterations, chosen to match the typical duration of T1 events and defect lifetimes. We identify two distinct subpopulations of cells: ‘slow’ (blue tones), with no correlation to the local density, and ‘fast’ (yellow tones), located where the local defect density is higher. The former corresponds to cells undergoing a T1 cycle and the latter participating in cell intercalation, hence to collective cell migration. (g) Temporal statistics of tissue remodeling events in multiphase field simulations. Average time between two intercalation events (orange) and average period of a T1 cycle (green) versus the rotational diffusion coefficient $D_{r}$ . The box plot in the inset shows the statistics of events analyzed for the case at $D_{r} = 4 \times 10^{- 5}$ . (Pairwise comparisons were performed with the two-sided t-test: $^{* * *} p < 10^{- 3}$ ). In the main graph, error bars are reported as the first (bottom bar) and third (upper bar) quartile of the dataset.

Figure 3—source data 1 Numerical data displayed in panels (e), (f) and (g).: https://cdn.elifesciences.org/articles/105397/elife-105397-fig3-data1-v1.zip
Download elife-105397-fig3-data1-v1.zip

To track topological defects, we first polygonize the cells by detecting their contour and marking the vertices, as shown in Figure 3c. For each cell, we then compute the shape function

γ_{6} = \frac{\sum_{v = 1}^{V} | r_{v} |^{6} e^{6 i ϕ_{v}}}{\sum_{v = 1}^{V} | r_{v} |^{6}},

where $𝒓_{v}$ , with $v = 1, 2 \dots V$ is the positions of the vth vertex with respect to the centroid of the polygon and $ϕ_{v} = \arctan (y_{v} / x_{v})$ its orientation with respect to the x-axis. As demonstrated in Armengol-Collado et al., 2023b, the phase $A r g (γ_{6}) / 6$ of the shape function identifies the sixfold orientation of a cell and is represented as white six-legged stars in Figure 3d. The shape function is then coarse-grained over the length scale $R = 1.5 R_{c e l l}$ to reconstruct the shape parameter $Γ_{6} = ⟨ γ_{6} ⟩_{R}$ . Finally, topological defects are identified as singularities in the otherwise smoothly varying field $Γ_{6} = Γ_{6} (r)$ . T1 processes, conversely, are readily detected by tracking those events leading to a recombination of the change of neighbors of a given cell (see Figure 3c).

For each $D_{r}$ value, we reconstruct the probability of finding a T1 at a given distance from a defect (see Figure 3c and d) and compare it with that of an arbitrary cell. Both distributions are approximately Gaussian and have vanishing mean if expressed in terms of the signed distance $Δ x$ (see Figure 3e). Furthermore, for larger $D_{r}$ values, where the short correlation time renders the system more disordered and dense of defects, the two probability distributions are hardly distinguishable, with the probability of finding a defect in proximity of a T1 only slightly larger than that associated with an arbitrary cell. As $D_{r}$ is decreased, also the density of defects decreases, and the latter probability distribution becomes flatter and flatter, while the former remains unchanged, thereby confirming that T1 processes are de facto a realization of hexatic defect unbinding. To correlate structure and dynamics, we show in Figure 3f the mean squared displacement of T1 cells as a function of local defect density. These measurements were taken over a time interval $Δ t = 2.5 \times 10^{4}$ , which aligns with the characteristic duration of T1 events in our simulations (see ‘Multiphase field model of epithelial tissues’ for further details). Our data confirm the existence of two different classes: that is, ‘slow’ and ‘fast’ cells, respectively denoted by purple and orange tones. While for both classes of cells, the mean squared displacement is roughly uniform for all values of the local defect density, fast cells only appear where the density is higher, thus providing an alternative signature of cell intercalation and T1 cycles. Cells undergoing a T1 cycle oscillate about their initial position but do not participate in collective migration and, therefore, exhibit small mean square displacement. By contrast, intercalating cells drive collective migration and, based on the correspondence identified here, can only be found in regions of high defect densities. Note that the difference in the mean square displacement of cells undergoing intercalation is at least one order of magnitude larger than that of cells affected by an isolated T1 cycle (Figure 3f). Hence, although isolated internal T1s (and by extension T1 cycles) can have small long-ranged effects, those effects are negligible in comparison to the collective long-ranged motion induced by a full cell intercalation.

Finally, following Das et al., 2021; Erdemci-Tandogan and Manning, 2021 we study the temporal correlations of T1 cycles and cell intercalation. Figure 3g shows a plot of the average time between two intercalation events and the average period of T1 cycles, for varying rotational diffusion coefficients $D_{r}$ . The former displays a decreasing trend with $D_{r}$ , indicating that rotational noise improves the performance of cell intercalation, thus favoring a faster collective migration. Conversely, the dynamics of T1 cycles is essentially unaffected by rotational diffusion. Both behaviors can be readily rationalized on the basis of our hydrodynamic description. As previously observed in relation to Figure 3e, rotational noise favors the proliferation of hexatic defects, hence the intercalation rate. Conversely, being defect annihilation driven by the passive attractive forces resulting from the entropic elasticity of the hexatic phase, once a T1 cycle is initiated, its duration is essentially independent of the noise strength.

Before concluding, we discuss three especially striking aspects of collective cell migration highlighted by our approach and amenable to experimental scrutiny. First, from the active flow given by Equation 1 and neglecting irrelevant inertial effects, one can estimate the magnitude of the active forces at play: that is, $F_{a c t i v e} \sim α_{6} / ℓ^{3}$ , where $ℓ$ denotes the range of the distortion caused by the unbinding dislocations. Such a scaling form results primarily from the sixfold structure of cellular forces under the assumption – still to be verified in hexatic epithelial layers, but consistent with previous observations on nematic cell cultures (Duclos et al., 2017; Saw et al., 2017; Kawaguchi et al., 2017; Blanch Mercader et al., 2018; Yashunsky et al., 2022) – that $α_{6}$ is, at least approximately, spatially uniform. Second, in order for these forces to prompt defect unbinding, thereby triggering cell intercalation, they must overcome the passive Coulomb-like forces driving their annihilation. At the length scale of the primary cluster, the magnitude of the latter is roughly given by $F_{p a s s i v e} \sim L_{6} / ℓ$ , with $L_{6}$ the orientational stiffness of the hexatic phase (see ‘Methods’). Equating $F_{a c t i v e}$ and $F_{p a s s i v e}$ provides an estimate of the typical size of the primary cluster at which the cellular layer becomes unstable to collective migration: that is, $ℓ_{6} = \sqrt{| α_{6} | / L_{6}}$ . This active hexatic length scale was identified in Armengol-Collado et al., 2023b and is believed to play a role analog to that of the active nematic length scale $ℓ_{2} = \sqrt{L_{2} / | α_{2} |}$ (Giomi, 2015). The latter is the fundamental parameter controlling the collective behavior of active nematic liquid crystals and, depending on how it compares with other extrinsic and intrinsic length scales, determines the hydrodynamic stability (Duclos et al., 2018) of active nematics, the distribution of the vortex area in chaotic cytoskeletal flows (Guillamat et al., 2017; Lemma et al., 2019), the size of nematic domains in bacterial colonies (You et al., 2018) and in vitro cultures of spindle-like cells (Blanch Mercader et al., 2018), etc. Similarly, we expect that confining epithelia at length scales smaller than $ℓ_{6}$ has the effect of suppressing cell intercalation and possibly rendering necessary different locomotion strategies and possibly favoring a switch to mesenchymal phenotypes. Finally, although actomyosin networks render individual cells contractile (Schwarz and Safran, 2002), at a collective level, epithelial (Saw et al., 2017; Blanch Mercader et al., 2018; Balasubramaniam et al., 2021) and neural progenitor monolayers behave as an extensile system (Kawaguchi et al., 2017). In agreement with this recent experimental evidence (Saw et al., 2017; Blanch Mercader et al., 2018; Balasubramaniam et al., 2021), our analysis shows that only an extensile activity can fuel collective migration in confluent epithelia and, leveraging on the familiar language of topological defects, it further provides a simple key to rationalizing the mechanical advantage of such a biological strategy: that is, only extensile forces can provide the type of repulsive interactions that are necessary for defects to unbind, thus to trigger cell intercalation.

Discussion

In conclusion, we have investigated the physical mechanisms behind cell intercalation in confluent epithelial layers using analytics, and a combination of continuum and discrete modeling. After having established that cell intercalation, hence collective migration, requires a burst of correlated T1 processes, we demonstrated how the latter originates from the unbinding of neutral quadrupoles of hexatic disclinations. As in the KTHNY melting scenario (Nelson and Halperin, 1979), the latter are spontaneously generated by fluctuations and can either annihilate, thereby restoring the original configuration (T1 cycles), or further unbind in two pairs of $\pm 1 / 6$ disclinations, decreasing the translational and orientation order and increasing its fluidity. As these pairs move away from each other by activity, they stir and shear the cellular network, thereby initiating further T1 processes, which, cooperatively, give rise to cell migration. To assess the significance of our predictions, we have tested them against numerical simulations of the MPF (Loewe et al., 2020; Monfared et al., 2023; Camley et al., 2014; Palmieri et al., 2015; Peyret et al., 2019; Alert and Trepat, 2020), finding good quantitative agreement. We stress, however, that the theoretical framework employed here is not specifically tailored to the MPF or other discrete models of epithelial layers – a choice that would entail the risk of constructing a ‘model of a model’, rather than a model of the physical system itself. Conversely, by taking a more generic top-down approach, our theory sheds light on the structure of the cellular forces driving collective migration in epithelia, suggests the possibility of a confinement-induced switch to mesenchymal phenotypes, and provides an explanation of the observed extensile activity of in vitro epithelial layers (Saw et al., 2017; Blanch Mercader et al., 2018; Balasubramaniam et al., 2021).

Methods

Hydrodynamic equations of epithelial layers

Our hydrodynamic equations of epithelial layers have been given in Armengol-Collado et al., 2023b and account for both hexatic and nematic order, with the former being dominant at short and the latter at long length scales. Since cell intercalation occurs at small length scales, here we can ignore nematic order and focus solely on the hydrodynamics of the hexatic phase. In addition to the standard density ρ and velocity $v$ , this can be described in terms of the sixfold order parameter tensor $𝑸_{6} = 4 ⟦ ⟨ 𝝂^{\otimes 6} ⟩ ⟧ = 4 | Ψ_{6} | ⟦ 𝒏^{\otimes 6} ⟧$ . In this expression, $^{\otimes 6}$ indicates the sixfold tensorial product, $Q_{6} = 4 [[⟨ ν^{\otimes 6} ⟩]] = 4 | Ψ_{6} | [[n^{\otimes 6}]]$ and $n = \cos θ e_{x} + \sin θ e_{y}$ are respectively the fluctuating and average orientation of the hexatic building blocks and the operator $⟦ \dots ⟧$ renders its argument traceless and symmetric (see Giomi et al., 2022a; Giomi et al., 2022b for a general introduction to p-atic hydrodynamics). The short-scale hydrodynamic equations of the epithelial layer are then given, in the most generic form, by

\frac{D ρ}{D t} + ρ \nabla \cdot v = (k_{d} - k_{a}) ρ,

ρ \frac{D v}{D t} = \nabla \cdot σ + f,

\frac{D Q_{6}}{D t} = Γ_{6} H_{6} + 6 [Q_{6} \cdot ω] + λ_{6} [\nabla^{\otimes 4} u] + \bar{λ_{6}} t r (u) Q_{6} + v_{6} [u^{\otimes 3}] .

with $ν = \cos ϑ e_{x} + \sin ϑ e_{y}$ is the material derivative. In Equation 3a, $k_{d}$ and $k_{a}$ the cell division and apoptosis rates, here assumed equal. For simplicity, we also assume uniform density throughout the system so that Equation 3a reduces to the standard incompressibility condition $\nabla \cdot v = 0$ . In Equation 3b, $σ$ is the total stress tensor and $[[\dots]]$ an external body force. In Equation 3c, $𝒖 = [\nabla 𝒗 + {(\nabla 𝒗)}^{⊺}] / 2$ and $ω = [\nabla v - (\nabla v)^{⊺}] / 2$ , with $⊺$ indicating transposition, are respectively the strain rate and vorticity tensors and entail the coupling between hexatic order and flow, with $λ_{6}$ and $ν_{6}$ material constants and ${(\nabla^{\otimes n})}_{i_{1} i_{2} . . . i_{n}} = \partial_{i_{1}} \partial_{i_{2}} \dots \partial_{i_{n}}$ . Because of incompressibility, the term ${\bar{λ}}_{6} tr (𝒖) 𝑸_{6}$ in Equation 3c vanishes. The tensor $ω = [\nabla v - (\nabla v)^{⊺}] / 2$ is the hexatic analog of the molecular tensor, dictating the relaxation dynamics of the order parameter tensor toward the minimum of the orientational free energy

F = \int d A (\frac{L_{6}}{2} | \nabla Q_{6} |^{2} + \frac{A_{6}}{2} | Q_{6} |^{2} + \frac{B_{6}}{4} | Q_{6} |^{4}),

where $| \dots |^{2}$ is the Euclidean norm and is such that $| Q_{6} |^{2} = | Ψ_{6} |^{2} / 2$ . The constant $L_{6}$ is the order parameter stiffness, while $A_{6}$ and $B_{6}$ are phenomenological constants setting the magnitude of the coarse-grained complex order parameter at equilibrium: $| Ψ_{6} | = | Ψ_{6}^{(0)} | = \sqrt{- 2 A_{6} / B_{6}}$ , when $H_{i_{1} i_{2} \dots i_{6}} = 0$ .

The stress tensor figuring in Equation 3b can be customarily decomposed into a passive and an active contribution: that is, $𝝈 = 𝝈^{(p)} + 𝝈^{(a)}$ . The passive stress, in turn, can be expressed as $σ^{(p)} = - P 1 + σ^{(v)} + σ^{(e)} + σ^{(r)}$ , where $P$ is the pressure and $𝝈^{(v)} = 2 η ⟦ 𝒖 ⟧$ the viscous stress, with η the shear viscosity. The tensor $σ_{i j}^{(e)} = - L_{6} \partial_{i} Q_{6} ⊙ \partial_{j} Q_{6}$ is the elastic stress, arising in response to a static deformation of a fluid patch and the symbol $⊙$ indicates a contraction of all matching indices of the two operands yielding a tensor whose rank equates the number of unmatched indices (two in this case). Finally $𝝈^{(r)} = - λ_{6} \nabla^{\otimes 4} ⊙ 𝑯_{6} + 3 (𝑸_{6} \cdot 𝑯_{6} - 𝑯_{6} \cdot 𝑸_{6})$ is the reactive stress tensor, which embodies the conservative forces arising in response to flow-induced distortions of the hexatic orientation. The active stress tensor $σ^{(a)}$ was introduced in Armengol-Collado et al., 2023b on the basis of phenomenological and microscopic arguments and is given by $𝝈^{(a)} = α_{6} \nabla^{\otimes 4} ⊙ 𝑸_{6}$ , with $α_{6}$ a constant.

Multiphase field model for epithelial tissues

The model

The multiphase field model is a cell-resolved model where each cell is described in a two-dimensional space by a concentration field $φ_{c} = φ_{c} (r)$ , with $c = 1, 2 \dots N_{c e l l}$ , and $N_{c e l l}$ the total number of cells in the systems (Loewe et al., 2020; Monfared et al., 2023). The equilibrium state is defined by the free energy $F = \int d A f$ , where the free energy density $f$ is

f = \frac{α}{4} \sum_{c} φ_{c}^{2} (φ_{c} - φ_{0})^{2} + \frac{k_{φ}}{2} \sum_{c} | \nabla φ_{c} |^{2} + ϵ \sum_{c < c^{'}} φ_{c}^{2} φ_{c^{'}}^{2} + k_{a d} \sum_{c < c^{'}} \nabla φ_{c} \cdot \nabla φ_{c^{'}} + \sum_{c} λ {(1 - \frac{1}{π φ_{0}^{2} R_{c e l l}^{2}} \int d A φ_{c}^{2})}^{2} .

The first two terms in the free energy represent a $ϕ^{4}$ theory for phase separation. The parameters $α > 0$ and $k_{φ} > 0$ control the nominal surface tension $σ = \sqrt{8 k_{φ} α}$ and interfacial thickness $ξ = \sqrt{2 k_{φ} / α}$ , stabilizing spherical shapes in an isolated environment. This ensures the concentration field $φ_{c}$ remains close to $φ_{0}$ inside the cell and vanishes outside. The term proportional to $ϵ > 0$ enforces volume exclusion, preventing cell overlap, while the $k_{a d} > 0$ term models cell–cell adhesion, promoting tissue confluency in dense environments (Monfared et al., 2023). Additionally, the $λ > 0$ term constrains cell area near its nominal value $π R_{c e l l}^{2}$ , where $R_{c e l l}$ is the preferred cell radius.

The phase field $φ_{c}$ evolves via the Allen–Cahn equation:

\partial_{t} φ_{c} + v_{c} \cdot \nabla φ_{c} = - M \frac{δ F}{δ φ_{c}},

where non-equilibrium effects arise from the non-equilibrium advection term involving $𝒗_{c} = v_{0} (\cos θ_{c} 𝒆_{x} + \sin θ_{c} 𝒆_{y})$ . Here, $v_{0}$ is the constant self-propulsion speed of the cth cell, and $θ_{c}$ defines its migration direction. The angle $θ_{c}$ , in turn, follows a stochastic dynamics:

\frac{d θ_{c}}{d t} = η_{c},

with noise $η_{c}$ satisfying $⟨ η_{c} (t) η_{c^{'}} (t^{'}) ⟩ = 2 D_{r} δ_{c c^{'}} δ (t - t^{'})$ , where $D_{r}$ sets the noise diffusivity. The mobility $M$ in Equation 8 governs the relative strength of thermodynamic relaxation versus non-equilibrium migration.

We stress that, while multiphase field models can incorporate non-equilibrium effects in various ways, here the sole source of non-equilibrium is the self-propulsion velocity $v_{c}$ .

Here, $α > 0$ and $k_{φ} > 0$ are material parameters which can be used to tune the surface tension $σ = \sqrt{8 k_{φ} α}$ and the interfacial thickness $ξ = \sqrt{2 k_{φ} / α}$ of isolated cells and thermodynamically favor spherical cell shapes, so that the concentration field is large (i.e., $φ_{c} \approx φ_{0}$ ) inside the cells and zero outside. The repulsive bulk term proportional to $ϵ > 0$ captures the fact that cells cannot overlap, while the term proportional to $k_{a d} > 0$ models the interfacial adhesion between cells, ultimately favoring tissue confluency in a crowded environment (Monfared et al., 2023). The term proportional to $λ > 0$ forces the cells’ area around its nominal value $π R_{c e l l}^{2}$ , with $R_{c e l l}$ the preferential cell radius. The phase field $φ_{c}$ evolves according to the Allen–Cahn equation

\partial_{t} φ_{c} + v_{c} \cdot \nabla φ_{c} = - M \frac{δ F}{δ φ_{c}},

where $𝒗_{c} = v_{0} (\cos θ_{c} 𝒆_{x} + \sin θ_{c} 𝒆_{y})$ is the velocity at which the cth cell self-propels, with $v_{0}$ a constant speed and $θ_{c}$ the angle defining the nominal direction of cell migration. The latter evolves according to the stochastic equation

\frac{d θ_{c}}{d t} = η_{c},

where $η_{c}$ is a noise term with correlation function $⟨ η_{c} (t) η_{c^{'}} (t^{'}) ⟩ = 2 D_{r} δ_{c c^{'}} δ (t - t^{'})$ and $D_{r}$ a constant controlling noise diffusivity. The constant $M$ in Equation 8 is the mobility measuring the relevance of thermodynamic relaxation with respect to non-equilibrium cell migration. Equation 8 is solved with a finite-difference approach through a predictor-corrector finite difference Euler scheme implementing second order stencil for space derivatives (Carenza et al., 2019).

We have integrated the dynamical equations with a number of cells $N_{c e l l s} = 440$ in a system of size 384 × 403 with periodic boundary conditions. The model parameters values used in the simulations are as follows: $α = 0.2$ , $k_{φ} = 0.2$ , $ϵ = 0.1$ , $k_{a d} = 0.005$ , $λ = 600$ , $R_{c e l l} = 11.5$ , $v_{0} = 0.006$ . The noise variance $D_{r}$ was varied in the range $10^{- 6} \leq D_{r} < 8 \times 10^{- 5}$ .

Cell segmentation

In order to find the sixfold orientation of each cell, we proceeded to the cell segmentation of the simulated configuration. This procedure consists of the following steps. First, we define the thresholded density of the whole cell layer as

Φ = \sum_{c = 1}^{N_{c e l l}} ϑ_{H} (φ_{c} - φ_{t h}),

where $ϑ_{H}$ is the Heaviside theta function such that $ϑ_{H} (x) = 1$ if $x \geq 0$ and $ϑ_{H} (x) = 0$ otherwise. Here, $φ_{t h}$ is the threshold marking the boundary of each cell. In particular, we choose $φ_{t h} = φ_{0} / 2$ . The resulting field is $Φ = 1$ inside each cell; $Φ = 2$ at the interface of two cells; $Φ = 3$ or $Φ = 4$ on the vertices. As at the interface the field $φ_{c}$ of each cell smoothly changes from $φ_{0}$ to 0, therefore dropping below the threshold $φ_{t h}$ , it is possible to find pixels where $Φ = 0$ . These spurious features are adjusted by replacing Equation 10 with an average over the pixels neighboring that where $Φ$ vanishes. That is,

Φ = \sum_{⟨ x, y ⟩} \sum_{c = 1}^{N_{c e l l}} \frac{ϑ_{H} (φ_{c})}{max {\sum_{⟨ x, y ⟩} ϑ_{H} (φ_{c}), 1}},

where $\sum_{⟨ x, y ⟩}$ stands for a sum over the nearest neighbors of the pixel where $Φ = 0$ . The procedure is reiterated until no change occurs in two consecutive iterations. Finally, upon segmenting the thresholded density, we identify the cell’s vertices ${r_{v}}_{c}$ , as those points where $φ_{t o t} > 2$ . Tissue rearrangement events are identified by tracking changes in the list of neighbors of each cell.

Cell orientation, coarse-graining, and topological defects

The sixfold orientation of a cell is computed via the shape function $γ_{6}$ defined in Equation 2, introduced in Armengol-Collado et al., 2023a. This construction is schematically illustrated in Figure 4. The single-cell orientation can then be coarse-grained to construct a continuous description of the cellular tissue (Armengol-Collado et al., 2023a). To do so, we use the shape order parameter $Γ_{6} = Γ_{6} (𝒓)$ , constructed upon averaging the shape function $γ_{6}$ of the segmented cells whose center of mass, $r_{c}$ , lies within a disk of radius $R$ centered in $𝒓$ . That is,

Γ_{6} (r) = \frac{1}{N_{d i s k}} \sum_{c = 1}^{N_{c e l l}} γ_{6} (r_{c}) ϑ_{H} (R - | r - r_{c} |),

Figure 4

Download asset Open asset

Shape function.

On the left, we see a graphical representation of the sixfold shape function $γ_{6}$ (see Equation 2 for more details) for a generic irregular polygon. On the right (black six-legged star) the phase and magnitude of $γ_{6}$ for the same cell.

where $N_{d i s k} = \sum_{c} ϑ_{H} (R - | r - r_{c} |)$ is the number of cells whose centers lie within the disk, and the coarse-graining radius is fixed to be $R = 1.5 R_{c e l l}$ .

We choose to sample the shape function on a square grid with grid spacing equal to the nominal cell radius $R_{c e l l}$ . Topological defects are then identified by computing the winding number along each unit cell:

s = \frac{1}{2 π} \oint_{◻} d θ = \frac{1}{2 π} \sum_{n = 1}^{4} [θ (r_{n + 1}) - θ (r_{n})] m o d \frac{2 π}{6},

where the symbol $◻$ denotes a square unit cell in the interpolation grid and $θ = A r g (Γ_{6}) / 6$ the phase of the shape order parameter.

Nematic and hexatic stress in MPF simulations

The combined effect of confluency and non-equilibrium cellular migration leads to the onset of internal stresses in the tissue. These can be, in turn, classified based on their symmetry properties with respect to the cellular shape parameter. To characterize the properties of intercellular stresses in MPF simulations, we show in Figure 3a and b the longitudinal hexatic, nematic stress in the tissue. This is computed by contracting the nematic and hexatic concentration gradient tensor for each cell $c$ ( $𝝈_{c}^{(nem)} = - k_{φ} {(\nabla ϕ_{c})}^{\otimes 2}$ , $σ_{c}^{(h e x)} = - k_{φ} (\nabla φ_{c})^{\otimes 6}$ ) with the cellular shape tensor of the corresponding order ( $𝑸_{c, p} = 𝒏_{c, p}^{\otimes p}$ , with $p = 2, 6$ ). Here, $n_{c, p} = (\cos θ_{c, p}, \sin θ_{c, p})$ with $θ_{c, p} = A r g (γ_{c, p}) / p$ the angular orientation of the complex shape function $γ_{p}$ relative to the cth cell (see Figure 4). Therefore, the explicit expression of the longitudinal nematic and hexatic stress is given by

σ_{∥}^{(n e m)} = \sum_{c = 1}^{N_{c e l l}} σ_{c}^{(n e m)} : Q_{c, 2} = - k_{φ} \sum_{c = 1}^{N_{c e l l}} \sum_{i, j} (\partial_{i} φ_{c}) (\partial_{j} φ_{c}) (n_{2, c})_{i} (n_{2, c})_{j},

σ_{∥}^{(h e x)} = \sum_{c = 1}^{N_{c e l l}} σ_{c}^{(h e x)} : Q_{c, 6} = - k_{φ} \sum_{c = 1}^{N_{c e l l}} \sum_{i 1, \dots i_{6}} (\partial_{i_{1}} φ_{c}) \dots (\partial_{i_{6}} φ_{c}) (n_{6, c})_{i_{1}} \dots (n_{6, c})_{i_{6}} .

Notice that positive (negative) values of the longitudinal stress correspond to contractile (extensile) stresses for both nematic and hexatic contributions.

Correlating cellular mean square displacement and defect density

The scatter plot in Figure 3f, which correlates cellular mean square displacement with defect density, was constructed as follows. The system was divided into square subregions of size $Δ ℓ = 35$ , each containing approximately four cells. For each subregion, we analyzed a time window of $Δ t = 25 \times 10^{3}$ iterations, measuring both the normalized mean square displacemenof cells (relative to the subregion area $Δ ℓ^{2}$ ) and the average defect density. The normalized displacement is calculated as $m.s.d. = \sum_{t = t^{*}}^{t^{*} + Δ t} {| 𝒓_{c} (t) - 𝒓_{c} (t - 1) |}^{2} / Δ ℓ^{2}$ , where $t^{*}$ denotes the start time of the observation window. These measurements were collected for all subregions and time windows to generate the scatter plot.

The subregion size $Δ ℓ$ was selected to match the hexanematic crossover length scale, ensuring each area contains about four cells – the fundamental unit of tissue remodeling. This size balances the need to resolve defect-rich and defect-poor regions while avoiding excessive enlargement that would disguise local variations. Similarly, the observation window $Δ t$ was chosen to correspond to the typical time required for a cell to traverse a subregion of size $Δ ℓ$ . Shorter windows might miss remodeling events, while longer windows would average over different dynamical states.

Analysis of characteristic time of cell intercalation and T1 cycles

To analyze the temporal statistics of remodeling events, we start differentiating these between T1 cycles and cell intercalation. T1 cycles are identified as events where cells first change their neighbors at time $t$ and then they return to their initial configuration at a later time $t + T$ . Analogously, for cell intercalation events – leading to a permanent tissue remodeling – we measure the time interval $T$ between two consecutive intercalation events. The results of the analysis at varying the rotational diffusion $D_{r}$ are shown in Figure 3g. For each value of $D_{r}$ considered, we first identify cell intercalations and T1 cycles for each cell in the system, then we average the measured time intervals, obtaining the distribution of mean values, explicitly shown in the box plot in the inset of Figure 3g. The statistical analysis shows that the two populations (T1 cycles and intercalation events) are significantly different, with T1 cycles occurring faster than cell intercalation. We repeat this analysis at varying $D_{r}$ , and we plot in Figure 3g the characteristic time obtained as the mean value of the cell intercalations and T1 cycles distributions. Importantly, we find that, while the time interval between cell intercalation events sensibly depends on the rotational diffusion, the typical timescale of T1 cycles does not.

Active flow of a hexatic defect quadrupole

Scalar order parameter

In this section, we provide a derivation of Equation 1. To this end, let $Ψ_{6} = | Ψ_{6} | e^{6 i θ}$ be the hexatic complex order parameter and consider a quadrupole of $\pm 1 / 6$ disclinations equidistantly placed from the center of the primary cluster. The phase $θ = θ (r)$ is then given by the convolution of the average orientation in the surrounding of each defect, that is,

θ = - \frac{1}{6} \arctan (\frac{y}{x - ℓ}) - \frac{1}{6} \arctan (\frac{y}{x + ℓ}) + \frac{1}{6} \arctan (\frac{y + ℓ}{x}) + \frac{1}{6} \arctan (\frac{y - ℓ}{x}),

where $ℓ$ is distance from the center. The quadrupolar distance $ℓ$ is by definition taken to be small compared to the size of the system. Thus, expanding Equation 16 for $| 𝒓 | / ℓ ≫ 1$ , we obtain the simpler expression

θ = - \frac{2 ℓ^{2} \sin 2 ϕ}{3 | r |^{2}} + O (| ℓ / r |^{6}) .

Notice that the Taylor expansion features only the quadrupolar term of order $ℓ^{2}$ and is exact up to sixth order in $ℓ / | r |$ ; that is, the dipolar term, of order $𝒪 (ℓ / | 𝒓 |)$ , and all other terms up to the sixth order vanish identically.

This result is extremely robust and can be derived in a number of ways. For instance, Equation 16 can be obtained from the solution of the Poisson equation

\nabla^{2} φ = ρ_{d}

where φ is a dual field such that $\partial_{i} θ = - ϵ_{i j} \partial_{j} φ$ and the right-hand side of the Equation 18 is analogous to the electrostatic charge density (Chaikin et al., 1995). At large distance from the defects, Equation 18 can be solved by multipole expansion (Jackson, 1999), that is,

φ = a_{0} \log \frac{r_{0}}{| r |} + \sum_{n = 1}^{\infty} \frac{a_{n} \cos n θ + b_{n} \sin n θ}{| r |^{n}},

where $r_{0}$ is an irrelevant length scale and $a_{n}$ and $b_{n}$ are coefficients given by

a_{n} = \frac{1}{n} \int d A | r |^{n} \cos (n ϕ) ρ_{d},

b_{n} = \frac{1}{n} \int d A | r |^{n} \sin (n ϕ) ρ_{d} .

Thus, up to the quadrupole term, the expansion of φ is given by

φ = a_{0} \log \frac{r_{0}}{| r |} + \frac{a_{1} \cos ϕ + b_{1} \sin ϕ}{| r |} + \frac{a_{2} \cos 2 ϕ + b_{2} \sin 2 ϕ}{| r |^{2}} + \dots

As in electrostatics, the density $ρ_{d}$ is given by

ρ_{d} = \frac{1}{6} [- δ (r - ℓ e_{x}) - δ (r + ℓ e_{x}) + δ (r - ℓ e_{y}) + δ (r + ℓ e_{y})],

where $ℓ$ is again the distance from the center of the primary cluster. Now, because the defect quadrupole has, by construction, vanishing total strength and dipole moment, $a_{0} = 0$ and $a_{1} = b_{1} = 0$ . Of the quadrupolar terms, on the other hand, $a_{2} = - ℓ^{2} / 3$ and $b_{2} = 0$ , thus

φ = - \frac{ℓ^{2} \cos 2 ϕ}{3 | r |^{2}} .

Finally, going from φ to the original field θ one finds

θ = - \frac{2 ℓ^{2} \sin 2 ϕ}{3 | r |^{2}},

thus confirming the expression given in Equation 17.

Active force

To shed light on the structure of the cellular flow triggered by a T1 process, we solve the Stokes equation in the presence of an active force of the form $f^{(a)} = \nabla \cdot σ^{(a)}$ , where

σ^{(a)} = α_{6} \nabla^{\otimes 4} ⊙ Q_{6},

is the active hexatic stress tensor introduced in Armengol-Collado et al., 2023b. Calculating the divergence gives

\begin{array}{ll} f^{(a)} & = 960 \frac{α_{6} ℓ^{2}}{| r |^{7}} {[- 3 \cos 7 ϕ + \frac{ℓ^{2}}{| r |^{2}} (3 \cos 5 ϕ - 14 \cos 9 ϕ)] \\ e_{x} + [3 \sin 7 ϕ - \frac{ℓ^{2}}{| r |^{2}} (3 \sin 5 ϕ - 14 \sin 9 ϕ)] e_{y}}, \end{array}

up to correction of order $O (| ℓ / r |^{6})$ . A plot of the force field is shown in Figure 5a.

Figure 5

Download asset Open asset

Active hexatic defect quadrupole: convergent extension analytics.

(a) Force field: stream density plot of the force field Eq. (C11). It exhibits a clear, local, convergent-extension pattern in the vicinity of the quadrupolar radius $ℓ$ . (b) Velocity field: stream density plot of the velocity field Eq. (C17). It exhibits a clear, local, convergent-extension flow pattern in the vicinity of the quadrupolar radius $ℓ$ . (c) Velocity field approximated close to defect core: Stream density plot of the velocity field Equation 1. It exhibits a clear, local, convergent-extension flow pattern in the vicinity of the quadrupolar radius $ℓ$ . In all plots, the black disk corresponds to the radius of the quadrupole. Our analytical solution is valid outside the disk.

Flow field

To reconstruct the cellular motion generated by a T1 process, we solve the incompressible Stokes equation for the flow sourced the active force $𝒇^{(a)}$ : that is,

η \nabla^{2} v - \nabla P + f^{(a)} = 0,

\nabla \cdot 𝒗 = 0,

where η is the shear viscosity and $P$ the pressure. To this end, we turn to the Oseen formal solution

v (r) = \int_{0}^{2 π} d ϕ^{'} \int_{ℓ}^{R} d r^{'} r^{'} G (r - r^{'}) \cdot f^{(a)} (r^{'}),

where

G (r) = \frac{1}{4 π η} [(\log \frac{L}{| r |} - 1) 1 + \frac{r \otimes r}{| r |^{2}}],

is the two-dimensional Oseen tensor (see, e.g., Giomi, 2015), with $L$ a constant, and $R$ is a large distance cut-off. Without loss of generality, one can set $L = R \sqrt{e}$ in Equation 29. To calculate the integrals in Equation 28, we make use of the logarithmic expansion

\log \frac{| r - r^{'} |}{L} = \log \frac{r_{>}}{L} - \sum_{1}^{\infty} \frac{1}{m} {(\frac{r_{>}}{r_{<}})}^{m} \cos [m (ϕ - ϕ^{'})],

with $r_{≷}$ the maximum (minimum) between $| r |$ and $| 𝒓^{'} |$ , and of the orthogonality of trigonometric functions

\int_{0}^{2 π} d ϕ^{'} \cos [m (ϕ - ϕ^{'})] \cos n ϕ^{'} = π \cos n ϕ δ_{m n} .

The resulting flow field surrounding the defect quadrupole is then given by

\begin{array}{ll} \frac{v}{2 α_{6} ℓ^{2} / η} & = - [\frac{30 {(ℓ^{2} - {| r |}^{2})}^{2}}{{| r |}^{12}} (3 {| r |}^{2} \cos 8 ϕ + 14 ℓ^{2} \cos 10 ϕ) + \frac{60}{{| r |}^{8}} \cos 6 ϕ (3 {| r |}^{2} + 6 ℓ^{2} \log \frac{| r |}{ℓ} - 4 ℓ^{2})] r \\ + 6 (\frac{6}{{| r |}^{5}} - \frac{5 ℓ^{2}}{{| r |}^{7}}) (\cos 5 ϕ e_{x} - \sin 5 ϕ e_{y}) + \frac{30}{7} (\frac{6 ℓ^{2}}{{| r |}^{7}} - \frac{7}{{| r |}^{5}}) (\cos 7 ϕ e_{x} - \sin 7 ϕ e_{y}) \\ + \frac{35}{3} ℓ^{2} (\frac{8 ℓ^{2}}{{| r |}^{9}} - \frac{9}{{| r |}^{7}}) (\cos 9 ϕ e_{x} - \sin 9 ϕ e_{y}) . \end{array}

Figure 5b shows a plot of this flow, while Figure 5c shows a plot of the short distance approximation given in Equation 1 of the main text.

Numerical simulations of defect annihilation and unbinding

Numerical model and validation

The time-dependent flows in Figure 2ai and Figure 2bi as well as the orientational field color maps in Figure 2aii–iv and Figure 2bii–iv are obtained by numerically integrating Equation 3a using a vorticity-stream function finite difference scheme. All equations are discretized on a two-dimensional square grid of sizes 256 × 256 and 1024 × 1024 with periodic boundary conditions. For both grid sizes, the grid spacing is $Δ x = Δ y = 1$ and the time stepping $Δ t = 0.1$ . The validity of this numerical approach is benchmarked by many numerical studies on liquid crystals and active matter (see for example Krommydas et al., 2023; Giomi et al., 2022a; Giomi et al., 2022b; Giomi et al., 2014). In all simulations, we set $ρ = 1$ , $η = 1$ , $L_{6} = 0.5$ , $A_{6} = - 0.2$ , $B_{6} = 0.4$ , $Γ_{6} = 1$ and $λ_{6} = 1.11$ . All parameters are expressed in the arbitrary units used in the numerical simulations.

Defect annihilation and unbinding

To construct the initial configuration of $Ψ_{6}$ , we set $ℓ = 7$ , $| Ψ_{6} | = 1$ and take θ as given in Equation 16 inside a disk of radius $R_{D} = 28$ and random outside. We then thermalize this configuration by keeping the orientation of the order parameter in the disk fixed and relaxing $Ψ_{6}$ everywhere else. This allows us to obtain a defect-free configuration where $| Ψ_{6} | \approx 1$ everywhere, except that close to the defect cores where $| Ψ_{6} | \approx 0$ . Notice that, on a doubly periodical domain, $\sum_{i} s_{i} = 0$ . Therefore, no other topological defect is found at the end of such relaxation procedure. Simulations are carried out until the total free energy relative variation drops under 0.1% with respect to two consecutive iterations. This corresponds to a state where defects have annihilated, and the hexatic liquid crystal has achieved a smooth configuration everywhere in the simulation box. For both annihilation and unbinding numerical experiments, we scan $α_{6}$ for a wide range of positive and negative values. For any negative values of activity, we obtain increasingly sheared versions of the flow pattern in Figure 2ai. Similarly, for positive values of $α_{6}$ we obtain increasingly sheared versions of the same flow pattern, but with the direction of the flow inverted.

Active defect dipole annihilation: The origin of the unbinding

In this section, we provide a brief account of the annihilation dynamics of a pair of $\pm 1 / 6$ active hexatic defects, in which it is possible to recognize the fundamental mechanism driving defect unbinding. To this end, we place the defects on the x-axis at a distance of $Δ x = 64$ and construct the initial configuration of the hexatic order parameter $Ψ_{6} = e^{6 i θ}$ by setting $θ = \pm \arctan [y / (x \pm Δ x / 2)]$ inside a disk of radius $R_{D} = 5$ centered at the defect cores, and random outside. We thermalize this configuration by keeping the phase of the order parameter in the two disks fixed, while relaxing $Ψ_{6}$ everywhere else. As before, this procedure allows us to obtain a state where $| Ψ_{6} | \approx 1$ everywhere, except that close to the two defect cores where $| Ψ_{6} | \approx 0$ . We use this as the initial state for our annihilation experiment. Simulations are carried out until defects have annihilated and the total free energy relative variation drops under 0.1% with respect to two consecutive iterations. The model parameters, expressed in lattice units, are again: $Δ t = 1$ , $ρ = 1$ , $η = 1$ , $L_{6} = 0.5$ , $A_{6} = - 0.2$ , $B_{6} = 0.4$ , $Γ_{6} = 1$ , and $λ_{6} = 1.11$ .

Figure 6 shows the trajectories of the positive (red) and negative (blue) defects during annihilation, for three realizations of the activity parameter $α_{6}$ , that is, $α_{6} = 0.1$ (contractile), $α_{6} = 0$ (passive), and $α_{6} = - 0.1$ (extensile). For contractile activity, the backflow sourced by the active stress, Equation 25, annihilation is sped up with respect to the passive case. By contrast, for extensile activity, annihilation is delayed. The same effects lead to the breakup of the quadrupole into two defect pairs, provided the repulsive forces introduced by the active flow overcome the attractive Coulomb-like forces between defects.

Figure 6

Download asset Open asset

Trajectories of annihilating defects in time The red lines are the trajectories of positive and blue negative defects, respectively.

Defects are sped up by positive activity ( $α_{6} = 0.1$ ), but they are slowed down instead by negative activity ( $α_{6} = - 0.1$ ).

Figure 6—source data 1 All positions of topological defects for all iterations needed to reproduce the defect annihilation trajectory plots in this figure.: https://cdn.elifesciences.org/articles/105397/elife-105397-fig6-data1-v1.zip
Download elife-105397-fig6-data1-v1.zip
Figure 6—source data 2 All positions of topological defects for all iterations needed to reproduce the defect annihilation trajectory plots in this figure.: https://cdn.elifesciences.org/articles/105397/elife-105397-fig6-data2-v1.zip
Download elife-105397-fig6-data2-v1.zip
Figure 6—source data 3 All positions of topological defects for all iterations needed to reproduce the defect annihilation trajectory plots in this figure.: https://cdn.elifesciences.org/articles/105397/elife-105397-fig6-data3-v1.zip
Download elife-105397-fig6-data3-v1.zip
Figure 6—source data 4 All positions of topological defects for all iterations needed to reproduce the defect annihilation trajectory plots in this figure.: https://cdn.elifesciences.org/articles/105397/elife-105397-fig6-data4-v1.zip
Download elife-105397-fig6-data4-v1.zip
Figure 6—source data 5 All positions of topological defects for all iterations needed to reproduce the defect annihilation trajectory plots in this figure.: https://cdn.elifesciences.org/articles/105397/elife-105397-fig6-data5-v1.zip
Download elife-105397-fig6-data5-v1.zip
Figure 6—source data 6 All positions of topological defects for all iterations needed to reproduce the defect annihilation trajectory plots in this figure.: https://cdn.elifesciences.org/articles/105397/elife-105397-fig6-data6-v1.zip
Download elife-105397-fig6-data6-v1.zip
Figure 6—source data 7 All positions of topological defects for all iterations needed to reproduce the defect annihilation trajectory plots in this figure.: https://cdn.elifesciences.org/articles/105397/elife-105397-fig6-data7-v1.zip
Download elife-105397-fig6-data7-v1.zip
Figure 6—source data 8 All positions of topological defects for all iterations needed to reproduce the defect annihilation trajectory plots in this figure.: https://cdn.elifesciences.org/articles/105397/elife-105397-fig6-data8-v1.zip
Download elife-105397-fig6-data8-v1.zip
Figure 6—source data 9 All positions of topological defects for all iterations needed to reproduce the defect annihilation trajectory plots in this figure.: https://cdn.elifesciences.org/articles/105397/elife-105397-fig6-data9-v1.zip
Download elife-105397-fig6-data9-v1.zip
Figure 6—source data 10 All positions of topological defects for all iterations needed to reproduce the defect annihilation trajectory plots in this figure.: https://cdn.elifesciences.org/articles/105397/elife-105397-fig6-data10-v1.zip
Download elife-105397-fig6-data10-v1.zip
Figure 6—source data 11 All positions of topological defects for all iterations needed to reproduce the defect annihilation trajectory plots in this figure.: https://cdn.elifesciences.org/articles/105397/elife-105397-fig6-data11-v1.zip
Download elife-105397-fig6-data11-v1.zip
Figure 6—source data 12 All positions of topological defects for all iterations needed to reproduce the defect annihilation trajectory plots in this figure.: https://cdn.elifesciences.org/articles/105397/elife-105397-fig6-data12-v1.zip
Download elife-105397-fig6-data12-v1.zip
Figure 6—source data 13 All positions of topological defects for all iterations needed to reproduce the defect annihilation trajectory plots in this figure.: https://cdn.elifesciences.org/articles/105397/elife-105397-fig6-data13-v1.zip
Download elife-105397-fig6-data13-v1.zip
Figure 6—source data 14 All positions of topological defects for all iterations needed to reproduce the defect annihilation trajectory plots in this figure.: https://cdn.elifesciences.org/articles/105397/elife-105397-fig6-data14-v1.zip
Download elife-105397-fig6-data14-v1.zip
Figure 6—source data 15 All positions of topological defects for all iterations needed to reproduce the defect annihilation trajectory plots in this figure.: https://cdn.elifesciences.org/articles/105397/elife-105397-fig6-data15-v1.zip
Download elife-105397-fig6-data15-v1.zip
Figure 6—source data 16 All positions of topological defects for all iterations needed to reproduce the defect annihilation trajectory plots in this figure.: https://cdn.elifesciences.org/articles/105397/elife-105397-fig6-data16-v1.zip
Download elife-105397-fig6-data16-v1.zip
Figure 6—source data 17 All positions of topological defects for all iterations needed to reproduce the defect annihilation trajectory plots in this figure.: https://cdn.elifesciences.org/articles/105397/elife-105397-fig6-data17-v1.zip
Download elife-105397-fig6-data17-v1.zip
Figure 6—source data 18 All positions of topological defects for all iterations needed to reproduce the defect annihilation trajectory plots in this figure.: https://cdn.elifesciences.org/articles/105397/elife-105397-fig6-data18-v1.zip
Download elife-105397-fig6-data18-v1.zip

Data availability

The current manuscript is a theoretical study. Figure 2—source data 1, Figure 3—source data 1 and Figure 6—source data 1–18 contain the numerical data used to generate the Figures 2, 3 and 6. The source code used to generate numerical data is enclosed as Source code 1.

References

1. Alert R
2. Trepat X
(2020) Physical models of collective cell migration
Annual Review of Condensed Matter Physics 11:77–101.

https://doi.org/10.1146/annurev-conmatphys-031218-013516
- Google Scholar
(2023a) Epithelia are multiscale active liquid crystals
Nature Physics 19:1773–1779.

https://doi.org/10.1038/s41567-023-02179-0
- Google Scholar
(2023b) Hydrodynamics and multiscale order in confluent epithelia
eLife 13:e86400.

https://doi.org/10.7554/eLife.86400
- PubMed
- Google Scholar
1. Balasubramaniam L
2. Doostmohammadi A
3. Saw TB
4. Narayana G
5. Mueller R
6. Dang T
7. Thomas M
8. Gupta S
9. Sonam S
10. Yap AS
11. Toyama Y
12. Mège RM
13. Yeomans JM
14. Ladoux B
(2021) Investigating the nature of active forces in tissues reveals how contractile cells can form extensile monolayers
Nature Materials 20:1156–1166.

https://doi.org/10.1038/s41563-021-00919-2
- PubMed
- Google Scholar
(2004) Myosin-dependent junction remodelling controls planar cell intercalation and axis elongation
Nature 429:667–671.

https://doi.org/10.1038/nature02590
- PubMed
- Google Scholar
1. Blanchard GB
(2017) Taking the strain: quantifying the contributions of all cell behaviours to changes in epithelial shape
Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences 372:20150513.

https://doi.org/10.1098/rstb.2015.0513
- PubMed
- Google Scholar
(2018) Turbulent dynamics of epithelial cell cultures
Physical Review Letters 120:208101.

https://doi.org/10.1103/PhysRevLett.120.208101
- PubMed
- Google Scholar
1. Buckley CD
2. Tan J
3. Anderson KL
4. Hanein D
5. Volkmann N
6. Weis WI
7. Nelson WJ
8. Dunn AR
(2014) Cell adhesion the minimal cadherin-catenin complex binds to actin filaments under force
Science 346:1254211.

https://doi.org/10.1126/science.1254211
- PubMed
- Google Scholar
1. Camley BA
2. Zhang Y
3. Zhao Y
4. Li B
5. Ben-Jacob E
6. Levine H
7. Rappel WJ
(2014) Polarity mechanisms such as contact inhibition of locomotion regulate persistent rotational motion of mammalian cells on micropatterns
PNAS 111:14770–14775.

https://doi.org/10.1073/pnas.1414498111
- PubMed
- Google Scholar
(2019) Lattice Boltzmann methods and active fluids
The European Physical Journal. E, Soft Matter 42:81.

https://doi.org/10.1140/epje/i2019-11843-6
- PubMed
- Google Scholar
1. Cavey M
2. Rauzi M
3. Lenne P-F
4. Lecuit T
(2008) A two-tiered mechanism for stabilization and immobilization of E-cadherin
Nature 453:751–756.

https://doi.org/10.1038/nature06953
- PubMed
- Google Scholar
Book
(1995)
Principles of Condensed Matter Physics

Cambridge University Press.
- Google Scholar
1. Cheung KJ
2. Ewald AJ
(2016) A collective route to metastasis: Seeding by tumor cell clusters
Science 352:167–169.

https://doi.org/10.1126/science.aaf6546
- PubMed
- Google Scholar
1. Cislo DJ
2. Yang F
3. Qin H
4. Pavlopoulos A
5. Bowick MJ
6. Streichan SJ
(2023) Active cell divisions generate fourfold orientationally ordered phase in living tissue
Nature Physics 19:1201–1210.

https://doi.org/10.1038/s41567-023-02025-3
- PubMed
- Google Scholar
1. Das A
2. Sastry S
3. Bi D
(2021) Controlled neighbor exchanges drive glassy behavior, intermittency, and cell streaming in epithelial tissues
Physical Review X 11:041037.

https://doi.org/10.1103/PhysRevX.11.041037
- Google Scholar
(2017) Topological defects in confined populations of spindle-shaped cells
Nature Physics 13:58–62.

https://doi.org/10.1038/nphys3876
- Google Scholar
1. Duclos G
2. Blanch-Mercader C
3. Yashunsky V
4. Salbreux G
5. Joanny JF
6. Prost J
7. Pd S
(2018) Spontaneous shear flow in confined cellular nematics
Nature Physics 14:728–732.

https://doi.org/10.1038/s41567-018-0099-7
- Google Scholar
(2022) Active T1 transitions in cellular networks
The European Physical Journal. E, Soft Matter 45:29.

https://doi.org/10.1140/epje/s10189-022-00175-5
- PubMed
- Google Scholar
1. Durand M
2. Heu J
(2019) Thermally driven order-disorder transition in two-dimensional soft cellular systems
Physical Review Letters 123:188001.

https://doi.org/10.1103/PhysRevLett.123.188001
- PubMed
- Google Scholar
1. Eckert J
2. Ladoux B
3. Mège R-M
4. Giomi L
5. Schmidt T
(2023) Hexanematic crossover in epithelial monolayers depends on cell adhesion and cell density
Nature Communications 14:5762.

https://doi.org/10.1038/s41467-023-41449-6
- PubMed
- Google Scholar
1. Engl W
2. Arasi B
3. Yap LL
4. Thiery JP
5. Viasnoff V
(2014) Actin dynamics modulate mechanosensitive immobilization of E-cadherin at adherens junctions
Nature Cell Biology 16:587–594.

https://doi.org/10.1038/ncb2973
- PubMed
- Google Scholar
1. Erdemci-Tandogan G
2. Manning ML
(2021) Effect of cellular rearrangement time delays on the rheology of vertex models for confluent tissues
PLOS Computational Biology 17:e1009049.

https://doi.org/10.1371/journal.pcbi.1009049
- PubMed
- Google Scholar
(2014) Vertex models of epithelial morphogenesis
Biophysical Journal 106:2291–2304.

https://doi.org/10.1016/j.bpj.2013.11.4498
- PubMed
- Google Scholar
1. Friedl P
2. Gilmour D
(2009) Collective cell migration in morphogenesis, regeneration and cancer
Nature Reviews. Molecular Cell Biology 10:445–457.

https://doi.org/10.1038/nrm2720
- PubMed
- Google Scholar
Book
1. Gennes PGD
2. Prost J
(1993) The Physics of Liquid Crystals. International Series of Monographs on Physics
Clarendon Press.

https://doi.org/10.1093/oso/9780198520245.001.0001
- Google Scholar
(2014) Defect dynamics in active nematics
Philosophical Transactions 372:20130365.

https://doi.org/10.1098/rsta.2013.0365
- PubMed
- Google Scholar
1. Giomi L
(2015) Geometry and topology of turbulence in active nematics
Physical Review X 5:031003.

https://doi.org/10.1103/PhysRevX.5.031003
- Google Scholar
1. Giomi L
2. Toner J
3. Sarkar N
(2022a) Hydrodynamic theory of p-atic liquid crystals
Physical Review. E 106:024701.

https://doi.org/10.1103/PhysRevE.106.024701
- PubMed
- Google Scholar
1. Giomi L
2. Toner J
3. Sarkar N
(2022b) Long-ranged order and flow alignment in sheared p-atic liquid crystals
Physical Review Letters 129:067801.

https://doi.org/10.1103/PhysRevLett.129.067801
- PubMed
- Google Scholar
(2008) Discrete rearranging disordered patterns, part I: robust statistical tools in two or three dimensions
The European Physical Journal. E, Soft Matter 25:349–369.

https://doi.org/10.1140/epje/i2007-10298-8
- PubMed
- Google Scholar
(2017) Taming active turbulence with patterned soft interfaces
Nature Communications 8:564.

https://doi.org/10.1038/s41467-017-00617-1
- PubMed
- Google Scholar
1. Haeger A
2. Alexander S
3. Vullings M
4. Kaiser FMP
5. Veelken C
6. Flucke U
7. Koehl GE
8. Hirschberg M
9. Flentje M
10. Hoffman RM
11. Geissler EK
12. Kissler S
13. Friedl P
(2020) Collective cancer invasion forms an integrin-dependent radioresistant niche
Journal of Experimental Medicine 217:e20181184.

https://doi.org/10.1084/jem.20181184
- Google Scholar
(2020) Chiral stresses in nematic cell monolayers
Soft Matter 16:764–774.

https://doi.org/10.1039/c9sm01851d
- PubMed
- Google Scholar
1. Honda H
2. Eguchi G
(1980) How much does the cell boundary contract in a monolayered cell sheet?
Journal of Theoretical Biology 84:575–588.

https://doi.org/10.1016/s0022-5193(80)80021-x
- PubMed
- Google Scholar
(2023) Mechanochemical active feedback generates convergence extension in epithelial tissue
Physical Review Letters 131:238301.

https://doi.org/10.1103/PhysRevLett.131.238301
- PubMed
- Google Scholar
1. Irvine KD
2. Wieschaus E
(1994) Cell intercalation during Drosophila germband extension and its regulation by pair-rule segmentation genes
Development 120:827–841.

https://doi.org/10.1242/dev.120.4.827
- PubMed
- Google Scholar
Book
1. Jackson JD
(1999)
Classical Electrodynamics

Wiley.
- Google Scholar
(2023) Robust statistical properties of T1 transitions in a multi-phase field model of cell monolayers
Scientific Reports 13:10096.

https://doi.org/10.1038/s41598-023-37064-6
- PubMed
- Google Scholar
(2017) Topological defects control collective dynamics in neural progenitor cell cultures
Nature 545:327–331.

https://doi.org/10.1038/nature22321
- PubMed
- Google Scholar
1. Keller R
2. Davidson L
3. Edlund A
4. Elul T
5. Ezin M
6. Shook D
7. Skoglund P
(2000) Mechanisms of convergence and extension by cell intercalation
Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences 355:897–922.

https://doi.org/10.1098/rstb.2000.0626
- PubMed
- Google Scholar
(2008) The forces that shape embryos: physical aspects of convergent extension by cell intercalation
Physical Biology 5:015007.

https://doi.org/10.1088/1478-3975/5/1/015007
- PubMed
- Google Scholar
1. Kosterlitz JM
2. Thouless DJ
(1972) Long range order and metastability in two dimensional solids and superfluids
Journal of Physics C: Solid State Physics 5:L124.

https://doi.org/10.1088/0022-3719/5/11/002
- Google Scholar
1. Kosterlitz JM
2. Thouless DJ
(1973) Long range order and metastability in two dimensional solids and superfluids
Journal of Physics C: Solid State Physics 6:1181.

https://doi.org/10.1088/0022-3719/6/7/010
- Google Scholar
1. Kosterlitz JM
(2016) Kosterlitz-Thouless physics: a review of key issues
Reports on Progress in Physics 79:026001.

https://doi.org/10.1088/0034-4885/79/2/026001
- PubMed
- Google Scholar
(2023) Hydrodynamic enhancement of p-atic defect dynamics
Physical Review Letters 130:098101.

https://doi.org/10.1103/PhysRevLett.130.098101
- PubMed
- Google Scholar
1. Lemma LM
2. DeCamp SJ
3. You Z
4. Giomi L
5. Dogic Z
(2019) Statistical properties of autonomous flows in 2D active nematics
Soft Matter 15:3264–3272.

https://doi.org/10.1039/c8sm01877d
- PubMed
- Google Scholar
1. Li YW
2. Ciamarra MP
(2018) Role of cell deformability in the two-dimensional melting of biological tissues
Physical Review Materials 2:045602.

https://doi.org/10.1103/PhysRevMaterials.2.045602
- Google Scholar
(2020) Solid-liquid transition of deformable and overlapping active particles
Physical Review Letters 125:038003.

https://doi.org/10.1103/PhysRevLett.125.038003
- PubMed
- Google Scholar
(2008) Discrete rearranging disordered patterns, part II: 2D plasticity, elasticity and flow of a foam
The European Physical Journal. E, Soft Matter 25:371–384.

https://doi.org/10.1140/epje/i2007-10300-7
- PubMed
- Google Scholar
(2023) Mechanical basis and topological routes to cell elimination
eLife 12:e82435.

https://doi.org/10.7554/eLife.82435
- PubMed
- Google Scholar
(2024) Biophysical control of plasticity and patterning in regeneration and cancer
Cellular and Molecular Life Sciences 81:9.

https://doi.org/10.1007/s00018-023-05054-6
- Google Scholar
1. Nelson DR
2. Halperin BI
(1979) Dislocation-mediated melting in two dimensions
Physical Review B 19:2457–2484.

https://doi.org/10.1103/PhysRevB.19.2457
- Google Scholar
1. Palmieri B
2. Bresler Y
3. Wirtz D
4. Grant M
(2015) Multiple scale model for cell migration in monolayers: Elastic mismatch between cells enhances motility
Scientific Reports 5:11745.

https://doi.org/10.1038/srep11745
- PubMed
- Google Scholar
1. Paré AC
2. Zallen JA
(2020) Cellular, molecular, and biophysical control of epithelial cell intercalation
Current Topics in Developmental Biology 136:167–193.

https://doi.org/10.1016/bs.ctdb.2019.11.014
- PubMed
- Google Scholar
(2020) Hexatic phase in a model of active biological tissues
Soft Matter 16:3914–3920.

https://doi.org/10.1039/d0sm00109k
- PubMed
- Google Scholar
(2019) Sustained oscillations of epithelial cell sheets
Biophysical Journal 117:464–478.

https://doi.org/10.1016/j.bpj.2019.06.013
- PubMed
- Google Scholar
1. Rauzi M
(2020) Cell intercalation in a simple epithelium
Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences 375:20190552.

https://doi.org/10.1098/rstb.2019.0552
- PubMed
- Google Scholar
(2018) Radially patterned cell behaviours during tube budding from an epithelium
eLife 7:e35717.

https://doi.org/10.7554/eLife.35717
- PubMed
- Google Scholar
1. Saw TB
2. Doostmohammadi A
3. Nier V
4. Kocgozlu L
5. Thampi S
6. Toyama Y
7. Marcq P
8. Lim CT
9. Yeomans JM
10. Ladoux B
(2017) Topological defects in epithelia govern cell death and extrusion
Nature 544:212–216.

https://doi.org/10.1038/nature21718
- PubMed
- Google Scholar
1. Schwarz US
2. Safran SA
(2002) Elastic interactions of cells
Physical Review Letters 88:048102.

https://doi.org/10.1103/PhysRevLett.88.048102
- PubMed
- Google Scholar
1. Serra-Picamal X
2. Conte V
3. Vincent R
4. Anon E
5. Tambe DT
6. Bazellieres E
7. Butler JP
8. Fredberg JJ
9. Trepat X
(2012) Mechanical waves during tissue expansion
Nature Physics 8:628–634.

https://doi.org/10.1038/nphys2355
- Google Scholar
(2023) Generating active T1 transitions through mechanochemical feedback
eLife 12:e79862.

https://doi.org/10.7554/eLife.79862
- PubMed
- Google Scholar
1. Staple DB
2. Farhadifar R
3. Röper JC
4. Aigouy B
5. Eaton S
6. Jülicher F
(2010) Mechanics and remodelling of cell packings in epithelia
The European Physical Journal. E, Soft Matter 33:117–127.

https://doi.org/10.1140/epje/i2010-10677-0
- PubMed
- Google Scholar
(2016) Unipolar distributions of junctional Myosin II identify cell stripe boundaries that drive cell intercalation throughout Drosophila axis extension
eLife 5:e12094.

https://doi.org/10.7554/eLife.12094
- PubMed
- Google Scholar
1. Tetley RJ
2. Mao Y
(2018) The same but different: cell intercalation as a driver of tissue deformation and fluidity
Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences 373:20170328.

https://doi.org/10.1098/rstb.2017.0328
- PubMed
- Google Scholar
1. Walck-Shannon E
2. Hardin J
(2014) Cell intercalation from top to bottom
Nature Reviews. Molecular Cell Biology 15:34–48.

https://doi.org/10.1038/nrm3723
- PubMed
- Google Scholar
(2020) Anisotropy links cell shapes to tissue flow during convergent extension
PNAS 117:13541–13551.

https://doi.org/10.1073/pnas.1916418117
- PubMed
- Google Scholar
Book
1. Weaire DL
2. Hutzler S
(1999) The Physics of Foams
Oxford, UK: Oxford University Press.

https://doi.org/10.1093/oso/9780198505518.001.0001
- Google Scholar
1. Yamada S
2. Pokutta S
3. Drees F
4. Weis WI
5. Nelson WJ
(2005) Deconstructing the cadherin-catenin-actin complex
Cell 123:889–901.

https://doi.org/10.1016/j.cell.2005.09.020
- PubMed
- Google Scholar
(2022) Chiral edge current in nematic cell monolayers
Physical Review X 12:041017.

https://doi.org/10.1103/PhysRevX.12.041017
- Google Scholar
(2010) alpha-Catenin as a tension transducer that induces adherens junction development
Nature Cell Biology 12:533–542.

https://doi.org/10.1038/ncb2055
- PubMed
- Google Scholar
1. You Z
2. Pearce DJG
3. Sengupta A
4. Giomi L
(2018) Geometry and mechanics of microdomains in growing bacterial colonies
Physical Review X 8:031065.

https://doi.org/10.1103/PhysRevX.8.031065
- Google Scholar
1. Young AP
(1979) Melting and the vector Coulomb gas in two dimensions
Physical Review B 19:1855–1866.

https://doi.org/10.1103/PhysRevB.19.1855
- Google Scholar
1. Zallen JA
2. Wieschaus E
(2004) Patterned gene expression directs bipolar planar polarity in Drosophila
Developmental Cell 6:343–355.

https://doi.org/10.1016/s1534-5807(04)00060-7
- PubMed
- Google Scholar

Article and author information

Author details

Dimitrios Krommydas

Instituut-Lorentz, Universiteit Leiden, Leiden, Netherlands

Contribution
Conceptualization, Data curation, Formal analysis, Investigation, Writing – original draft, Writing – review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0001-5936-0949
Livio N Carenza
1. Instituut-Lorentz, Universiteit Leiden, Leiden, Netherlands
2. Physics Department, College of Sciences, Koç University, Istanbul, Turkiye
Contribution
Conceptualization, Data curation, Formal analysis, Investigation, Writing – original draft, Writing – review and editing

Competing interests
No competing interests declared
Luca Giomi

Instituut-Lorentz, Universiteit Leiden, Leiden, Netherlands

Contribution
Conceptualization, Formal analysis, Supervision, Funding acquisition, Investigation, Writing – original draft, Writing – review and editing

For correspondence
giomi@lorentz.leidenuniv.nl

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0001-7740-5960

Funding

European Research Council (HexaTissue)

Dimitrios Krommydas

European Molecular Biology Organization (ALTF 353-2023)

Livio N Carenza

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

This work was supported by the European Union via the ERC-CoGgrant HexaTissue and by Netherlands Organization for Scientific Research (NWO/OCW). LNC acknowledges the support of the Postdoctoral EMBO Fellowship ALTF 353-2023. Part of this work was carried out on the Dutch national e-infrastructure with the support of SURF through the Grant 2021.028 for computational time. The computer codes used to produce the numerical data presented in the article are available in the manuscript.

Version history

Sent for peer review: December 24, 2024
Preprint posted: December 27, 2024
Reviewed Preprint version 1: May 13, 2025
Reviewed Preprint version 2: November 20, 2025
Version of Record published: May 5, 2026

Cite all versions

You can cite all versions using the DOI https://doi.org/10.7554/eLife.105397. This DOI represents all versions, and will always resolve to the latest one.

Copyright

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.