Hydrodynamics and multiscale order in confluent epithelia

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.


,
(S1) with ϕ v = Arg(r v ) the angle between r v and the x−axis of a Cartesian frame.A schematic representation of these elements in an arbitrary irregular polygon is shown in Fig. S1a.Unlike the complex function ψ p = e ipϑ , 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  Figs.S1b and S1c shows examples of the functions γ 2 and γ 6 for a typical configuration of the SPV.We emphasize that γ p , which, as show in Ref. [20], arises from a p−fold generalization of the classic shape tensor [34], is solely determined by the positions of the vertices of an individual polygon and, therefore, does not depend on the spatial organization of the neighbouring 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: where the r c is the position of the c−th cell, Θ is the Heaviside step function, such that Θ(x) = 1 for x > 0 and 0 otherwise, and is the number of cells within a distance ℓ from r c .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 Fig. S2 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 θ e x + sin θ e y .

PASSIVE STRESSES
As explained in the main text, the passive contribution to the stress tensor is given by σ (p) = −P 1+σ (e) +σ (r) +σ (v) , where, as demonstrated in [31,32] 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. |Q p | 2 = |Ψ p | 2 /2 = const, and taking λ p = 0, Eq. (S4) reduces to the expression derived in [29,30].That is where the first term in Eq. ( S4b) has incorporated into the pressure P and K p denotes the orientational stiffness of the p−atic phase, related to the order parameter stiffness by and ε is the two-dimensional antisymmetric tensor, with ε xy = −ε yx = 1 and ε xx = ε yy = 0.

LINEAR FLUCTUATING HYDRODYNAMICS
To compute the structure factor, we follow [48] and augment Eqs. ( 2b) and (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 where δϑ, δφ and δv indicate a small departure from the homogeneous and stationary configurations of the fields ϑ, φ and v, D p = Γ p L p , χ p = Γ p χ 2,6 , and ξ (ϑ) and ξ (φ) are short-ranged correlated noise fields: i.e.
The velocity field δv, on the other hand, is found from the Stokes limit of Eq. (2b) in the main text, which, at the linear order in all fluctuating fields, takes the form where 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 with {i, j} ∈ {x, y} and the case of noiseless translational dynamics, corresponding to Fig.
(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 with c s the speed of sound.Together with the expression for the active stress given in Eq. (3) of the main text, this gives Now, in Fourier space Eq.(S9) can be cast in the form of the following linear algebraic equation where the hat denotes Fourier transformation.Next, using and solving Eq. ( S13) and incorporating the resulting velocity field in Eqs.(S7) gives, after several algebraic manipulation where the matrix M is given by , and the functions η (α) , with α ∈ {ρ, ϑ, φ}, are effective noise fields whose correlation functions are given by η where the functions Ĥ(α) = Ĥ(α) (q) are given by 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 Eqs. ( S7b) and (S7c), 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 where the dynamic structure factor S(q, ω), can be calculated from the correlation function To compute the left-hand side of Eq. (S20) one can solve Eq. (S15) with respect to δ ρ, δ θ and δ φ.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 S → ρ 0 Ξ (ρ) /(ςc 2 s ) when |q| → 0. The effect of the active currents is instead accounted for by the second and third term on the right-hand side of Eq. (S22), which can be cast in the general form where 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 where ω 1 , ω 2 and ω 3 are given by The integrand on the right-hand side of Eq. (S24) 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 where we have set ) ) Now, although the individual elements of the matrix M depend on the individual components of the wave vector -i.e.q x and q y -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 Eqs.
(2) about a generic pair of angles, ϑ 0 and φ 0 , and then use these to calculate a circular average -i.e.S(|q|) = 1/(2π) 2 dϑ 0 dφ 0 S(q) -or, more simply, by orienting q so to cancel the directional dependence.Thus, taking q x = q y = |q|/ √ 2 gives a simpler expression of the matrix M .That is Using the elements of this matrix in combination with Eqs.(S22), (S24), (S27) and (S29) yields the curves plotted in Fig.
(3).Finally, asymptotically expanding Eq. (S22) allows one, after lengthy algebraic manipulations, to calculate the coefficients s −2 and s 4 in Eq. ( 8).That is Notice that, while both orientational and translation noise affect the amplitude of density fluctuations at small length scales, where S(|q|) ∼ s 4 |q| 4 , translational noise becomes unimportant at the large scale, where S(|q|) ∼ s −2 /|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|) ∼ s 6 |q| 6 in the case of purely rotational noise and S(|q|) ∼ s 10 |q| 10 in the presence of rototranslational noise.The coefficients s 6 and s 10 can be computed as in the viscous case, to give NUMERICAL METHODS

The Voronoi model
In the self-propelled Voronoi model (SVM) [8] a confluent cell layer is approximated as a Voronoi tessellation of the plane.Each cell is characterized by the position r c of its center, with c = 1 , 2 . . .N , and a velocity v c = v 0 (cos θ c e x + sin θ c e y ), with v 0 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 remodelling events that underlie the tissue's collective dynamics.That is: where µ is the mobility coefficient and E = E(r 1 , r 2 . . .r N ) is an energy function involving exclusively geometrical quantities, such as the area A c and the perimeter P c of each cell: i.e.
with K A , K P , A 0 and P 0 constants.The first term in Eq. (S34) embodies a combination of cells' volumetric incompressibility and monolayer resistance to thickness fluctuations.The second term results from the cytoskeletal contractility (quadratic in P c ) and the effective interfacial tension caused by the cell-cell adhesion and the cortical tension (both linear in P c ) [6].The constants A 0 and P 0 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 with D r a rotational diffusion coefficient.To make progress, we next introduce the following dimensionless numbers: the shape index p 0 = P 0 / √ A 0 , which accounts for the spontaneous degree of acircularity of individual cells [8], and the Péclet number Pe = v 0 /(D r √ A 0 ), which quantifies the persistence of directed cellular motion in front of their diffusivity.
To obtain the plots in Fig.
(3), we numerically integrate Eqs.(S33) in a domain of size L g with periodic boundary conditions.At t = 0, the centroids r c 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 p 0 = 3.85, µK A A 0 /D r = 1, µK P /D r = 1, and D r ∆t = 5 × 10 −3 , where ∆t is the time-step used for the integration, and the average density of particles N A 0 /L 2 g = 1.We vary the Péclet number in the range 0.1 ≤ Pe ≤ 2.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 30 ≤ L g ≤ 200 at constant density.The density structure factor (light green circles) in Fig. (3) a was obtained, in particular, with Pe = 1.5.

The Multiphase field model
The multiphase field (MPF) model is a continuous model where each cell is described by a concentration field φ c = φ c (r) with c = 1, 2 . . .N and N the total number of cells.This model has been used to study the dynamics of confluent cell monolayers [11] and the mechanics of cell extrusion [12].Equilibrium configurations are obtained upon relaxing the free energy F = dA f , where the free energy density f is given by 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 where v c has the same meaning as in the SPV model described in the previous section and its dynamics is also governed by Eq. ( S33b).The constant M in Eq. ( S37) 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 = v 0 /(2D r R φ ) and the cell deformability d = ϵ/α.
The system of partial differential equations, Eq. (S37), is solved with a finite-difference approach through a predictorcorrector finite difference Euler scheme implementing second order stencil for space derivatives [51].The C-code implemented for numerical integration is parallelized by means of MPI.We consider systems of N = 361 cells in a square domain of L g = 380 grid points.Model parameters in simulation units are as follows: R ϕ = 11, φ 0 = 2.0, M α = 0.006, M k φ = 0.006, M ϵ = 0.01, M λ = 600, M γ = 0.008, D r ∆t = 10 −4 , being ∆t the time-step used to integrate Eq. (S37).We vary the speed of self-propulsion in the range 0.0 ≤ v 0 ≤ 0.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 v 0 /(M γ) ranges between 0 and 0.625.Moreover, the nominal packing fraction is N (πR 2 φ )/L 2 g = 0.95, while the ratio between the interface thickness and the nominal radius ξ/R φ = 0.12.The density structure factor (dark green triangles) in Fig. (3)a was obtained with Pe = 1.38.

NUMERICAL METHOD FOR INTEGRATION OF THE HYDRODYNAMIC EQUATIONS
Eqs. [2] have been integrated by means of a hybrid Lattice Boltzmann (LB) method, in which Eq. ( 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 [51].The code has been parallelized by means of Message Passage Interface (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 256 2 and 512 2 , for at least 1.5 × 10 7 lattice Boltzmann iterations (corresponding to ∼ 21 days 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: The coherence length of the nematic and hexatic liquid crystal can be expressed as the (L p /A p ) 1/2 = ∆x LB for both p = 2, 6, where ∆x LB is the grid spacing of the lattice Boltzmann algorithm.The active lengthscale as defined in the main text is given for the active nematics as ℓ 2 and ranges between 10∆x LB for α 2 = −0.0005and 1.5∆x LB for α 2 = −0.02.Conversely, for hexatics ℓ 6 and ranges up to 3.5∆x LB for |α 6 | = 0.05.To compare the results of the hydrodynamics simulations with the discrete models in Fig. (3a), we choose 2∆x LB = √ A 0 and 2∆x LB = R φ ∆x MP , with ∆x MP the grid spacing used to integrate Eq. (S37).

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 [44,45,46].To quantify the interplay between nematic and hexatic order, here we focus on the function C 26 (r) given in Eq. [9], reflecting the amount of cross-correlation in their fluctuations.
Here ψ 2 = e 2iϑ and ψ 6 = e 6iφ , while the fluctuating fields ϑ and φ represents 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 = ⟨e 2iϑ ⟩ = |Ψ 2 |e 2iθ and Ψ 6 = ⟨e 6iφ ⟩ = |Ψ 6 |e 6iϕ , 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 ) This approximation holds when the relative fluctuation of the p−atic phase Arg(ψ p ) is sufficiently small, so that consistent with the standard definition of p−atic order parameter.Thus, in particular, θ = ⟨ϑ⟩ and |Ψ 2 | = ⟨cos 2(ϑ − θ)⟩, whereas ϕ = ⟨φ⟩ and |Ψ 6 | = ⟨cos 6(φ − ϕ)⟩.This allows to write C 26 (r), as given by Eq. [9], in the form At equilibrium, both nematic and hexatic order can be approximated as uniform, so that and the problem reduces to calculating the connected correlation function Notice that Eq. ( S42) 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 Eqs.
(2), are negligible for the purpose of this calculation.To compute C ϑφ (r), one can take the passive limit of Eqs.(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 Eq. (6b).For χ 2,6 < 0, the hexatic and nematic directors are energetically favored to be parallel, so that ϑ ≈ φ.

Fig
Fig. S1.(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 Eq. (S1) 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.
Fig. S2.(a) Coarse-grained nematic orientation θ obtained from averaging the local shape of cells over domains of size 30 ℓ cell , 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.