Two-dimensional p-atic liquid crystals are traditionally described in terms of the orientation field ${\displaystyle {\psi }_p=e^{ip\vartheta }}$, with ${\displaystyle \vartheta }$ 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, ${\displaystyle {\mathit{Q}}_p=Q_{i_1{i}_2\cdots i_p}{\mathit{e}}_{i_1}\otimes {\mathit{e}}_{i_2}\otimes \cdots \otimes {\mathit{e}}_{i_p}}$ with ${\displaystyle i_n=\{x,y\}}$ and ${\displaystyle n=1,2\dots p}$, constructed upon averaging the pth tensorial power of the local orientation ${\displaystyle \mathit{\nu }=\cos \vartheta {\mathit{e}}_x+\sin \vartheta {\mathit{e}}_y}$. That is

where ${\displaystyle ⟨\cdots ⟩}$ denotes the ensemble average and the operator ${\displaystyle \Vert ...\Vert }$ has the effect of rendering an arbitrary tensor traceless and symmetric (Hess, 2015). The vector ${\displaystyle \mathit{n}=\cos \theta {\mathit{e}}_x+\sin \theta {\mathit{e}}_y}$ 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 ${\displaystyle 2\pi /p}$. The fields ${\displaystyle |{\mathrm{\Psi }}_p|}$ and ${\displaystyle \theta }$ represent, respectively, the magnitude and phase of the complex p-atic order parameter ${\displaystyle {\mathrm{\Psi }}_p=⟨{\psi }_p⟩}$, while the normalization factor is chosen so that ${\displaystyle |{\mathit{Q}}_p{|}^2=|{\mathrm{\Psi }}_p{|}^2/2}$ for all ${\displaystyle p}$ values. For ${\displaystyle p=2}$, Equation (1) readily gives the standard nematic order parameter tensor: i.e. ${\displaystyle {\mathit{Q}}_2=|{\mathrm{\Psi }}_2|(\mathit{n}\otimes \mathit{n}-\mathbb{1}/2)}$, with 1 the identity tensor. In practice, if a cell’s planar projection consists of a regular p-sided polygon, the microscopic orientation ${\displaystyle \vartheta }$ equates that of any of the vertices of the polygon. In the more realistic case of an irregular polygon, on the other hand, ${\displaystyle \vartheta }$ is given by the phase of the complex function ${\displaystyle {\gamma }_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 ${\mathit{Q}}_p$, the mass density ${\displaystyle \rho }$, and the momentum density ${\displaystyle \rho \mathit{v}}$, with ${\displaystyle \mathit{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):

where ${\displaystyle D/Dt={\mathrm{\partial }}_t+\mathit{v}\cdot \mathrm{\nabla }}$. Equation (2a) and Equation 2b are the mass and momentum conservation equations, with ${\displaystyle {\displaystyle k_{\mathrm{d}}}}$ and ${\displaystyle k_{\mathrm{a}}}$ rates of cell division and apoptosis, $\mathit{\sigma }$ the stress tensor and ${\displaystyle \mathit{f}}$ an arbitrary external force per unit area. In Equation (2c), ${\displaystyle {\mathrm{\Gamma }}_p^{-1}}$ is a rotational viscosity and ${\displaystyle {\mathit{H}}_p=-\delta F/\delta {\mathit{Q}}_p}$ is the molecular tensor describing the relaxation of the p-atic phase toward the minimum of the free energy ${\displaystyle F}$ (see Methods). The rank-2 tensors ${\displaystyle \mathit{\omega }=[\mathrm{\nabla }\mathit{v}-(\mathrm{\nabla }\mathit{v}{)}^{⊺}]/2}$ and ${\displaystyle \mathit{u}=[\mathrm{\nabla }\mathit{v}+(\mathrm{\nabla }\mathit{v}{)}^{⊺}]/2}$, with ${\displaystyle {\displaystyle ⊺}}$ 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 ${\mathit{Q}}_p$ with one of $\mathit{\omega }$: i.e. ${({\mathit{Q}}_p\cdot \mathit{\omega })}_{i_1{i}_2\mathrm{\cdots }{i}_p}=Q_{i_i{i}_2\mathrm{\cdots }j}{\omega }_{j{i}_p}$. On the second line ${\displaystyle ({\mathrm{\nabla }}^{\otimes n}{)}_{i_1{i}_2\cdots i_n}={\mathrm{\partial }}_{i_1}{\mathrm{\partial }}_{i_2}\cdots {\mathrm{\partial }}_{i_n}}$, while ${\displaystyle \lfloor \dots \rfloor }$ denotes the floor function and ${\displaystyle p\mathrm{mod}2=p-2\lfloor p/2\rfloor }$ is zero for even ${\displaystyle p}$ values and one for odd ${\displaystyle p}$ values. Finally, ${\displaystyle {\overline{\lambda }}_p}$, ${\displaystyle {\lambda }_p}$, and ${\displaystyle {\nu }_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 ${\displaystyle \mathit{f}}$ in Equation 2b and the stress tensor $\mathit{\sigma }$. 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: ${\displaystyle \mathit{f}=-\varsigma \mathit{v}}$, with ${\displaystyle {\displaystyle \varsigma }}$ 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. $\mathit{\sigma }={\mathit{\sigma }}^{(\mathrm{p})}+{\mathit{\sigma }}^{(\mathrm{a})}$. The passive stress tensor is in turn expressed as ${\displaystyle {\mathit{\sigma }}^{(\mathrm{p})}=-P\mathbb{1}+{\mathit{\sigma }}^{(\mathrm{e})}+{\mathit{\sigma }}^{(\mathrm{r})}+{\mathit{\sigma }}^{(\mathrm{v})}}$, where ${\displaystyle P}$ is the pressure, ${\mathit{\sigma }}^{(\mathrm{e})}$ is the elastic stress, arising in response of a static deformation of a fluid patch, and ${\mathit{\sigma }}^{(\mathrm{r})}$ and ${\mathit{\sigma }}^{(\mathrm{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 ${\mathit{\sigma }}^{(\mathrm{p})}$ was derived in Giomi et al., 2022b and is reported in Methods.

The active stress ${\mathit{\sigma }}^{(\mathrm{a})}$, on the other hand, can be constructed phenomenologically for arbitrary ${\displaystyle p}$ values in the form

where the symbol ${\displaystyle \odot }$ 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 ${\displaystyle {\mathit{A}}_p}$ and ${\displaystyle {\mathit{B}}_q}$ be two generic tensors of rank ${\displaystyle {\displaystyle p<q}}$, then ${\displaystyle ({\mathit{A}}_p\odot {\mathit{B}}_q{)}_{i_1{i}_2\cdots i_{q-p}}=A_{j_1{j}_2\cdots j_p}{B}_{j_1{j}_2\cdots j_p{i}_1{i}_2\cdots i_{q-p}}}$. The sum over ${\displaystyle 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. ${\displaystyle {\mathit{F}}_c=\sum _{k=1}^p{\mathit{F}}_k\delta (\mathit{r}-{\mathit{r}}_c-a{\mathit{\nu }}_k)}$ with ${\displaystyle {\mathit{F}}_k}$ the force exerted by a cell at each vertex and originating from the imbalance of the tensions ${\displaystyle {\mathit{T}}_{kl}}$, driven by the active contraction of the cellular junctions, converging at the kth vertex: i.e. ${\displaystyle {\mathit{F}}_k=\sum _l{\mathit{T}}_{kl}}$ (see Figure 1c). The quantities ${\displaystyle {\mathit{r}}_c}$ and ${\displaystyle a}$ are the cell’s centroid and circumradius, respectively, while ${\displaystyle {\mathit{\nu }}_k=\cos (\vartheta +2\pi k/p){\mathit{e}}_x+\sin (\vartheta +2\pi k/p){\mathit{e}}_y}$. We stress that, while the individual tensions acting along the junctions are exclusively contractile, the resulting vertex forces can be either contractile (i.e. ${\displaystyle {\mathit{F}}_k\cdot {\mathit{\nu }}_k<0}$) or extensile (${\displaystyle {\mathit{F}}_k\cdot {\mathit{\nu }}_k>0}$), depending on the overall tension distribution and the geometry of the cellular network. Next, assuming ${\displaystyle {\mathit{F}}_k=f{\mathit{\nu }}_k}$ and expanding the delta function about ${\displaystyle a=0}$ yields ${\displaystyle {\mathit{F}}_c=\sum _{m=0}^{\mathrm{\infty }}{\mathit{f}}_m}$, where

Because of the p-fold symmetry of the force complexion ${\displaystyle {\mathit{f}}_m=\mathbf{0}}$ for all even ${\displaystyle {\displaystyle m}}$ values, unless ${\displaystyle {\displaystyle m=p-1}}$, whereas odd ${\displaystyle m}$ values yields, up to symmetrization, ${\displaystyle \sum _{k=1}^p{\mathit{\nu }}_k^{\otimes (m+1)}\sim {\mathbb{1}}^{\otimes (m+1)/2}}$. Thus, after some algebraic manipulation, one finds ${\displaystyle {\mathit{F}}_c\approx -apf/2\mathrm{\nabla }[(1+a^2/8{\mathrm{\nabla }}^2+\cdots )\delta (\mathit{r}-{\mathit{r}}_c)]+{\mathit{f}}_{p-1}}$. Finally, taking ${\displaystyle ⟨\sum _c{\mathit{F}}_c⟩=-P^{(\mathrm{a})}\mathbb{1}+{\mathit{\sigma }}^{(\mathrm{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

where ${\displaystyle n=⟨\sum _c\delta (\mathit{r}-{\mathit{r}}_c)⟩}$ is the cell number density. From Equation 5b, one finds the following expression for the phenomenological parameter ${\displaystyle {\displaystyle {\alpha }_p}}$ in Equation 3: i.e. ${\displaystyle {\alpha }_p=(-a{)}^{p-1}pnf/[\sqrt{2^{p-2}}(p-1)!]}$. Notice that both constants ${\displaystyle a}$ and ${\displaystyle 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 ${\mathit{\sigma }}^{(\mathrm{p})}$, cannot be derived from equilibrium considerations. Other terms constructed by contracting ${\mathit{Q}}_p$ with ${\displaystyle {\mathrm{\nabla }}^{\otimes 2}}$ can be expressed as linear combinations of this and ${\mathit{\sigma }}^{(\mathrm{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 ${\displaystyle Q_{k_1{k}_2\dots k_p}{\mathrm{\partial }}_i{\mathrm{\partial }}_j{Q}_{k_1{k}_2\cdots k_p}}$, ${\displaystyle Q_{ij{k}_3\cdots k_p}{\mathrm{\partial }}_{l_1}{\mathrm{\partial }}_{l_2}{Q}_{l_1{l}_2{k}_3\cdots k_p}}$, etc., are mechanically equivalent to the terms ${\displaystyle {\mathrm{\partial }}_i{Q}_{k_1{k}_2\cdots k_p}{\mathrm{\partial }}_j{Q}_{k_1{k}_2\cdots k_p}}$ and ${\displaystyle Q_{k_1{k}_2\cdots i}{H}_{k_1{k}_2\cdots j}-H_{k_1{k}_2\cdots i}{Q}_{k_1{k}_2\cdots j}}$ 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 ${\displaystyle p}$ value is considered. Such a crossover is expected at length scales larger than ${\displaystyle \ell =({\alpha }_p/{\beta }_p{)}^{1/(p-4)}}$, 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 ${\mathit{Q}}_2$ and ${\mathit{Q}}_6$. 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 ${\displaystyle F=\int \mathrm{d}A(f_2+f_6+f_{2,6})}$, where

Here, ${\displaystyle A_p}$ and ${\displaystyle B_p}$ 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 ${\displaystyle {\displaystyle {\kappa }_{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 ${\displaystyle {\chi }_{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 ${\displaystyle f_{2,6}}$ can further be augmented with several additional terms of higher differential order: e.g. $({\mathit{Q}}_2\odot \nabla {\mathit{Q}}_2)\cdot ({\mathit{Q}}_6\odot \nabla {\mathit{Q}}_6)$, ${|\nabla ({\mathit{Q}}_2^{\otimes 3}\odot {\mathit{Q}}_6)|}^2$, ${\nabla }^2({\mathit{Q}}_2^{\otimes 3}\odot {\mathit{Q}}_6)$, 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:

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, ${\displaystyle {\ell }_2}$ and ${\displaystyle {\ell }_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 ${\displaystyle {\ell }_2\approx {\ell }_6}$. 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 ${\displaystyle {\ell }_{\chi ,2}=\sqrt{L_2/{\chi }_{2,6}}}$ and ${\displaystyle {\ell }_{\chi ,6}=\sqrt{L_6/{\chi }_{2,6}}}$. Their role will be discussed in the following section, in the framework of fluctuating hydrodynamics.

Finally, in the passive limit, when ${\displaystyle {\alpha }_2=0}$ and ${\displaystyle {\alpha }_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 ${\displaystyle \ell \ll {\ell }_6}$ and nematic order at length scales ${\displaystyle \ell \gg {\ell }_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. ${\displaystyle {\ell }_6\ll \ell \ll {\ell }_2}$, 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.

To elucidate the multiscale organization of the system, we next compute the structure factor ${\displaystyle S(|\mathit{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 ${\displaystyle k_{\mathrm{d}}=k_{\mathrm{a}}}$ and ${\displaystyle {\lambda }_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⁻⁷–10⁻⁶, we neglect inertial effects: i.e. ${\displaystyle \rho D\mathit{v}/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 ${\displaystyle \mathit{v}=\mathbf{0}}$, and augmented with noise terms to give the following exact asymptotic expansion

The first term entails the typical giant number density fluctuations associated with the active nematic behavior at the large scale, with ${\displaystyle s_{-2}\sim {\alpha }_2^2}$. This effect is overestimated at the linear order, leading to an inverse quadratic dependence on the wave number ${\displaystyle |\mathit{q}|}$ (Ramaswamy et al., 2003), but is generally renormalized by nonlinearities, so that ${\displaystyle \underset{|\mathit{q}|\to 0}{lim}S(|\mathit{q}|)\sim |\mathit{q}{|}^{-\alpha }}$, with ${\displaystyle 1<\alpha <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 ${\displaystyle s_{\beta }\sim {\alpha }_6^2}$ and the exponent ${\displaystyle \beta }$ 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, ${\displaystyle \beta =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, ${\displaystyle \beta }$ differs depending on whether noise affects both cells’ orientational and translational dynamics, or only the former. Specifically, when only orientational noise is considered, ${\displaystyle \beta =6}$. In contrast, ${\displaystyle \beta =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.

Figure 2

Download asset Open asset

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 ${\displaystyle v_0}$ and whose direction of motion undergoes rotational diffusion with diffusion coefficient ${\displaystyle D_{\mathrm{r}}}$ (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 ${\displaystyle S(|\mathit{q}|)\sim |\mathit{q}{|}^{\beta }}$, with ${\displaystyle \beta }$ monotonically decreasing from 6 to 4 upon increasing the Péclet number ${\displaystyle \mathrm{P}\mathrm{e}={\xi }_0/a}$ expressing the ratio between cells’ persistence length ${\displaystyle {\xi }_0=v_0/D_{\mathrm{r}}}$ and their typical size ${\displaystyle a}$ (see Figure 3b).

Figure 3

Download asset Open asset

Figure 3—source data 1

https://cdn.elifesciences.org/articles/86400/elife-86400-fig3-data1-v2.xlsx

Download elife-86400-fig3-data1-v2.xlsx

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 ${\displaystyle \mathrm{P}\mathrm{e}\to 0}$, 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, ${\displaystyle \beta =6}$. Increasing ${\displaystyle \mathrm{P}\mathrm{e}}$ 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 ${\displaystyle \beta =6}$ to ${\displaystyle \beta =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 ${\displaystyle \mathrm{P}\mathrm{e}}$ 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

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 ${\displaystyle |\mathit{r}|}$ compares to the length scales ${\displaystyle {\ell }_{\chi ,2}}$ and ${\displaystyle {\ell }_{\chi ,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 ${\displaystyle {\ell }_{\chi ,2}={\ell }_{\chi ,6}={\ell }_{\chi }}$, fluctuations dominate at short distances and the hexatic and nematic orientations are uncorrelated. Thus, ${\displaystyle C_{26}(\mathit{r})}$ is approximatively constant for ${\displaystyle |\mathit{r}|\ll {\ell }_{\chi }}$. The picture is reversed for ${\displaystyle |\mathit{r}|\gg {\ell }_{\chi }}$. In this range, the hexatic and nematic orientations are ‘locked’ in a parallel configuration, i.e. ${\displaystyle \mathrm{Arg}({\psi }_2)/2\approx \mathrm{Arg}({\psi }_6)/6}$, or tilted by ${\displaystyle \pi /6}$ with respect to each other, depending on the sign of the constant ${\displaystyle {\chi }_{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. ${\displaystyle C_{26}(\mathit{r})\sim (|\mathit{r}|/{\ell }_{\chi }{)}^{-{\eta }_{26}}}$, with ${\displaystyle {\eta }_{26}}$ a specific instance of the generic non-universal exponent ${\displaystyle {\eta }_{26}=6k_{B}T/(\pi K)}$, with ${\displaystyle K}$ the orientational stiffness of both phases (proportional to ${\displaystyle L_2=L_6}$). An analytical treatment of this simple case is reported in Methods. In the more generic case, in which ${\displaystyle {\ell }_{\chi ,2}\ne {\ell }_{\chi ,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 ${\displaystyle {\alpha }_2=0}$ and ${\displaystyle {\alpha }_6=0}$, is shown in Figure 4a. The curves in Figure 4b correspond instead to simulated configurations of the cross-correlation function of ${\displaystyle C_{26}(\mathit{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.

Figure 4

Download asset Open asset

Figure 4—source data 1

https://cdn.elifesciences.org/articles/86400/elife-86400-fig4-data1-v2.xlsx

Download elife-86400-fig4-data1-v2.xlsx

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. ${\displaystyle p=0}$), polar (i.e. ${\displaystyle p=1}$), or nematic (i.e. ${\displaystyle 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. ${\displaystyle 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.

Following Armengol-Collado et al., 2023, we use the shape function${\displaystyle {\gamma }_p}$ to quantify the amount of p-fold symmetry of an arbitrary cell. Denoting ${\displaystyle {\mathit{r}}_v}$ with ${\displaystyle v=1,2\dots V}$, the positions of its vertices with respect to the cell’s center of mass, one has

with ${\displaystyle {\varphi }_v=\mathrm{Arg}({\mathit{r}}_v)}$ the angle between ${\displaystyle {\mathit{r}}_v}$ 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 ${\displaystyle {\psi }_p=e^{ip\vartheta }}$, which has unit magnitude by construction, the magnitude ${\displaystyle |{\gamma }_p|}$ quantify the resemblance of a generic polygon with a regular p-sided polygon of the same size, while the phase ${\displaystyle \vartheta =\mathrm{Arg}({\gamma }_p)/p}$ marks the orientation of the polygon. For regular V-sided polygons, ${\displaystyle |{\gamma }_p|=1}$ provided ${\displaystyle p}$ is an integer multiple of ${\displaystyle V}$ and ${\displaystyle |{\gamma }_p|\approx 0}$ otherwise. Furthermore, from ${\displaystyle {\gamma }_p}$ one can readily compute.

Figure 5

Download asset Open asset

Figure 5b, c shows examples of the functions ${\displaystyle {\gamma }_2}$ and ${\displaystyle {\gamma }_6}$ for a typical configuration of the SPV. We emphasize that ${\displaystyle {\gamma }_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 ${\displaystyle {\gamma }_p}$ can then be coarse grained at the length scale ${\displaystyle \ell }$ to construct the shape parameter:

where the ${\displaystyle {\mathit{r}}_c}$ is the position of the cth cell, ${\displaystyle \mathrm{\Theta }}$ is the Heaviside step function, such that ${\displaystyle \mathrm{\Theta }(x)=1}$ for ${\displaystyle x>0}$ and 0 otherwise, and ${\displaystyle N_{\ell }=\sum _c\mathrm{\Theta }(\ell -|\mathit{r}-{\mathit{r}}_c|)}$ is the number of cells within a distance ${\displaystyle \ell }$ from ${\displaystyle {\mathit{r}}_c}$. As in the case of ${\displaystyle {\gamma }_p}$, the magnitude of ${\displaystyle {\mathrm{\Gamma }}_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 ${\displaystyle p=2}$. The different patches in panel (a) are regions with uniform ${\displaystyle \theta =\mathrm{Arg}({\mathrm{\Gamma }}_2)/2}$, while in panel (b), there are plotted streamlines showing the orientation of the director ${\displaystyle \mathit{n}=\cos \theta {\mathit{e}}_x+\sin \theta {\mathit{e}}_y}$.

Figure 6

Download asset Open asset

As explained in the main text, the passive contribution to the stress tensor is given by ${\displaystyle {\mathit{\sigma }}^{(\mathrm{p})}=-P\mathbb{1}+{\mathit{\sigma }}^{(\mathrm{e})}+{\mathit{\sigma }}^{(\mathrm{r})}+{\mathit{\sigma }}^{(\mathrm{v})}}$, where, as demonstrated in Giomi et al., 2022b

where ${\displaystyle \eta }$ and ${\displaystyle \zeta }$ 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. ${|{\mathit{Q}}_p|}^2={|{\mathrm{\Psi }}_p|}^2/2=\mathrm{const}$, and taking ${\displaystyle {\lambda }_p=0}$, Equation 13a reduces to the expression derived in Zippelius, 1980; Zippelius et al., 1980. That is

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

and $\mathit{varepsilon}$ is the two-dimensional antisymmetric tensor, with ${\displaystyle {\epsilon }_{xy}=-{\epsilon }_{yx}=1}$ and ${\displaystyle {\epsilon }_{xx}={\epsilon }_{yy}=0}$.

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 ${\displaystyle \vartheta }$ and ${\displaystyle \phi }$ the nematic and hexatic fluctuating orientation fields and linearizing the hydrodynamic equations about the homogeneous and stationary solutions, ${\displaystyle \vartheta =\phi =0}$ and ${\displaystyle \mathit{v}=\mathbf{0}}$, gives

where ${\displaystyle \delta \vartheta }$, ${\displaystyle \delta \phi }$, and ${\displaystyle \delta \mathit{v}}$ indicate a small departure from the homogeneous and stationary configurations of the fields ${\displaystyle {\displaystyle \vartheta }}$, ${\displaystyle {\displaystyle \phi }}$, and ${\displaystyle \mathit{v}}$, ${\displaystyle {\mathcal{D}}_p={\mathrm{\Gamma }}_p{L}_p}$, ${\displaystyle {\chi }_p={\mathrm{\Gamma }}_p{\chi }_{2,6}}$, and ${\displaystyle {\xi }^{(\vartheta )}}$ and ${\displaystyle {\xi }^{(\phi )}}$ are short-ranged correlated noise fields: i.e.

The velocity field ${\displaystyle \delta \mathit{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

where ${\displaystyle {\mathit{f}}^{(\mathrm{p})}=\mathrm{\nabla }\cdot {\mathit{\sigma }}^{(\mathrm{p})}}$ and ${\displaystyle {\mathit{f}}^{(\mathrm{a})}=\mathrm{\nabla }\cdot {\mathit{\sigma }}^{(\mathrm{a})}}$ are the body forces resulting from the passive and active stresses, respectively. The quantity ${\mathit{\xi }}^{(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 ${\mathit{\xi }}^{(v)}$ is either conservative or null, from which

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

with ${\displaystyle c_{\mathrm{s}}}$ the speed of sound. Together with the expression for the active stress given in Equation 3 of the main text, this gives

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

where the hat denotes Fourier transformation. Next, using

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

where the matrix ${\displaystyle \hat{\mathit{M}}}$ is given by

and the functions ${\displaystyle {\eta }^{(\alpha )}}$, with ${\displaystyle \alpha \in \{\rho ,\vartheta ,\phi \}}$, are effective noise fields whose correlation functions are given by

where the functions ${\displaystyle {\hat{H}}^{(\alpha )}={\hat{H}}^{(\alpha )}(\mathit{q})}$ are given by

Notice that, while hydrodynamic flow has the effect of coloring the orientational noise embodied in the stochastic fields ${\displaystyle {\xi }^{(\vartheta )}}$ and ${\displaystyle {\xi }^{(\phi )}}$, via the vorticity field on the right-hand side of Equation 16b, Equation 16c, this effect disappears at the small (i.e. ${\displaystyle |\mathit{q}|\to \mathrm{\infty }}$) and large (i.e. ${\displaystyle |\mathit{q}|\to 0}$) scale, as long as both viscous and frictional dissipation are present.

The static structure factor can be expressed in integral form as

where the dynamic structure factor ${\displaystyle S(\mathit{q},\omega )}$, can be calculated from the correlation function

To compute the left-hand side of Equation 29 one can solve Equation 24 with respect to ${\displaystyle \delta \hat{\rho }}$, ${\displaystyle \delta \hat{\vartheta }}$, and ${\displaystyle \delta \hat{\phi }}$. This gives

The static structure factor can then be expressed as

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

indicating that, if driven solely by pressure fluctuations, the system would relax toward a structureless homogeneous state with ${\displaystyle S\to {\rho }_0{\mathrm{\Xi }}^{(\rho )}/(\varsigma c_{\mathrm{s}}^2)}$ when ${\displaystyle |\mathit{q}|\to 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

where

The integral over ${\displaystyle \omega }$ can be derived using the residue theorem upon computing the roots of the complex third-order polynomial ${\displaystyle h}$. To make progress, we express

where ${\displaystyle {\omega }_1}$, ${\displaystyle {\omega }_2}$, and ${\displaystyle {\omega }_3}$ are given by

The integrand on the right-hand side of Equation 33 has, therefore, three pairs of purely imaginary poles: i.e. ${\displaystyle \pm i|{\omega }_1|}$, ${\displaystyle \pm i|{\omega }_2|}$, and ${\displaystyle \pm i|{\omega }_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

where we have set