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, orT1 for brevity, through which the vertices ofa honeycomb network merge and then split, thereby leading to a remodelling 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 cell (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 of 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 approximatively two and half time its initial length. In this process, individual cells persistently move across the anterior-posterior axis more than ten times their body length. In Drosophila, intercalation is also the main process driving the salivary gland tube formation, in which cells radially converge towards 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 with various stages of the so-called metastatic cascade: i.e. 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; 3) colonization ofa 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., 2019; 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 on recent advances toward deciphering orientational order in epithelial layers (Li and Ciamarra, 2018; Pasu-palak et al., 2020; Durand and Heu, 2019; Armengol-Collado et al., 2023a,b; 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,b; 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: i.e. point-like singularities in the otherwise regular configuration of a continuous p-atic order parameter - i.e. Ψ_p = 〈e^ipϑ〉, with ϑ the orientation of the individual building blocks and (〈…〉) the ensemble average (Giomi et al., 2022b,a) - around which the average cellular orientation rotates by 2πs, with s = ±1/p, ±2/p … 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, 1973; Nelson and Halperin, 1979; Young, 1979; Kosterlitz, 2016). According to this, the hierarchical unbinding of neutral defect complexes - i.e. for which - 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 Fig. 1a for a cluster of four cells, hereafter referred to as primary cell cluster. AT1 occurs when the internal 3-fold coordinated vertices shrink until merging into a 4-fold vertex, and then split 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., Refs. (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.

A full cell intercalation then consists of a burst of internal and a peripheral T1 processes, as is schematically summarized in Fig. 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 3-fold vertices shrink to a 4-fold vertex, the cluster may reverse its dynamics and return to its original configuration (see Fig. 1b). In the following, we refer to this scenario as T1 cycle: i.e. 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.

Figure 1

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 Fig. 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 ±1/6 disclinations. As in two-dimensional melting, such a defect-structure can arise spontaneously, as consequence of spatiotemporal fluctuations of physical or biological nature. The two complementary remodelling events following the intermediate configuration, i.e. cell intercalation (Fig. 1a) and the T1 cycle (Fig. 1b), correspond instead to the two possible “fates” of this initial excitation. Upon unbinding [Figs. 1a-(i) and 1b-(i)], the four defects comprising the quadrupole can further unbind into two ±1/6 disclination dipoles [Fig. 1a-(ii)], or annihilate [Fig. 1b-(ii)], 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 n/6 with respect to initial configuration. 2) Both the primary cluster and its surrounding cells, must perform a local convergent extension: i.e. 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 Ref. (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, Ψ₆|Ψ₆|e^6iθ, is θ = 0, and construct a configuration featuring a quadrupole of ±1/6 disclinations [see Figs. 2a-(i) and 2b-(i)]. Along one full loop encircling each of these defects, θ changes by ±π/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 ≤ θ ≥ π/6 around the quadruple, as indicated by the alternating blue and red tones in Fig. 2a-(ii) and 2b-(ii). 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 Fig. 2a-(iii)], while leaving the orientation of the cells in the exterior essentially undistorted. Thus, the unbinding ofa ±1/6 defect quadrupole from a defect-free configuration and its break up into two dipoles drives a π /6 rotation of the cells between the dipoles [see Fig. 2a-(iv)], consistently with our first requirement. Conversely, if defects annihilate [see Fig. 2b-(iii)], the cells’ initial orientation is restored after a transient orientational fluctuation [see Fig. 2b-(iv)].

In order to address the second and third requirements, we look at the configuration of the velocity field υ, 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, this is given by

Figure 2

where α₆ a constant, with 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 the primary cluster, where |r|/ℓ≧ 1, Eq. (1) gives υ_x(x, 0) ~ α₆/(ηℓ⁴)x and υ_y(y, 0) ~ −α₆/(ηℓ⁴)y. In agreement with our second requirement, this is 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, Eq. (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 - i.e. the sign of υ-depends solely on α₆, 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 remodelling 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 Figs. 2a and 2b, for two different α₆ values. When α₆ is large and negative, the quadrupole splits into two ±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 (de Gennes and Prost, 1993). This biases the unbinding dynamics of the defect dipoles, thereby setting, in concert with the passive Coulomblike forces at play, the direction along which the dipoles move away from each other (Krommydas et al., 2023). Finally, switching off active stresses - i.e. α₆ = 0, shown in Fig. 2b - suppresses both defect unbinding and shear. Consistently, inverting the direction of active forces from extensile to contractile - i.e. α₆ > 0, see Methods - results in a speed up of defects annihilation, hence of the T1 cycle, but never lead to the unbinding of ±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 Eq. (1) represents a central outcome, is further corroborated by recent numerical work by Erdemci-Tandogan et al. (Erdemci-Tandogan and Manning, 2021) and Das et al. (Das et al., 2021), who, using a different cell-resolved model of epithelia - i.e. the Vertex model (Honda and Eguchi, 1980) - showed that the rate ofT1 processes enhances the fluidity of tissues. Consistently, the speed of the stagnation flow sourced by cell intercalation increases like η⁻¹, 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,b; 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 υ₀, equal for all cells, and a fluctuating direction of motion, with rotational diffusion coefficient D_r. Adetailed 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 υ₀ and the total number of cells and vary the persistence length of the cells trajectories by tuning the rotational diffusion coefficient D_r. Fig. 3a,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).

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

where r_υ, with υ = 1, 2… V is the positions of the υ-th vertex with respect to the centroid of the polygon and ϕ_υ = arctan(y_υ/x_υ) its orientation with respect to the x-axis. As demonstrated in Ref. (Armengol-Collado et al., 2023b), the phase Arg(γ₆)/6 of the shape function identifies the 6-fold orientation of a cell and is represented as white six-legged stars in Fig. 3d. The shape function is then coarse-grained over the length scale R = 1.5R_cell to reconstruct the shape parameter Γ₆ = (〈γ₆〉_R. Finally, topological defects are identified as singularities in the otherwise smoothly varying field Γ₆ = Γ₆(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 Figs. 3c).

For each D_r value, we reconstruct the probability of finding a T1 at a given distance from a defect (see Figs. 3c and 3d) 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 Fig. 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 Fig. 3f, the mean squared displacement of T1 cells versus the local defect density. Our data confirm the existence of two different classes: i.e. “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 to 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 and order of magnitude larger than that of cells affected by an isolated T1 cycle (Fig. 3f). Hence, although isolated internal T1’s (and by extension T1 cycles) can have small long-ranged effects, those effects are negligible in comparison the collective long-ranged motion induced by a full cell intercalation.

Finally, following Refs. (Das et al., 2021; Erdemci-Tandogan and Manning, 2021) we study the temporal correlations ofT1 cycles and cell intercalation. Fig. 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 is a monotonically decreasing function of 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 Fig. 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 on 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 Eq. (1) and neglecting irrelevant inertial effects, one can estimate the magnitude of the active forces at play: i.e. F_active ~ α₆/ℓ³, where ℓ denotes the range of the distortion caused by the unbinding dislocations. Such a scaling form results primarily from the 6-fold 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., 2016; Saw et al., 2017; Kawaguchi et al., 2017; Blanch-Mercader et al., 2018; Yashunsky et al., 2022) - that a6 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_passive ~ L₆/ℓ, with L₆ the orientational stiffness of the hexatic phase (see Methods). Equating F_active and F_passive provides an estimate of the typical size of the primary cluster at which the cellular layer becomes unstable to collective migration: i.e. . This active hexatic length scale was identified in Ref. (Armengol-Collado et al., 2023b) and is believed to play a role analog to that ofthe active nematic length scale (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 ℓ₆ has the effect of suppressing cell intercalation and possibly render necessary different locomotion strategies and possibly favoring a switch to mesenchymal phenotypes. Lastly, in agreement with 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: i.e. only extensile forces can provide the type of repulsive interactions that are necessary for defects to unbind, thus to trigger cell intercalation.

Figure 3

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 modelling. 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 ±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. 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 ofin 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 Ref. (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 υ, this can be described in terms of the 6-fold order parameter tensor Q₆ = 4〚〈v^⊗6〉〛 = 4|Ψ₆| 〚n^⊗6〛, where v = cos ϑe_x + sin ϑe_y and n = cos θe_x + sin θe_y are respectively the fluctuating and average orientation of the hexatic building blocks and the operator 〚…〛 renders its argument traceless and symmetric (see Refs. (Giomi et al., 2022b,a) 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

with D/Dt = ∂_t + υ ⋅∇ is the material derivative. In Eq. (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 to that Eq. (3a) reduces to the standard incompressibility condition ∇⋅υ = 0. In Eq. (3b) σ is the total stress tensor and f an external body force. In Eq. (3c), u = [∇〈〉+(∇〈〉)^T]/2 and ω = [∇υ−(∇υ)^T]/2, with T indicating transposition, are respectively the strain rate and vorticity tensors and entail the coupling between hexatic order and flow, with λ₆ and υ₆ material constants and (. Because of incompressibility, the term in Eq. (3c) vanishes. The tensor H₆ = −δF/δQ₆ 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

where |…|² is the Euclidean norm and is such that |Q₆|² = |Ψ₆|²/2. The constant L₆ is the order parameter stiffness, while A₆ and B₆ are phenomenological constants setting the magnitude of the coarse-grained complex order parameter at equilibrium: , when = 0.

The stress tensor figuring in Eq. (3b) can be customarily decomposed in a passive and an active contributions: i.e. σ = σ^(p) + σ^(a). The passive stress, in turn, can be expressed as σ^(p) = −P1 + σ^(v) + σ^(e) + σ^(r), where P is the pressure and σ^(v) = 2η〚u〛 the viscous stress, with η the shear viscosity. The tensor is the elastic stress, arising in response of 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) = −λ₆∇^⊗4⊙H₆ + 3 (Q₆ ⋅ H₆ − H₆ ⋅ Q₆) 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 Ref. (Armengol-Collado et al., 2023b) on the basis of phenomenological and microscopic arguments and is given by σ^(a) = α₆∇^⊗4⊙Q₆, with α₆ 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… N_cell, and N_cell 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 ℱ = ∫ dAf, where the free energy density f is

Here α > 0 and k_φ > 0 are material parameters which can be used to tune the surface tension and the interfacial thickness of isolated cells and thermodynamically favor spherical cell shapes, so that the concentration field is large (i.e. φ_c ≈ φ₀) 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 π_ad > 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 , with R_cell the preferential cell radius. The phase field φ_c evolves according to the Allen-Cahn equation

where υ_c = υ₀(cos θ_ce_x + sin θ_ce_y) is the velocity at which the c-th cell self-propels, with υ₀ a constant speed and θ_c the angle defining the nominal direction of cell migration. The latter evolves according to the stochastic equation

where η_c is a noise term with correlation function and D_r a constant controlling noise diffusivity. The constant M in Eq. (6) is the mobility measuring the relevance of thermodynamic relaxation with respect to non-equilibrium cell migration. Eq. (6) 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_cells = 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_ad = 0.005, λ = 600, R_cell = 11.5, υ₀ = 0.006. The noise variance D_r was varied in the range 10⁻⁶ ≤ D_r < 8 × 10⁻⁵.

Cell segmentation

In order to find the 6-fold 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

where ϑ_H is the Heaviside theta function such that ϑ_H(x) = 1 if x ≥ 0 and ϑ_H(x) = 0 otherwise. Here φ_th is the threshold marking the boundary of each cell. In particular we choose φ_th = φ₀/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 φ₀ to 0, therefore dropping below the threshold φ_th, it is possible to find pixels where Φ = 0. These spurious features are adjusted by replacing Eq. (8) with an average over the pixels neighboring that where Φ vanishes. That is

where Σ_(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_υ}_c’, as those points where φ_tot > 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 6-fold orientation of a cell is computed via the shape function /₆ defined in Eq. (2), introduced in Ref. (Armengol-Collado et al., 2023a). This construction is schematically illustrated in Fig. 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 Γ₆ = Γ₆(r), constructed upon averaging the shape function γ₆ of the segmented cells whose center of mass, r_c’, lies within a disk of radius R centered in r. That is

Figure 4

where N_disk = Σ_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.5R_cell. We choose to sample the shape function on a square grid with grid–spacing equal to the nominal cell radius R_cell. Topological defects are then identified computing the winding number along each unit cell:

where the symbol ◽ denotes a square unit cell in the interpolation grid and θ = Arg(Γ₆)/6 the phase of the shape order parameter.

Nematic and Hexatic stress in MPF simulations

To characterize the properties of intercellular stresses in MPF simulations, we show in Fig. 3a,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 with the cellular shape tensor of the corresponding order (, with p = 2, 6). Here, n_c,p = (cos θ_c,p, sin θ_c,p) with θ_c,p = Arg(γ_c,p)/p the angular orientation of the complex shape function γ_p relative to the c-th cell (see Fig. 4). Therefore, the explicit expression of the longitudinal nematic and hexatic stress is given by

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

To build the scatter plot shown in Fig. 3f correlating cellular mean square displacement and defect density, we proceed as follows. First, we divide the system in squared regions of size Δ = 35, containing roughly 3.5 cells. Then we observe each of this regions for a time window of 25 × 10³ iterations, tracking both the mean square displacement of the cells in each subregion and the mean defect density therein. Finally, we record these observation for each subregion and for each time window analyzed and we use these measurements as entries to build the scatter plot in Fig. 3f.

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 Fig. 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 Fig. 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 Fig. 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 Eq. (1). To this end, let Ψ₆ = |Ψ₆|e^6iθ be the hexatic complex order parameter and consider a quadrupole of ±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

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 Eq. (14) for |r|/ℓ ≫ 1, we obtain the simpler expression

Notice that the Taylor expansion features only the quadrupolar term of order ℓ² and is exact up to 6-th order in ℓ/|r|; i.e. the dipolar term, of order 𝒪(ℓ/|r|), and all other terms up to the 6-th order vanish identically.

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

where φ is a dual field such that ∂_iθ = -ε_ijε_jφ and the right-hand side of the Eq. (16) is analogous to the electrostatic charge density (Chaikin et al., 1995). At large distance from the defects, Eq. (16) can be solved by multipole expansion (Jackson, 1999), that is:

where r₀ is an irrelevant length scale and a_n and b_n are coefficients given by

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

As in electrostatics, the density ρ_d is given by

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 and a₁ = b₁ = 0. Of the quadrupolar terms, on the other hand, a₂ = -ℓ²/3 and b₂ = 0, thus

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

thus confirming the expression given in Eq. (15).

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) = ∇ · σ^(a), where

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

up to correction of order 𝒪(|ℓ/r|⁶). A plot of the force field is shown Fig. 5a.

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 f^(a): i.e.

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

where

is the two-dimensional Oseen tensor (see e.g. Ref. (Giomi, 2015)), with L a constant, and R is a large distance cut-off. Without loss of generality, one can set in Eq. (27). To calculate the integrals in Eq. (26), we make use of the logarithmic expansion

with r_≷ the maximum (minimum) between |r| and |r’|, and of the orthogonality of trigonometric functions

Figure 5

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

Fig. 5b shows a plot of this flow, while Fig. 5c shows a plot of the short distance approximation given in Eq. (1) of the main text.

Numerical simulations of defect annihilation and unbinding

Numerical model and validation

The time-dependent flows shown in Fig. 2a(i) and Fig. 2b(i) are obtained by numerically integrating Eqs. (3) 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 Refs. (Krommydas et al., 2023; Giomi et al., 2022b,a, 2014)). In all simulations, we set: ρ = 1, η = 1, L₆ = 0.5, A₆ = -0.2, B₆ = 0.4, Γ₆ = 1 and λ₆ = 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 Ψ₆, we set ℓ = 7, |Ψ₆| = 1 and take θ as given in Eq. (14) 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 Ψ₆ everywhere else. This allows us to obtain a defect-free configuration where |Ψ₆| ≈ 1 everywhere, except that close to the defect cores where |Ψ₆| ≈ 0. Notice that, on a doubly periodical domain, Σ_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 α₆ 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 Fig. 2a(i). Similarly, for positive values of α₆ we obtain increasingly sheared versions of the same flow pattern, but with the direction of the flow inverted.

Figure 6

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 ±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 Ψ₆ = e^6iθ by setting θ = ± arctan[y/(x±Δx/2)] inside a disk of radius R_D = 5 centred 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 Ψ₆ everywhere else. As before, this procedure allows us to obtain a state where |Ψ₆| ≈ 1 everywhere, except that close to the two defect cores where |Ψ₆| ≈ 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₆ = 0.5, A₆ = -0.2, B₆ = 0.4, Γ₆ = 1 and λ₆ = 1.11.

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

Acknowledgements

This work is supported by the European Union via the ERC-CoGgrant HexaTissue and by Netherlands Organization for Scientific Research (NWO/OCW). L. N. C. 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 produced the numerical data presented in the article are available upon request.