Hydrodynamics and multiscale order in confluent epithelia

  1. Josep-Maria Armengol-Collado
  2. Livio Nicola Carenza
  3. Luca Giomi  Is a corresponding author
  1. Instituut-Lorentz, Leiden University, Netherlands


We formulate a hydrodynamic theory of confluent epithelia: i.e. monolayers of epithelial cells adhering to each other without gaps. Taking advantage of recent progresses toward establishing a general hydrodynamic theory of p-atic liquid crystals, we demonstrate that collectively migrating epithelia feature both nematic (i.e. p = 2) and hexatic (i.e. p = 6) orders, with the former being dominant at large and the latter at small length scales. Such a remarkable multiscale liquid crystal order leaves a distinct signature in the system’s structure factor, which exhibits two different power-law scaling regimes, reflecting both the hexagonal geometry of small cells clusters and the uniaxial structure of the global cellular flow. We support these analytical predictions with two different cell-resolved models of epithelia – i.e. the self-propelled Voronoi model and the multiphase field model – and highlight how momentum dissipation and noise influence the range of fluctuations at small length scales, thereby affecting the degree of cooperativity between cells. Our construction provides a theoretical framework to conceptualize the recent observation of multiscale order in layers of Madin–Darby canine kidney cells and pave the way for further theoretical developments.

Editor's evaluation

This important work presents a hydrodynamic description of confluent epithelial monolayers that captures different forms of orientational order in a scale dependent fashion and couples this order with flows driven by active stresses. Solid evidence for the validity of this approach is provided by detailed numerical simulations of different model tissues. This work should be of interest to a broad range of biophysicists interested in tissue mechanics and active matter.



Collective cell migration – i.e. the ability of multicellular systems to cooperatively flow, even in the absence of a central control mechanism – has surged, in the past decade, as one of the central questions in cell biology and tissue biophysics (Friedl and Gilmour, 2009). Whether spreading on a synthetic substrate (Serra-Picamal et al., 2012) or invading the extracellular matrix (Haeger et al., 2020), multicellular systems can move coherently within their micro-environment and coordinate the dynamics of their actin cytoskeleton, while retaining cell–cell contacts. This ability lies at the heart of a myriad of processes that are instrumental for life, such as embryonic morphogenesis and wound healing, but also of life-threatening conditions, such as metastatic cancer.

Understanding the physical origin of this behavior inevitably demands reliable theoretical models, aimed at providing a conceptual framework for dissecting and deciphering the wealth of biophysical data stemming from in vitro experiments and in vivo observations. Following the pioneering works by Honda, 1978; Nagai and Honda, 2001; Farhadifar et al., 2007; Bi et al., 2015; Bi et al., 2016, and others (Boromand et al., 2018; Mueller et al., 2019; Loewe et al., 2020; Monfared et al., 2021), cell-resolved models have played so far the leading role in this endeavour. Taking inspiration from the physics of foams (Graner et al., 2008; Marmottant et al., 2008), these models portray a confluent tissue as a collection of adjacent or overlapping polygonal cells (Figure 1a, b), whose dynamics is assumed to be governed by a set of overdamped Langevin equations, expressing the interplay between cells’ autonomous motion and remodeling events, which change the local topology of the cellular networks.

Epithelia and hexanematic order.

(a) Example of multiscale hexanematic order in an in vitro layer of Madin–Darby canine kidney (MDCK) cells (b) and its computer-constructed segmentation. Both panels are adapted from Figure 3 of Armengol-Collado et al., 2023. The six-legged stars in the shaded region denote the sixfold orientation of the cells obtained using the approach summarized in Methods. The colored stripes mark the configuration of the nematic director at the length scale of the light-blue disk. (c) Schematic representation of the sixfold symmetric force complexion exerted by cells. The red arrows indicate the structure of the contractile forces acting within the cellular junctions.

Despite their conceptual simplicity, cell-resolved models agree remarkably well with experimental data on confluent monolayers (Park et al., 2015; Atia et al., 2018). In particular they account for a solid-to-liquid transition controlled by the cells velocity and their compliance to deformations (Bi et al., 2015; Bi et al., 2016; Loewe et al., 2020). Furthermore, as demonstrated by Pica Ciamarra and coworkers, the solid and isotropic liquid states of these model-epithelia are separated by an intermediate hexatic phase, in which the system exhibits the typical sixfold rotational symmetry of two-dimensional crystals and yet is able to flow (Li and Ciamarra, 2018; Pasupalak et al., 2020). Shortly after discovery, the same property has been recovered within the framework of the cellular Potts model, thereby strengthening the idea that hexatic order may in fact serve as a guiding principle to unravel the collective dynamics of confluent epithelia (Durand and Heu, 2019). Furthermore, recent in vitro studies of Madin–Darby canine kidney (MDCK) cell layers demonstrated that epithelial layers can in fact feature both nematic and hexatic orders, with the former being dominant at large and the latter at short length scales (see Figure 1a, b and Armengol-Collado et al., 2023; Eckert et al., 2023). This remarkable example of physical organization in biological matter, referred to as multiscale hexanematic order in Armengol-Collado et al., 2023, is believed to complement the complex network or regulatory pathways available to individual cells to achieve multicellular organization and select specific scale-dependent collective migration strategies.

Motivated by these recent discoveries, in this article we propose a continuum theory of confluent epithelia rooted in the hydrodynamics of liquid crystals with generic p-atic rotational symmetry (hereafter p-atic liquid crystals). Previous theories of epithelial hydrodynamics can be schematically grouped in two categories: (1) models based on (isotropic/polar/nematic) active gels (Ranft et al., 2010; Popović et al., 2017; Pérez-González et al., 2019); (2) models built around the so-called shape tensor (Ishihara et al., 2017; Czajkowski et al., 2018; Hernandez and Marchetti, 2021; Grossman and Joanny, 2022), i.e. a rank-2 tensor, similar to the inertia tensor in kinematics, that embodies the geometrical structure of the polygonal cells. Although both classes of models hold great heuristic value and represent a solid foundation for any future development, they suffer from the same limitation: being based on a tensorial order parameter whose rank is two or less, they can account at most for twofold rotational symmetry (i.e. nematic order), while leaving the small-scale hexatic order unresolved. To overcome this limitation, here we exploit recent advances toward extending the classic hydrodynamic theory of hexatic liquid crystals (Zippelius, 1980; Zippelius et al., 1980) to account for arbitrary p-fold rotational symmetry order (Giomi et al., 2022a; Giomi et al., 2022b), with p=2 and p=6 being the most relevant cases (but possibly not the only) in the context of epithelial dynamics. We demonstrate that multiscale order is inherent to active liquid crystals with coupled order parameters, because of the indissoluble connection between shape and forces characterizing this class of non-equilibrium systems. Using fluctuating hydrodynamics, we explicitly compute the structure factor of epithelial layers and unveil a fascinating interplay between the nature of momentum dissipation (i.e. viscosity or friction) and noise at short length scales, where hexatic order is dominant. Such a mechanism profoundly affects the range of density fluctuations and could be harnessed to control the degree of collectiveness of cellular motion. Finally, by testing predictions against two different microscopic models of epithelia we demonstrate the robustness of multiscale hexanematic order across the rich landscape of models of epithelia.

Results and discussion

The model

Two-dimensional p-atic liquid crystals are traditionally described in terms of the orientation field ψp=eipϑ, with ϑ the local orientation of the p-fold mesogens. A more general approach, proposed in Giomi et al., 2022a; Giomi et al., 2022b and especially suited for hydrodynamics, revolves instead around the rank-p tensor order parameter, Qp=Qi1i2ipei1ei2eip with in={x,y} and n=1,2p, constructed upon averaging the pth tensorial power of the local orientation ν=cosϑex+sinϑey. That is

(1) Qp=2p2[[νp]]=2p2|Ψp|[[np]],

where denotes the ensemble average and the operator ... has the effect of rendering an arbitrary tensor traceless and symmetric (Hess, 2015). The vector n=cosθex+sinθey is the analog of the director field in standard lexicon of nematic liquid crystals and marks the average cellular direction, which in turn is invariant under rotations of 2π/p. The fields |Ψp| and θ represent, respectively, the magnitude and phase of the complex p-atic order parameter Ψp=ψp, while the normalization factor is chosen so that |Qp|2=|Ψp|2/2 for all p values. For p=2, Equation (1) readily gives the standard nematic order parameter tensor: i.e. Q2=|Ψ2|(nn1/2), with 1 the identity tensor. In practice, if a cell’s planar projection consists of a regular p-sided polygon, the microscopic orientation ϑ equates that of any of the vertices of the polygon. In the more realistic case of an irregular polygon, on the other hand, ϑ is given by the phase of the complex function γp, arising form the p-fold generalization of the classic shape tensor (Aubouy et al., 2003). This function was introduced in Armengol-Collado et al., 2023 and is reviewed in Methods for sake of completeness.

The order parameter tensor Qp, the mass density ρ, and the momentum density ρv, with v the local velocity field, comprise the set of hydrodynamic variables describing the dynamics of a generic p-atic fluid, which in turn is governed by the following set of partial differential equations (Giomi et al., 2022a; Giomi et al., 2022b):

(2a) DρDt+ρv=(kdka)ρ,
(2b) ρDvDt=σ+f,
(2c) DQpDt=ΓpHp+p[[Qpω]]+λ¯ptr(u)Qp+λp[[(p2)u]]+νp[[(pmod2)up/2]]

where D/Dt=t+v. Equation (2a) and Equation 2b are the mass and momentum conservation equations, with kd and ka rates of cell division and apoptosis, σ the stress tensor and f an arbitrary external force per unit area. In Equation (2c), Γp1 is a rotational viscosity and Hp=δF/δQp is the molecular tensor describing the relaxation of the p-atic phase toward the minimum of the free energy F (see Methods). The rank-2 tensors ω=[v(v)]/2 and u=[v+(v)]/2, with indicating transposition, are the vorticity and strain rate tensors, respectively, whereas the dot product in the first line of the equation implies a contraction of one index of Qp with one of ω: i.e. (Qpω)i1i2ip=Qiii2jωjip. On the second line (n)i1i2in=i1i2in, while denotes the floor function and pmod2=p2p/2 is zero for even p values and one for odd p values. Finally, λ¯p, λp, and νp are material parameters expressing the strength of the coupling between p-atic order and flow.

Now, in order for Equation 2a to account for the dynamics of epithelial cell layers, we must specify the structure of the external force f in Equation 2b and the stress tensor σ. As cells collectively crawl on a substrate, at a speed of order 0.1–1 µm/min (Brugués et al., 2014; Angelini et al., 2011), the former can be model as a Stokesian drag: f=ςv, with ς a drag coefficient. A more realistic treatment of the interplay between the cells and the substrate would account for the traction forces exerted by the cells’ cryptic lamellipodium as well as for the compliance of the substrate (Trepat et al., 2009) and will be considered in the future. The stress tensor, on the other hand, is routinely decomposed into a passive and an active component: i.e. σ=σ(p)+σ(a). The passive stress tensor is in turn expressed as σ(p)=P1+σ(e)+σ(r)+σ(v), where P is the pressure, σ(e) is the elastic stress, arising in response of a static deformation of a fluid patch, and σ(r) and σ(v) are, respectively, the reactive (i.e. energy preserving) and viscous (i.e. energy dissipating) stresses originating from the reversible and irreversible couplings between p-atic order and flow. The generic expression of σ(p) was derived in Giomi et al., 2022b and is reported in Methods.

The active stress σ(a), on the other hand, can be constructed phenomenologically for arbitrary p values in the form

(3) σ(a)=p(αp(p2)Qp+βp[[2|Qp|2]]),

where the symbol denotes a contraction of all matching indices of the two operands and yields a tensor whose rank equates the number of unmatched indices: i.e. letting Ap and Bq be two generic tensors of rank p<q, then (ApBq)i1i2iqp=Aj1j2jpBj1j2jpi1i2iqp. The sum over p, finally, reflects the possibility of having not only one, but multiple types of p-atic order coexisting within the same system, as experiments on in vitro layers of MDCK cells have recently suggested (Armengol-Collado et al., 2023; Eckert et al., 2023).

Before exploring the consequences of the latter assumption, some comment about the physical interpretation of the terms featured in Equation 3 is in order. The first term on the right-hand side of Equation 3 is the stress resulting from the contractile or extensile forces exerted at the length scale of individual cells. To illustrate this concept one can assume each cell to exert a p-fold symmetric force complexion: i.e. Fc=k=1pFkδ(rrcaνk) with Fk the force exerted by a cell at each vertex and originating from the imbalance of the tensions Tkl, driven by the active contraction of the cellular junctions, converging at the kth vertex: i.e. Fk=lTkl (see Figure 1c). The quantities rc and a are the cell’s centroid and circumradius, respectively, while νk=cos(ϑ+2πk/p)ex+sin(ϑ+2πk/p)ey. We stress that, while the individual tensions acting along the junctions are exclusively contractile, the resulting vertex forces can be either contractile (i.e. Fkνk<0) or extensile (Fkνk>0), depending on the overall tension distribution and the geometry of the cellular network. Next, assuming Fk=fνk and expanding the delta function about a=0 yields Fc=m=0fm, where

(4) fm=m[(a)mfm!(k=1pνk(m+1))δ(rrc)].

Because of the p-fold symmetry of the force complexion fm=0 for all even m values, unless m=p1, whereas odd m values yields, up to symmetrization, k=1pνk(m+1)1(m+1)/2. Thus, after some algebraic manipulation, one finds Fcapf/2[(1+a2/82+)δ(rrc)]+fp1. Finally, taking cFc=P(a)1+σ(a) gives the following expression for contributions to the pressure and the deviatoric stress resulting from the active expansion and contraction of the cells. That is

(5a) P(a)=apf2(n+a282n+),
(5b) σ(a)=(-a)p-1pnf2p-2(p-1)!(p-2)Qp,

where n=cδ(rrc) is the cell number density. From Equation 5b, one finds the following expression for the phenomenological parameter αp in Equation 3: i.e. αp=(a)p1pnf/[2p2(p1)!]. Notice that both constants a and f involved in Equation 5a are, in general, order dependent. We will come back on this aspect in Conclusion.

The second term in Equation 3, in contrast, expresses the active stress resulting from the spatial variations of the p-atic order parameter and, although similar to other contributions to the passive stress σ(p), cannot be derived from equilibrium considerations. Other terms constructed by contracting Qp with 2 can be expressed as linear combinations of this and σ(p), thus lead to a mere renormalization of the material parameters. It must be noted that the stress tensor enters in Equation 2b only via its divergence. Thus, possible second-order active terms such as Qk1k2kpijQk1k2kp, Qijk3kpl1l2Ql1l2k3kp, etc., are mechanically equivalent to the terms iQk1k2kpjQk1k2kp and Qk1k2iHk1k2jHk1k2iQk1k2j arising from the passive stresses, as both sets of terms lead to the same body forces.

We observe that Equation 3 already entails a multiscale hydrodynamic behavior even when a single p value is considered. Such a crossover is expected at length scales larger than =(αp/βp)1/(p4), where the second term of the right-hand side of Equation 3 overweights the first term, reflecting the p-fold symmetry of the local active forces. In the presence of multiple types of p-atic order, the p-dependent structure of the active stress renders the multiscale nature of the system enormously more dramatic. To illustrate this crucial point, here we postulate the system to behave as a hexanematic liquid crystal. Formally, such a scenario can be accounted by simultaneously solving two variants of Equation 2c, for Q2 and Q6. In turn, the interplay between nematic and hexatic order results from a combination of dynamical and energetic effects. The former arise from active flow, which affects the local configuration of both tensor order parameters via the last four terms in Equation 2c. The latter, instead, can be embedded into the free energy F=dA(f2+f6+f2,6), where

(6a) fp=12Lp|Qp|2+12Ap|Qp|2+14Bp|Qp|4,
(6b) f2,6=κ2,6|Q2|2|Q6|2+χ2,6Q23Q6.

Here, Ap and Bp are constants setting the magnitude of the order parameter at the length scale of the short distance cut-off, here assumed to be of the order of the cell size, and κ2,6 determines the extent to which the magnitude of the hexatic order parameter is influenced by that of the nematic order parameter and vice versa. The constant χ2,6, on the other hand, is analogous to an inherent susceptibility, expressing the propensity of the nematic and hexatic directors toward mutual alignment. The free energy contribution f2,6 can further be augmented with several additional terms of higher differential order: e.g. (Q2Q2)(Q6Q6), |(Q23Q6)|2, 2(Q23Q6), etc. For simplicity, here we ignore these and higher-order couplings and focus on the zeroth order terms included in Equation 6b.

Crucially, Equation 3, Equation 6a entail two length scales, reflecting the distance at which the passive torques originating from the entropic elasticity of the nematic and hexatic phases counterbalance those arising from the active stresses:

(7) 2=L2|α2|,6=|α6|L6.

The former is the well-known active nematic length scale, dictating both the hydrodynamic stability (Voituriez et al., 2005) and the large-scale structure of spatiotemporal chaos in active nematics (Giomi, 2015) and whose signature in multicellular systems has been identified in both eukaryotes (Blanch-Mercader et al., 2018) and prokaryotes (You et al., 2018). The latter, on the other hand, sets the typical size of hexatic domains at the small length scale. Remarkably, 2 and 6 inversely depend on the magnitude of cellular forces (see Equation 5a). Thus, increasing activity has the effect of collapsing the multiscale structure of the system toward a single length scale, where 26. Two additional length scales, of purely passive nature, originate from the competition between rotational diffusion and the ordering dynamics driven by either liquid crystalline structure on the other one. These are given by χ,2=L2/χ2,6 and χ,6=L6/χ2,6. Their role will be discussed in the following section, in the framework of fluctuating hydrodynamics.

Finally, in the passive limit, when α2=0 and α6=0, Equations 2 and 6, reduce to those of a two-dimensional liquid crystal with coupled nematic and hexatic order parameter. The latter can be found, e.g., in free-standing liquid hexatic films (Dierker and Pindak, 1987; Sprunt and Litster, 1987), where molecules are either orthogonal to the mid-surface of the film or tilted by a fixed angle. In the latter case, the projection of the average molecular direction on the tangent plane of the mid-surface gives rise to in-plane nematic order, which is coupled to the sixfold bond-orientational order associated with the underlying hexatic phase (see e.g. Bruinsma and Aeppli, 1982; Selinger and Nelson, 1989; Selinger, 1991 for a theoretical account and Drouin-Touchette et al., 2022 for recent developments). As we will detail in the following, activity profoundly alters this scenario by acting as a mechanical bandpass filter, which renders hexatic order dominant at length scales 6 and nematic order at length scales 2. We stress that by dominant, here we intend able to drive morphological features, dynamical behaviors, and fluctuations reflecting the underlying orientational order. At intermediate length scales, i.e. 62, there is no dominant order and the system’s collective behavior is determined by the complex interplay of competing active and passive effects. To make progress, here we focus on the most dramatic hexatic- and nematic-dominated behaviors and treat intermediate length scales as simply as possible.

Multiscale order in epithelia

To elucidate the multiscale organization of the system, we next compute the structure factor S(|q|), using the classic framework of fluctuating hydrodynamics (see e.g. Ramaswamy et al., 2003). To this end, we assume both the nematic and the hexatic scalar order parameters to be uniform throughout the system and set kd=ka and λp=0 for simplicity. We stress that the validity of this approximation is strictly related with the present comparison between the hydrodynamic theory presented in this article and cell-resolved models. An assessment of the relevance of this and the other material parameters featured in Equation 2a can only be achieved via experimental scrutiny and is likely to depend on the specific cell type and environmental conditions. Furthermore, as the typical Reynolds number of collective epithelial flow is in the range 10−7–10−6, we neglect inertial effects: i.e. ρDv/Dt=0. With these simplifications, whose legitimacy will be assessed a posteriori, one can reduce Equation 2b to three coupled differential equations for the density and the phases of the hexatic and nematic order parameter tensors (see Methods). These equations, in turn, can be linearized about the trivial configuration, where all fields are spatially uniform and v=0, and augmented with noise terms to give the following exact asymptotic expansion

(8) S(|q|)s2|q|2+sβ|q|β.

The first term entails the typical giant number density fluctuations associated with the active nematic behavior at the large scale, with s2α22. This effect is overestimated at the linear order, leading to an inverse quadratic dependence on the wave number |q| (Ramaswamy et al., 2003), but is generally renormalized by nonlinearities, so that lim|q|0S(|q|)|q|α, with 1<α<2 (Shankar et al., 2018; Chaté, 2020).

The second term, on the other hand, reflects the sixfold symmetry characterizing the structure of epithelia at the small length scale, with sβα62 and the exponent β determined by the specific energy dissipation mechanism, as well as by the specific structure of the noise. As detailed in Methods, here we consider four alternative scenarios, obtained upon combining two different momentum dissipation mechanisms (i.e. viscosity and friction) with two different types of noise (i.e. rototranslational and purely rotational). In the presence of viscous dissipation, i.e. a regime referred to as ‘wet‘ in the jargon of active matter, β=4 irrespective of the nature of noise. Conversely, in the ‘dry‘ limit, when the shear and bulk viscosity vanish and momentum dissipation solely results from the frictional interactions with the substrate, β differs depending on whether noise affects both cells’ orientational and translational dynamics, or only the former. Specifically, when only orientational noise is considered, β=6. In contrast, β=10 in the presence of conservative rototranslational noise. We again stress that Equation (8) is an exact asymptotic expansion, as one could verify upon comparison with the full analytical solutions plotted in Figure 2, and not a truncated power series.

Structure factor S(|q|) obtained from the analytical solutions of the linearized hydrodynamic equations in the presence of two different noise fields: purely rotational (blue) and rototranslational (red).

The full analytical expression of S(|q|) is given in Methods, together with a derivation of the exact asymptotic expansions of Equation (8). (a) As long as viscous dissipation takes place (i.e. ‘wet’ regime), S(|q|)|q|4 in the limit |q|, irrespective of the type of noise. (b) On the other hand, when friction is the sole momentum dissipation mechanism at play (‘dry’ regime), S(|q|)|q|6 in case of rotational noise and S(|q|)|q|10 when noise is both rotational and translational. In both panels, the wave number |q| is rescaled by q¯=2π/¯, with ¯=(2+6)/2 and 2 and 6 as defined in Equation (7).

To test the significance of these predictions and connect the present hydrodynamic theory with the existing literature, in Figure 3a we compare the structure factor obtained from numerical simulations of two different cell-resolved models of epithelia – i.e. the self-propelled Voronoi (SPV) model (Bi et al., 2016) and the multiphase field (MPF) model (Loewe et al., 2020) (see the insets Figure 3b for typical configurations of the two models) – with that resulting from a numerical integration of Equation 2a (Carenza et al., 2019; Carenza et al., 2020), with none of the simplifications behind Equation (8). In both microscopic models, cells are treated as persistent random walkers, self-propelling at constant speed v0 and whose direction of motion undergoes rotational diffusion with diffusion coefficient Dr (see Methods for details). Noise is therefore expected to affect both the rotational and translational dynamics of the cell monolayer, although in a way that, unlike in our analytical treatment, cannot be trivially decoupled. Consistently with our linear analysis, both data sets exhibit two different power-law scaling regimes at small and large length scales. At small length scales, the structure factor scales like S(|q|)|q|β, with β monotonically decreasing from 6 to 4 upon increasing the Péclet number Pe=ξ0/a expressing the ratio between cells’ persistence length ξ0=v0/Dr and their typical size a (see Figure 3b).

Density fluctuations in model epithelia.

(a) Structure factor of model-epithelia calculated from a numerical integration of Equation 2a (black line) and from simulations of two different cell-resolved models: i.e. the self-propelled Voronoi model (SPV, red) and the multiphase field (MPF) model (blue), for a particular choice of parameters. The dashed diagonal lines mark the scaling regimes obtained analytically at the linear order, Equation (8), and the wave number |q| is rescaled by qcell=2π/ΔxLB, where ΔxLB is the grid size used by the Lattice Boltzmann integrator (see Methods for details). (b) The exponent β, as defined in Equation (8), versus the Péclet number Pe, reflecting the persistence of directed cellular motion in front of diffusion. Error bars calculated as standard error over $n=250$ configurations for both the SPV and the MPF models. Insets: typical configurations of the SPV (bottom left) and MPF (top right) models.

Conversely, at large length scales, the structure factor scales like an inverse power law, with exponent consistent with the large-scale behavior of active nematics (Chaté, 2020). These observations can be rationalized in the light of the previous fluctuating hydrodynamic analysis. In the limit Pe0, cells do not self-propel, noise is predominantly orientational and momentum propagates only at distances comparable to the average cell size. Under this circumstances, an in silico cell layer, whether modeled via the SPV or the MPF, behaves therefore as a ‘dry’ active system subject to purely rotational noise, for which, consistently with our analysis, β=6. Increasing Pe has the twofold effect of converting noise from purely rotational to rototranslational and, by stimulating cooperativity in the cellular motion, to increase the range of momentum propagation, thus driving a crossover of the cell layer from ‘dry’ to ‘wet’, hence from β=6 to β=4. The simple linear calculation, summarized in Methods, does not allow us to resolve the full crossover, but does provide a precise estimate of the upper and lower bounds. Finally, along the wet–dry crossover, viscosity must emerge from the cells’ lateral interactions. A precise understanding of this process is outside of the scope of the present work, but recent numerical work on the Vertex model has already highlighted the existence of a rich landscape of exotic rheological phenomena, resulting from the interplay between cellular motion, morphology, and adhesion (Tong et al., 2022; Hertaeg et al., 2022). The latter could possibly explain the non-monotonic behavior at small Pe values, as a crossover from a shear-thinning to the shear-thickening behavior (Hertaeg et al., 2022) for additional numerical evidence of this effect.

A different signature of multiscale hexanematic order can be identified in the structure of the cross-correlation function

(9) C26(r)=ψ2(r)ψ6(0)+ψ2(r)ψ6(0)2.

At equilibrium, and if deformations are sufficiently gentle to render backflow effects negligible, its behavior can be divided in two regimes, depending on how the distance |r| compares to the length scales χ,2 and χ,6 defined in the previous section and expressing the typical distance at which the mutual alignment rate of the hexatic and nematic orientations overcome that of rotational diffusion. In the simplest possible setting, when χ,2=χ,6=χ, fluctuations dominate at short distances and the hexatic and nematic orientations are uncorrelated. Thus, C26(r) is approximatively constant for |r|χ. The picture is reversed for |r|χ. In this range, the hexatic and nematic orientations are ‘locked’ in a parallel configuration, i.e. Arg(ψ2)/2Arg(ψ6)/6, or tilted by π/6 with respect to each other, depending on the sign of the constant χ2,6, and the cross-correlation function exhibits the standard power-law decay characterizing two-dimensional liquid crystals with a single-order parameter: i.e. C26(r)(|r|/χ)η26, with η26 a specific instance of the generic non-universal exponent η26=6kBT/(πK), with K the orientational stiffness of both phases (proportional to L2=L6). An analytical treatment of this simple case is reported in Methods. In the more generic case, in which χ,2χ,6 and the relaxation rates of the hexatic and nematic phase differ, the cross-correlation function has a less standard functional form, but still features a slow and fast decay regime at short and large distances, respectively. An example of such a scenario, obtained from a numerical integration of Equation 2a with α2=0 and α6=0, is shown in Figure 4a. The curves in Figure 4b correspond instead to simulated configurations of the cross-correlation function of C26(r) for finite hexatic and nematic activity. In this case, the cross-correlation function exhibits an oscillatory behavior at short distances and vanishes at a length scale that becomes progressively large as the hexatic activity is increased. Consistently with our previous analysis, this latter feature confirms the existence of a hierarchy of orientationally ordered structures nested into each other at different length scales.

Cross-correlation function C26(r), as defined in Equation 9, obtained from a numerical integration of Equation 2a augmented with rotational noise.

(a) In the passive case, when α2=0 and α6=0, the correlation function decays with |r| at a rate i.e. lower at short distances, where the dynamics of the hexatic and nematic orientations is dominated by fluctuations, and larger at long distances, where the orientations are ‘locked’ in a parallel configuration, or tilted by π/6 with respect to each other. (b) In the active case, conversely, the cross-correlation function has a damped oscillatory behavior. Consistently with Equation (7) and the related discussion, the range of the oscillations, corresponding to the distance at which these are fully damped, increases with the hexatic activity α6, indicating an enhancement of hexatic order at larger length scales. Shaded region corresponding to the standard deviation of n=500 configurations. Distance is expressed in terms of the grid size ΔxLB used by the Lattice Boltzmann integrator (see Methods for details).

Taken together, our calculations of the structure factor and the cross-correlation function demonstrate that the hydrodynamic theory embodied in Equations 2 and 6 is able to account for the multiscale hexanematic order observed in experiments (Armengol-Collado et al., 2023; Eckert et al., 2023) and harnesses it into a continuum mechanical framework. Whereas the origin of hexanematic order is still a matter of investigation, the current experimental and numerical evidence suggests that, similarly to granular materials (Majmudar and Behringer, 2005), large-scale nematic order could arise from the self-organization of the microscopic force hexapoles into force chains. The possibility of similarity between these two phenomena has also been in relation to the initial phase of Drosophila gastrulation, where linear arrays of cells simultaneously undergo apical constriction in the ventral furrow region (Jason Gao et al., 2016).


In conclusion, we have introduced a continuum model of collectively migrating layers of epithelial cells, built upon a recent generalization of the hydrodynamic theory of p-atic liquid crystals (Giomi et al., 2022a; Giomi et al., 2022b). This approach allows one to account for arbitrary discrete rotational symmetries, thereby going beyond existing hydrodynamic theories of epithelia (Ranft et al., 2010; Popović et al., 2017; Pérez-González et al., 2019; Ishihara et al., 2017; Czajkowski et al., 2018; Hernandez and Marchetti, 2021; Grossman and Joanny, 2022), where the algebraic structure of the hydrodynamic variables renders impossible to account for liquid crystal order other than isotropic (i.e. p=0), polar (i.e. p=1), or nematic (i.e. p=2). Upon computing the static structure factor and comparing this with the outcome of two different cell-resolved models – i.e. the SPV (Bi et al., 2016) and MPF (Loewe et al., 2020) models – we have shown that, consistently with recent experimental findings (Armengol-Collado et al., 2023; Eckert et al., 2023), epithelial layers may in fact comprise both nematic and hexatic (i.e. p=6) order, coexisting at different length scales. Although the consequences of such a remarkable versatility are yet to be explored, we expect hexatic order to be relevant for short-scale remodeling events, where the local nature of hexatic order, combined with the rich dynamics of hexatic defects (Zippelius et al., 1980; Amir and Nelson, 2012), may mediate processes such as cell intercalations and the rearrangement of multicellular rosettes (Blankenship et al., 2006; Rauzi, 2020). Such a local motion, in turn, may be coordianted at the large scale by the underlying nematic order, giving rise to a persistent unidirectional flow, such as that observed during wound healing and cancer progression (Friedl and Gilmour, 2009). Furthermore, the existence of multiscale liquid crystal order echoes the most recent understanding of phenotypic plasticity in tissues, according to which the epithelial (i.e. solid-like) and mesenchymal (i.e. liquid-like) states represent the two ends of a spectrum of intermediate phenotypes (Zhang and Weinberg, 2018). These intermediate states display distinctive cellular characteristics, including adhesion, motility, stemness and, in the case of cancer cells, invasiveness, drug resistance, etc. Can multiscale liquid crystal order help understanding how the biophysical properties of tissues vary along the epithelial–mesenchymal spectrum? This and related questions will be addressed in the near future.


Quantification of p-atic order in epithelial layers

Following Armengol-Collado et al., 2023, we use the shape functionγp to quantify the amount of p-fold symmetry of an arbitrary cell. Denoting rv with v=1,2V, the positions of its vertices with respect to the cell’s center of mass, one has

(10) γp=v=1V|rv|peipϕvv=1V|rv|p,

with ϕv=Arg(rv) the angle between rv and the x-axis of a Cartesian frame. A schematic representation of these elements in an arbitrary irregular polygon is shown in Figure 5a. Unlike the complex function ψp=eipϑ, which has unit magnitude by construction, the magnitude |γp| quantify the resemblance of a generic polygon with a regular p-sided polygon of the same size, while the phase ϑ=Arg(γp)/p marks the orientation of the polygon. For regular V-sided polygons, |γp|=1 provided p is an integer multiple of V and |γp|0 otherwise. Furthermore, from γp one can readily compute.

(11) ψp=γp|γp|.
Nematic and hexatic shape function.

(a) Irregular polygonal cell with a red cross marking its center of mass and rv and ϕv the radial vector and the angle to one of the six vertices, respectively. (b) and (c) show the same tessellation of the plane with cells of different shapes and the shape analysis using the function in Equation 10 for the nematic (p=2) and hexatic (p=6) case. Rods and stars are oriented according to the phase of γp and the color corresponds to its magnitude.

Figure 5b, c shows examples of the functions γ2 and γ6 for a typical configuration of the SPV. We emphasize that γp, which, as shown in Armengol-Collado et al., 2023, arises from a p-fold generalization of the classic shape tensor (Aubouy et al., 2003), is solely determined by the positions of the vertices of an individual polygon and, therefore, does not depend on the spatial organization of the neighboring cells. As a consequence, this approach establishes an orientation purely based on cellular shape, thereby eliminating the arbitrariness involved with associating a network of bonds to a planar tessellation, where the latter is not inherent.

The shape function γp can then be coarse grained at the length scale to construct the shape parameter:

(12) Γp(r)=1Nc=1γp(rc)Θ(|rrc|),

where the rc is the position of the cth cell, Θ is the Heaviside step function, such that Θ(x)=1 for x>0 and 0 otherwise, and N=cΘ(|rrc|) is the number of cells within a distance from rc. As in the case of γp, the magnitude of Γp reflects the resemblance between a multicelluar cluster and a regular p-sided polygon, while its phase marks the cluster’s global orientation. The outcome of an application of this method to the Voronoi model is illustrated in Figure 6 for p=2. The different patches in panel (a) are regions with uniform θ=Arg(Γ2)/2, while in panel (b), there are plotted streamlines showing the orientation of the director n=cosθex+sinθey.

Hexanematic symmetry of Voronoi model.

(a) Coarse-grained nematic orientation θ obtained from averaging the local shape of cells over domains of size 30cell, with cell the average size of individual cells. Regions with the same color represent domains of coherent nematic orientation. (b) Part of the system where we use Γ2 to characterize the nematic phase. Solid lines represent the nematic director and the color inditicates the magnitude of the nematic shape function. (c) Voronoi cell structure of a region where the nematic field is uniform. Polygons are colored according to |γ6| and the stars are oriented according to Arg(γ6)/6.

Passive stresses

As explained in the main text, the passive contribution to the stress tensor is given by σ(p)=P1+σ(e)+σ(r)+σ(v), where, as demonstrated in Giomi et al., 2022b

(13a) σij(e)=-LpiQpjQp,
(13b) σij(r)=λ¯pQpHpδij+(1)p1λpk1k2kp2p2Hk1k2ij+p2(Qk1k2iHk1k2jHk1k2iQk1k2j),
(13c) σij(v)=2η[[uij]]+ζtr(u)δij,

where η and ζ are, respectively, the shear and bulk viscosity and the other material parameters are defined in the main text. Under the assumptions of uniform order parameter, i.e. |Qp|2=|Ψp|2/2=const, and taking λp=0, Equation 13a reduces to the expression derived in Zippelius, 1980; Zippelius et al., 1980. That is

(14) σ(e)+σ(r)=P1+Kp2ε2θKpθθ,

where the first term in Equation 13b has incorporated into the pressure P and Kp denotes the orientational stiffness of the p-atic phase, related to the order parameter stiffness by

(15) Kp=p2|Ψp|22Lp

and varepsilon is the two-dimensional antisymmetric tensor, with εxy=εyx=1 and εxx=εyy=0.

Linear fluctuating hydrodynamics

To compute the structure factor, we follow Ramaswamy et al., 2003 and augment Equation 2b, Equation 2c with short-ranged correlated noise field. Then calling ϑ and φ the nematic and hexatic fluctuating orientation fields and linearizing the hydrodynamic equations about the homogeneous and stationary solutions, ϑ=φ=0 and v=0, gives

(16a) tδρ=ρ0δv,
(16b) tδϑ=D22δϑ+12ez(×δv)+94χ2(δϑδφ)+ξ(ϑ),
(16c) tδφ=D62δφ+12ez(×δv)+14χ6(δφδϑ)+ξ(φ),

where δϑ, δφ, and δv indicate a small departure from the homogeneous and stationary configurations of the fields ϑ, φ, and v, Dp=ΓpLp, χp=Γpχ2,6, and ξ(ϑ) and ξ(φ) are short-ranged correlated noise fields: i.e.

(17) ξ(α)(r,t)ξ(β)(r,t)=2(Ξ(ϑ)δαϑδβϑ+Ξ(φ)δαφδβφ)δ(rr)δ(tt).

The velocity field δv, on the other hand, is found from the Stokes limit of Equation 2b in the main text, which, at the linear order in all fluctuating fields, takes the form

(18) η2δv+ζ(δv)ςδv+f(p)+f(a)+ξ(v)=0.

where f(p)=σ(p) and f(a)=σ(a) are the body forces resulting from the passive and active stresses, respectively. The quantity ξ(v) is a translational noise field. In the absence of external stimuli, it is reasonable to assume that global momentum is neither created nor dissipated by translational fluctuations, but only redistributed across the cell layer. Thus ξ(v) is either conservative or null, from which

(19) ξi(v)(r,t)ξj(v)(r,t)=2Ξ(v)δij(2)δ(rr)δ(tt),

with {i,j}{x,y} and the case of noiseless translational dynamics, corresponding to Figure 3 in the main text, is recovered in the limit Ξ(v)0. The pressure P, in turn, can be related to the density by a linear equation of state of the form

(20) P=cs2ρ,

with cs the speed of sound. Together with the expression for the active stress given in Equation 3 of the main text, this gives

(21a) f(p)=(cs2xδρ+K22y2δϑ+K62y2δφ)ex(cs2yδρ+K22x2δϑ+K62x2δφ)ey,
(21b) f(a)=[α2yδϑ+32α6(y45x2y2+52x4)yδφ]ex+[α2xδϑ+32α6(x45x2y2+52y4)xδφ]ey.

Now, in Fourier space Equation 18 can be cast in the form of the following linear algebraic equation

(22) [(η|q|2+ς)1+ζqq]δv^=f^(p)+f^(a)+ξ^(v),

where the hat denotes Fourier transformation. Next, using

(23) [(η|q|2+ς)1+ζqq]1=[(η+ζ)|q|2+ς]1ζqq(η|q|2+ς)[(η+ζ)|q|2+ς],

and solving Equation 22 and incorporating the resulting velocity field in Equation 16a gives, after several algebraic manipulation

(24) iω[δρ^δϑ^δφ^]=M^[δρ^δϑ^δφ^]+[η^(ρ)η^(ϑ)η^(φ)],

where the matrix M^ is given by


and the functions η(α), with α{ρ,ϑ,φ}, are effective noise fields whose correlation functions are given by

(25) η^(α)(q,ω)η^(β)(q,ω)=(2π)32H^(α)(q)δαβδ(q+q)δ(ω+ω),

where the functions H^(α)=H^(α)(q) are given by

(26) H^(ρ)=ρ02|q|4[(η+ζ)|q|2+ς]2Ξ(v),
(27) H^(α)=Ξ(ϑ)δαϑ+Ξ(φ)δαφ+|q|44(η|q|2+ς)2Ξ(v).

Notice that, while hydrodynamic flow has the effect of coloring the orientational noise embodied in the stochastic fields ξ(ϑ) and ξ(φ), via the vorticity field on the right-hand side of Equation 16b, Equation 16c, this effect disappears at the small (i.e. |q|) and large (i.e. |q|0) scale, as long as both viscous and frictional dissipation are present.

Structure factor

The static structure factor can be expressed in integral form as

(28) S(q)=dω2πS(q,ω).

where the dynamic structure factor S(q,ω), can be calculated from the correlation function

(29) δρ^(q,ω)δρ^(q,ω)=(2π)3S(q,ω)δ(q+q)δ(ω+ω).

To compute the left-hand side of Equation 29 one can solve Equation 24 with respect to δρ^, δϑ^, and δφ^. This gives

(30a) δρ^=iη^(ρ)ωiM^ρρη^(ϑ)[M^ρϑ(ωiM^φφ)+iM^ρφM^φϑ]+η^(φ)[M^ρφ(ωiM^ϑϑ)+iM^ρϑM^ϑφ](ωiM^ρρ)[ω2iω(M^ϑϑ+M^φφ)M^ϑϑM^φφ+M^ϑφM^φϑ],
(30b) δϑ^=η^(ϑ)(iω+M^φφ)η^(φ)M^ϑφ[ω2iω(M^ϑϑ+M^φφ)M^ϑϑM^φφ+M^ϑφM^φϑ],
(30c) δφ^=η^(φ)(iω+M^ϑϑ)η^(ϑ)M^φϑ[ω2iω(M^ϑϑ+M^φφ)M^ϑϑM^φφ+M^ϑφM^φϑ].

The static structure factor can then be expressed as

(31) S=S(ρ)+S(ϑ)+S(φ).

The first term on the right-hand side can be readily calculated in the form

(32) S(ρ)=dωπH^(ρ)M^ρρ2+ω2=H^(ρ)|M^ρρ|=ρ0|q|2Ξ(v)cs2[(η+ζ)|q|2+ς],

indicating that, if driven solely by pressure fluctuations, the system would relax toward a structureless homogeneous state with Sρ0Ξ(ρ)/(ςcs2) when |q|0. The effect of the active currents is instead accounted for by the second and third terms on the right-hand side of Equation 31, which can be cast in the general form

(33) S(α)=H(α)dωπg(α)(ω)|h(ω)|2,α={ϑ,φ},


(34a) g(ϑ)(ω)=(M^ρϑω)2+(M^ρφM^φϑM^ρϑM^φφ)2,
(34b) g(φ)(ω)=(M^ρφω)2+(M^ρϑM^ϑφM^ρφM^ϑϑ)2,
(34c) h(ω)=(ωiM^ρρ)[ω2iω(M^ϑϑ+M^φφ)M^ϑϑM^φφ+M^ϑφM^φϑ].

The integral over ω can be derived using the residue theorem upon computing the roots of the complex third-order polynomial h. To make progress, we express

(35) |h(ω)|2=(ω2+ω12)(ω2+ω22)(ω2+ω32),

where ω1, ω2, and ω3 are given by

(36a) ω1=M^ρρ,
(36b) ω2=12(M^ϑϑ+M^φφ(M^ϑϑM^φφ)2+4M^ϑφM^φϑ),
(36c) ω3=12(M^ϑϑ+M^φφ+(M^ϑϑM^φφ)2+4M^ϑφM^φϑ.).

The integrand on the right-hand side of Equation 33 has, therefore, three pairs of purely imaginary poles: i.e. ±i|ω1|, ±i|ω2|, and ±i|ω3|. Next, turning the integration range to an infinite semicircular contour on the complex upper half-plane and summing the associated residues gives, after lengthy algebraic manipulations

(37a) S(ϑ)=H(ϑ)[Ω1M^ρϑ2+Ω2(M^ρφM^φϑM^ρϑM^φφ)2]Ω1Ω2Ω3Ω12,
(37b) S(φ)=H(φ)[Ω1M^ρφ2+Ω2(M^ρϑM^ϑφM^ρφM^ϑϑ)2]Ω1Ω2Ω3Ω12,

where we have set

(38a) Ω1=|ω1||ω2||ω3|,
(38b) Ω2=|ω1|+|ω2|+|ω3|,
(38c) Ω2=|ω1||ω2|+|ω1||ω3|+|ω2||ω3|.

Now, although the individual elements of the matrix M^ depend on the individual components of the wave vector – i.e. qx and qy – this is an artefact of linearizing the hydrodynamic equations about a specific orientation (i.e. ϑ=φ=0 in this case). Because of the lack of long-ranged order and of specific directions that could affect the spectrum of density fluctuations, the latter is expected to be isotropic, thus S=S(|q|). To remove the fictitious angular dependence, one can either linearize Equation 2a about a generic pair of angles, ϑ0 and φ0, and then use these to calculate a circular average – i.e. S(|q|)=1/(2π)2dϑ0dφ0S(q) – or, more simply, by orienting q so to cancel the directional dependence. Thus, taking qx=qy=|q|/2 gives a simpler expression of the matrix M^. That is

(39) M^=[ρ0cs2|q|2(η+ζ)|q|2+ςρ0α2|q|2(η+ζ)|q|2+ς3ρ0α6|q|64[(η+ζ)|q|2+ς]0D2|q|2K2|q|44(η|q|2+ς)+94χ2K6|q|44(η|q|2+ς)94χ20K2|q|44(η|q|2+ς)14χ6D6|q|2K6|q|44(η|q|2+ς)+14χ6].

Using the elements of this matrix in combination with Equations 31; 33, Equations 36a; 38a yields the curves plotted in Figure 3. Finally, asymptotically expanding Equation 31 allows one, after lengthy algebraic manipulations, to calculate the coefficients s2 and s4 in Equation 8. That is

(40a) s2=ρ0α22[(9χ2)2Ξφ+χ62Ξϑ]cs2(9χ2D6+χ6D2)[ρ0cs2(9χ2+χ6)+ς(9χ2D6+χ6D2)],
(40b) s4=72ρ0α62[(K22+8ηD2K2+8η2D22)Ξ(v)+K22Ξϑ+2η2(K2+4ηD2)2Ξφ]cs2(η+ζ)[K2+K6+4η(D2+D6)]4.

Notice that, while both orientational and translation noise affect the amplitude of density fluctuations at small length scales, where S(|q|)s4|q|4, translational noise becomes unimportant at the large scale, where S(|q|)s2/|q|2. Furthermore, as long as viscous dissipation is at play, switching off translational noise (i.e. Ξ(v)0) does not alter the scaling behavior of the structure factor at neither range of length scales. Taking the dry limit (i.e. η0 and ζ0) leaves the large-scale behavior unaltered, but does affect the scaling of density fluctuations at short length scales, where translational fluctuations are most prominent. Specifically, S(|q|)s6|q|6 in the case of purely rotational noise and S(|q|)s10|q|10 in the presence of rototranslational noise. The coefficients s6 and s10 can be computed as in the viscous case, to give

(41a) s6=(32)2ρ02α62Ξ(φ)ς2(D2+D6)3,
(41b) s10=(34)2ρ02α62Ξ(φ)ς4(D2+D6)3.

Numerical methods

The Voronoi model

In the self-propelled Voronoi model (Bi et al., 2016) a confluent cell layer is approximated as a Voronoi tessellation of the plane. Each cell is characterized by the position rc of its center, with c=1,2N, and a velocity vc=v0(cosθcex+sinθcey), with v0 a constant speed and θc an orientation. We stress that, in general, the center of a Voronoi polygon does not correspond to the polygon’s centroid (i.e. center of mass). The dynamics of these variables is governed by the following set of overdamped Langevin equations, expressing the interplay between cells’ autonomous motion and the remodeling events that underlie the tissue’s collective dynamics. That is:

(42a) drcdt=vcμrcE,
(42b) dθcdt=ηc,

where µ is the mobility coefficient and E=E(r1,r2rN) is an energy function involving exclusively geometrical quantities, such as the area Ac and the perimeter Pc of each cell: i.e.

(43) E=c[KA(AcA0)2+KP(PcP0)2],

with KA, KP, A0, and P0 constants. The first term in Equation 43 embodies a combination of cells’ volumetric incompressibility and monolayer resistance to thickness fluctuations. The second term results from the cytoskeletal contractility (quadratic in Pc) and the effective interfacial tension caused by the cell–cell adhesion and the cortical tension (both linear in Pc) (Farhadifar et al., 2007). The constants A0 and P0 represent, respectively, the preferred area and perimeter of each cell. The quantity ηc, on the other hand, is a random number with zero mean and correlation function

(44) ηc(t)ηc(t)=2Drδccδ(tt),

with Dr a rotational diffusion coefficient. To make progress, we next introduce the following dimensionless numbers: the shape index p0=P0/A0, which accounts for the spontaneous degree of acircularity of individual cells (Bi et al., 2016), and the Péclet number Pe=v0/(DrA0), which quantifies the persistence of directed cellular motion in front of their diffusivity.

To obtain the plots in Figure 3, we numerically integrate Equation 42a in a domain of size Lg with periodic boundary conditions. At t=0, the centroids rc are placed in a slightly perturbed hexagonal grid with a random initial velocity. After reaching the non-equilibrium steady state, we perform statistical averages of relevant observables.

In our numerical simulations, we set p0=3.85, μKAA0/Dr=1, μKP/Dr=1, and DrΔt=5×103, where Δt is the time-step used for the integration, and the average density of particles NA0/Lg2=1. We vary the Péclet number in the range 0.1Pe2.0. The results presented in Results are robust to the variation of the system size, as no qualitative difference was observed upon varying the domain size in the range 30Lg200 at constant density. The density structure factor (light green circles) in Figure 3a was obtained, in particular, with Pe=1.5.

The MPF model

The MPF model is a continuous model where each cell is described by a concentration field φc=φc(r) with c=1,2N and N the total number of cells. This model has been used to study the dynamics of confluent cell monolayers (Loewe et al., 2020) and the mechanics of cell extrusion (Monfared et al., 2021). Equilibrium configurations are obtained upon relaxing the free energy F=dAf, where the free energy density f is given by

(45) f=α4cφc2(φcφ0)2+kφ2c(φc)2+ϵc<cφc2φc2+cλ(11πϕ02Rφ2dAφc2)2.

Here, α and kϕ are material parameters which can be used to tune the surface tension γ=(8kφα/9)1/2 and the interfacial thickness ξ=(2kφ/α)1/2 of isolated cells and thermodynamically favor spherical cell shapes. The constant ϵ captures the repulsion between cells. The concentration field is large (i.e. φcϕ0) inside the cells and zero outside. The contribution proportional to λ in the free energy enforces cell incompressibility whose nominal radius is given by Rφ. The relaxational dynamics of the field φc is governed by the Allen–Cahn equation

(46) tφc+vcφc=MδFδφc,

where vc has the same meaning as in the SPV model described in the previous section and its dynamics is also governed by Equation 42b. The constant M in Equation 46 is the mobility measuring the relevance of thermodynamic relaxation with respect to non-equlibrium cell migration. The dimensionless parameters of the model are the Péclet number Pe=v0/(2DrRφ) and the cell deformability d=ϵ/α.

The system of partial differential equations, Equation 46, 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). The C-code implemented for numerical integration is parallelized by means of Message Passage Interface (MPI). We consider systems of N=361 cells in a square domain of Lg=380 grid points. Model parameters in simulation units are as follows: Rϕ=11, φ0=2.0, Mα=0.006, Mkφ=0.006, Mϵ=0.01, Mλ=600, Mγ=0.008, DrΔt=104, being Δt the time-step used to integrate Equation 46. We vary the speed of self-propulsion in the range 0.0v00.005. In terms of dimensionless parameters this corresponds to having d=1.66 and Pe ranging between 0 and 2.30. The timescale of cell motility with respect to the timescale of elastic relaxation driven by surface tension v0/(Mγ) ranges between 0 and 0.625. Moreover, the nominal packing fraction is N(πRφ2)/Lg2=0.95, while the ratio between the interface thickness and the nominal radius ξ/Rφ=0.12. The density structure factor (dark green triangles) in Figure 3a was obtained with Pe=1.38.

Numerical method for integration of the hydrodynamic equations

Equation 2a has been integrated by means of a hybrid lattice Boltzmann (LB) method, in which Equation (2b) is solved through a predictor–corrector LB algorithm and the remaining equations via a predictor–corrector finite-difference Euler approach, with a first-order upwind scheme and second-order accurate stencils for the computation of spacial derivatives (Carenza et al., 2019). The code has been parallelized by means of MPI, by dividing the computational domain in slices and by implementing the ghost-cell method to compute derivatives on the boundary of the computational subdomains. Runs have been performed using 64 CPUs in two-dimensional geometries, on a computational box of size 2562 and 5122, for at least 1.5×107 LB iterations (corresponding to ∼21 and ∼84 days of CPU-time, respectively, for the smaller and larger computational boxes). Periodic boundary conditions have been imposed. The director fields (for both p=2 and p=6) have been randomly initialized. The initial density field is assumed to be uniform with ρ=2.0 everywhere. The model parameters in simulations units are as follows: η=ζ=1.66, λ2=λ6=1.1, ν2=ν6=0.0, Γ2=0.4, A2=B2=0.04, L2=0.04, Γ6=0.4, A6=B6=0.004, L6=0.004, κ2,6=ξ2,6=0.004. Nematic activity α2 has been varied in the range 0.02α20.0005 and hexatic activity α6 in the range 0.050α60.050. We set the active parameters β2 and β6=0. The density structure factor (continuous black line) in Figure 3a was obtained with α2=2×103 and α6=2×102.

The coherence length of the nematic and hexatic liquid crystal can be expressed as the (Lp/Ap)1/2=ΔxLB for both p=2,6, where ΔxLB is the grid spacing of the LB algorithm. The active length scale as defined in the main text is given for the active nematics as 2 and ranges between 10ΔxLB for α2=0.0005 and 1.5ΔxLB for α2=0.02. Conversely, for hexatics 6 and ranges up to 3.5ΔxLB for |α6|=0.05. To compare the results of the hydrodynamics simulations with the discrete models in Figure 3a, we choose 2ΔxLB=A0 and 2ΔxLB=RφΔxMP, with ΔxMP the grid spacing used to integrate Equation 46.

Comparison with passive liquid crystals with coupled order parameters

In this section, we show how multiscale hexanematic order differs from previously reported examples of liquid crystal order with coupled order parameters (Bruinsma and Aeppli, 1982; Selinger and Nelson, 1989; Selinger, 1991). To quantify the interplay between nematic and hexatic order, here we focus on the function C26(r) given in Equation 9, reflecting the amount of cross-correlation in their fluctuations. Here, ψ2=e2iϑ and ψ6=e6iφ, while the fluctuating fields ϑ and φ represent again the local nematic and hexatic orientations, respectively. Averaging ψ2 and ψ6 over the scale of a volume element, yields the order complex parameters Ψ2=e2iϑ=|Ψ2|e2iθ and Ψ6=e6iφ=|Ψ6|e6iϕ, with θ and ϕ the average orientations. To make progress, we assume that, at the scale of a volume element, both microscopic orientations ϑ and φ are Gaussianly distributed about their mean values, so that, in general

(47) Ψp=ψpe12var[Arg(ψp)]+iArg(ψp),

from which

(48) |Ψp|e12var[Arg(ψp)],Arg(Ψp)=Arg(ψp).

This approximation holds when the relative fluctuation of the p-atic phase Arg(ψp) is sufficiently small, so that

(49) |Ψp|112[Arg(ψp)Arg(Ψp)]2cos[Arg(ψp)Arg(Ψp)],

consistent with the standard definition of p-atic order parameter. Thus, in particular, θ=ϑ and |Ψ2|=cos2(ϑθ), whereas ϕ=φ and |Ψ6|=cos6(φϕ). This allows to write C26(r), as given by Equation (9), in the form

(50) C26(r)=Ψ2(r)Ψ6(0)+Ψ2(r)Ψ6(0)2e12[ϑ(r)φ(0)ϑ(r)φ(0)].

At equilibrium, both nematic and hexatic order can be approximated as uniform, so that

(51) Ψ2(r)Ψ6(0)+Ψ2(r)Ψ6(0)2=|Ψ2||Ψ6|cos(2θ6ϕ)const,

and the problem reduces to calculating the connected correlation function

(52) Cϑφ(r)=ϑ(r)φ(0)ϑ(r)φ(0).

Notice that Equation (51) is not strictly valid for a quasi long-ranged ordered liquid crystal, where also θ and ϕ are expected to vary in space. These spatial variations, however, occur on length scales comparable with the system size and, as long as this is much larger than any of the intrinsic length scales entailed in Equation 2a, are negligible for the purpose of this calculation. To compute Cϑφ(r), one can take the passive limit of Equation 2c and linearize the resulting equations about the lowest free energy configuration. This, in turn, is determined by the sign of the constant χ2,6 in Equation 6b. For χ2,6<0, the hexatic and nematic directors are energetically favored to be parallel, so that ϑφ. Conversely, when χ2,6>0, the hexatic and nematic directors are preferentially tilted by π/6, hence ϑ=φ±π/6. For presentational clarity, here we focus on the former case and, at the end of this section, we show how the same behavior holds for positive χ2,6 values. Thus, assuming χ2,6<0 and expanding Equation 2c about ϑφ, gives

(53a) tϑ=D22ϑ94|χ2|(ϑφ)+ξ(ϑ),
(53b) tφ=D62φ14|χ6|(φϑ)+ξ(φ)

where, as in the previous sections, we have set Dp=ΓpLp and χp=Γpχ2,6 and introduced the Gaussian noise fields ξ(ϑ) and ξ(ϑ), having vanishing mean and finite variance. Unlike the active case, however, at equilibrium the latter is related to the environmental temperature by the fluctuation–dissipation theorem. This implies

(54) ξ(α)(r,t)ξ(β)(r,t)=2kBT(δαϑδβϑγ2+δαφδβφγ6)δ(rr)δ(tt),

where γp=Kp/Dp, with Kp the orientational stiffness defined in Equation 15, is the rotational viscosity of the associated p-atic phase. Equation 53a can now be decoupled and used to compute the correlation function Cϑφ(r). For simplicity, here we set D2=D6=D, γ2=γ6=γ, and 9χ2=χ6=2χ. With this choice, taking

(55a) φ+=12(φ+ϑ),
(55b) φ=12(φϑ),

gives, after simple algebraic manipulations

(56a) tφ+=D2φ++ξ+,
(56b) tφ=D2φ|χ|φ+ξ,

where ξ+=(ξ(φ)+ξ(ϑ))/2 and ξ=(ξ(φ)ξ(ϑ))/2. Moreover, using Equation (54), one finds

(57) ξn(r,t)ξm(r,t)=2kBTγδnmδ(rr)δ(tt),

where {n,m}={+,}. Equation 56a can now be solved in Fourier space and real time to give

(58) φ^n(q,t)=eSn(q,t)[φ^n(q,0)+0tdteSn(q,t)ξ^n(q,t)],

where the hat indicates Fourier transformation and

(59) Sn(q,t)=Dt(|q|2+mn2),

where m+=0 and m2=χ2=D/|χ|. The calculation of the cross-correlation function Cϑφ(r) is now reduced to calculating the autocorrelation functions of the fields φ+ and φ. Specifically

(60) Cϑφ(r)=C++(r)C(r),


(61) Cnm(r)=φn(r)φm(0)φn(r)φm(0),

and we have made use of Equation (54) to demonstrate that C+(r)=C+(r)=0. The non-vanishing correlation functions, on the other hand, can be expressed as

(62) Cnn(r)=limt0<|q|<Λd2q(2π)2eiqr|φ^n(q,t)|2,

where Λ=2π/a is a short-distance cut-off and |φ^n(q,t)|2 is the finite-time orientational structure factor defined from the relation

(63) φ^n(q,t)φ^n(q,t)=(2π)2|φ^n(q,t)|2δ(q+q)δ(tt).

After standard algebraic manipulations one finds

(64) |φ^n(q)|2=limt|φ^n(q,t)|2=kBTK1|q|2+mn2.

from which Equation (62) can be calculated in the form

(65) Cnn(r)=kBTK0<|q|<Λd2q(2π)2eiqr|q|2+mn2.

Evidently, Equation (65) is equivalent to that obtained in a purely static setting from the Hamiltonian

(66) H=12d2r[K|φ+|2+K|φ|2+m2φ2],

of the non-interacting scalar fields φ+ and φ. Now, in the case of the ‘massive’ field φ, the Fourier integral in Equation (65) converges to

(67) C(r)=kBT2πKK0(|r|χ),

in the range |r|a. Here, K0 is a modified Bessel function of the second kind, whose asymptotic expansion at short and long distances is given by

(68) K0(z){γEMlogz20<z1,π2zezz1,

with γEM the Euler–Mascheroni constant. In the case of the ‘massless‘ field φ+, on the other hand, the Fourier integral diverges in the infrared, but the correlation function C++(r) can still be computed as the Laplacian Green function on an infinite domain punctured by a hole of radius a at the origin. Thus

(69) C++(r)=kBT2πKlog|r|a.

Combining this with Equations (67) and (69) yields the following expression for the correlation function

(70) Cϑφ(r)=kBT2πK[log|r|a+K0(|r|χ)],

where |r|a. Finally, using Equation (50) and the asymptotic expansions of Equation (68) gives the following expression for the cross-correlation function

(71) C26(r){const.|r|χ(|r|a)η26|r|χ,

where η26 is an instance of the generic non-universal exponent

(72) ηpp=ppkBT2πK,

in the specific case p=2 and p=6. Lastly, when χ2,6>0, the same procedure can be carried out by expanding Equation (2c) about ϑ=φ±π/6 and taking φ+=(φ+ϑ)/2 and φ=(φϑ±π/6)/2, from which one finds again Equation 72.

Data availability

Figures 2–4 contain the numerical data used to generate the figures.


    1. Rauzi M
    (2020) Cell intercalation in a simple epithelium
    Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences 375:20190552.

Article and author information

Author details

  1. Josep-Maria Armengol-Collado

    Instituut-Lorentz, Leiden University, Leiden, Netherlands
    Conceptualization, Data curation, Formal analysis, Visualization, Writing - original draft
    Contributed equally with
    Livio Nicola Carenza
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-0740-3040
  2. Livio Nicola Carenza

    Instituut-Lorentz, Leiden University, Leiden, Netherlands
    Conceptualization, Data curation, Formal analysis, Supervision, Visualization, Writing - original draft
    Contributed equally with
    Josep-Maria Armengol-Collado
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-5996-331X
  3. Luca Giomi

    Instituut-Lorentz, Leiden University, Leiden, Netherlands
    Conceptualization, Formal analysis, Supervision, Methodology, Writing - original draft, Project administration
    For correspondence
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-7740-5960


European Research Council

  • Livio Nicola Carenza
  • Luca Giomi

Nederlandse Organisatie voor Wetenschappelijk Onderzoek

  • Josep-Maria Armengol-Collado

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


We are indebted with Massimo Pica Ciamarra and David Nelson for insightful discussions. This work is supported by the ERC-CoG grant HexaTissue (LNC and LG) and by Netherlands Organization for Scientific Research (NWO/OCW) as part of the research program 'The active matter physics of collective metastasis' with project number Science-XL 2019.022 (J-M A-C). Part of this work was carried out on the Dutch national e-infrastructure with the support of SURF through the NWO Grant 2021.028 for computational time.

Version history

  1. Preprint posted: February 1, 2022 (view preprint)
  2. Received: January 24, 2023
  3. Accepted: January 5, 2024
  4. Accepted Manuscript published: January 8, 2024 (version 1)
  5. Version of Record published: March 25, 2024 (version 2)


© 2024, Armengol-Collado, Carenza et al.

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.


  • 648
  • 175
  • 1

Views, downloads and citations are aggregated across all versions of this paper published by eLife.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Josep-Maria Armengol-Collado
  2. Livio Nicola Carenza
  3. Luca Giomi
Hydrodynamics and multiscale order in confluent epithelia
eLife 13:e86400.

Share this article


Further reading

    1. Cell Biology
    2. Physics of Living Systems
    Ivan Castello-Serrano, Frederick A Heberle ... Ilya Levental
    Research Article

    The organelles of eukaryotic cells maintain distinct protein and lipid compositions required for their specific functions. The mechanisms by which many of these components are sorted to their specific locations remain unknown. While some motifs mediating subcellular protein localization have been identified, many membrane proteins and most membrane lipids lack known sorting determinants. A putative mechanism for sorting of membrane components is based on membrane domains known as lipid rafts, which are laterally segregated nanoscopic assemblies of specific lipids and proteins. To assess the role of such domains in the secretory pathway, we applied a robust tool for synchronized secretory protein traffic (RUSH, Retention Using Selective Hooks) to protein constructs with defined affinity for raft phases. These constructs consist solely of single-pass transmembrane domains (TMDs) and, lacking other sorting determinants, constitute probes for membrane domain-mediated trafficking. We find that while raft affinity can be sufficient for steady-state PM localization, it is not sufficient for rapid exit from the endoplasmic reticulum (ER), which is instead mediated by a short cytosolic peptide motif. In contrast, we find that Golgi exit kinetics are highly dependent on raft affinity, with raft preferring probes exiting the Golgi ~2.5-fold faster than probes with minimal raft affinity. We rationalize these observations with a kinetic model of secretory trafficking, wherein Golgi export can be facilitated by protein association with raft domains. These observations support a role for raft-like membrane domains in the secretory pathway and establish an experimental paradigm for dissecting its underlying machinery.

    1. Microbiology and Infectious Disease
    2. Physics of Living Systems
    Chi Zhang, Rongjing Zhang, Junhua Yuan
    Research Article

    Bacteria in biofilms secrete potassium ions to attract free swimming cells. However, the basis of chemotaxis to potassium remains poorly understood. Here, using a microfluidic device, we found that Escherichia coli can rapidly accumulate in regions of high potassium concentration on the order of millimoles. Using a bead assay, we measured the dynamic response of individual flagellar motors to stepwise changes in potassium concentration, finding that the response resulted from the chemotaxis signaling pathway. To characterize the chemotactic response to potassium, we measured the dose–response curve and adaptation kinetics via an Förster resonance energy transfer (FRET) assay, finding that the chemotaxis pathway exhibited a sensitive response and fast adaptation to potassium. We further found that the two major chemoreceptors Tar and Tsr respond differently to potassium. Tar receptors exhibit a biphasic response, whereas Tsr receptors respond to potassium as an attractant. These different responses were consistent with the responses of the two receptors to intracellular pH changes. The sensitive response and fast adaptation allow bacteria to sense and localize small changes in potassium concentration. The differential responses of Tar and Tsr receptors to potassium suggest that cells at different growth stages respond differently to potassium and may have different requirements for potassium.