Hydrodynamics and multiscale order in confluent epithelia
Abstract
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 patic 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 powerlaw 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 cellresolved models of epithelia – i.e. the selfpropelled 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.
https://doi.org/10.7554/eLife.86400.sa0Introduction
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 (SerraPicamal et al., 2012) or invading the extracellular matrix (Haeger et al., 2020), multicellular systems can move coherently within their microenvironment 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 lifethreatening 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), cellresolved 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.
Despite their conceptual simplicity, cellresolved models agree remarkably well with experimental data on confluent monolayers (Park et al., 2015; Atia et al., 2018). In particular they account for a solidtoliquid 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 modelepithelia are separated by an intermediate hexatic phase, in which the system exhibits the typical sixfold rotational symmetry of twodimensional 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 ArmengolCollado et al., 2023; Eckert et al., 2023). This remarkable example of physical organization in biological matter, referred to as multiscale hexanematic order in ArmengolCollado et al., 2023, is believed to complement the complex network or regulatory pathways available to individual cells to achieve multicellular organization and select specific scaledependent 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 patic rotational symmetry (hereafter patic 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érezGonzález et al., 2019); (2) models built around the socalled shape tensor (Ishihara et al., 2017; Czajkowski et al., 2018; Hernandez and Marchetti, 2021; Grossman and Joanny, 2022), i.e. a rank2 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 smallscale 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 pfold 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 nonequilibrium 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
Twodimensional patic liquid crystals are traditionally described in terms of the orientation field $\psi}_{p}={e}^{ip\vartheta$, with $\vartheta$ the local orientation of the pfold mesogens. A more general approach, proposed in Giomi et al., 2022a; Giomi et al., 2022b and especially suited for hydrodynamics, revolves instead around the rankp tensor order parameter, $\mathit{Q}}_{p}={Q}_{{i}_{1}{i}_{2}\cdots \phantom{\rule{thinmathspace}{0ex}}{i}_{p}}{\mathit{e}}_{{i}_{1}}\otimes {\mathit{e}}_{{i}_{2}}\otimes \cdots \otimes {\mathit{e}}_{{i}_{p}$ with $i}_{n}=\{x,y\$ and $n=1,\phantom{\rule{thinmathspace}{0ex}}2\dots \phantom{\rule{thinmathspace}{0ex}}p$, constructed upon averaging the pth tensorial power of the local orientation $\mathit{\nu}=\mathrm{cos}\vartheta \phantom{\rule{thinmathspace}{0ex}}{\mathit{e}}_{x}+\mathrm{sin}\vartheta \phantom{\rule{thinmathspace}{0ex}}{\mathit{e}}_{y}$. That is
where $\u27e8\cdots \u27e9$ denotes the ensemble average and the operator $\Vert ...\Vert$ has the effect of rendering an arbitrary tensor traceless and symmetric (Hess, 2015). The vector $\mathit{n}=\mathrm{cos}\theta \phantom{\rule{thinmathspace}{0ex}}{\mathit{e}}_{x}+\mathrm{sin}\theta \phantom{\rule{thinmathspace}{0ex}}{\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 $2\pi /p$. The fields ${\mathrm{\Psi}}_{p}$ and $\theta$ represent, respectively, the magnitude and phase of the complex patic order parameter ${\mathrm{\Psi}}_{p}=\u27e8{\psi}_{p}\u27e9$, while the normalization factor is chosen so that ${\mathit{Q}}_{p}{}^{2}={\mathrm{\Psi}}_{p}{}^{2}/2$ for all $p$ values. For $p=2$, Equation (1) readily gives the standard nematic order parameter tensor: i.e. ${\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 psided polygon, the microscopic orientation $\vartheta$ equates that of any of the vertices of the polygon. In the more realistic case of an irregular polygon, on the other hand, $\vartheta$ is given by the phase of the complex function $\gamma}_{p$, arising form the pfold generalization of the classic shape tensor (Aubouy et al., 2003). This function was introduced in ArmengolCollado et al., 2023 and is reviewed in Methods for sake of completeness.
The order parameter tensor ${\mathit{Q}}_{p}$, the mass density $\rho$, and the momentum density $\rho \mathit{v}$, with $\mathit{v}$ the local velocity field, comprise the set of hydrodynamic variables describing the dynamics of a generic patic fluid, which in turn is governed by the following set of partial differential equations (Giomi et al., 2022a; Giomi et al., 2022b):
where $D/Dt={\mathrm{\partial}}_{t}+\mathit{v}\cdot \mathrm{\nabla}$. Equation (2a) and Equation 2b are the mass and momentum conservation equations, with $\displaystyle {k}_{\mathrm{d}}$ and $k}_{\mathrm{a}$ rates of cell division and apoptosis, $\mathit{\sigma}$ the stress tensor and $\mathit{f}$ an arbitrary external force per unit area. In Equation (2c), $\mathrm{\Gamma}}_{p}^{1$ is a rotational viscosity and $\mathit{H}}_{p}=\delta F/\delta {\mathit{Q}}_{p$ is the molecular tensor describing the relaxation of the patic phase toward the minimum of the free energy $F$ (see Methods). The rank2 tensors $\mathit{\omega}=[\mathrm{\nabla}\mathit{v}(\mathrm{\nabla}\mathit{v}{)}^{\u22ba}]/2$ and $\mathit{u}=[\mathrm{\nabla}\mathit{v}+(\mathrm{\nabla}\mathit{v}{)}^{\u22ba}]/2$, with $\displaystyle \u22ba$ 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 $({\mathrm{\nabla}}^{\otimes n}{)}_{{i}_{1}{i}_{2}\cdots \phantom{\rule{thinmathspace}{0ex}}{i}_{n}}={\mathrm{\partial}}_{{i}_{1}}{\mathrm{\partial}}_{{i}_{2}}\cdots \phantom{\rule{thinmathspace}{0ex}}{\mathrm{\partial}}_{{i}_{n}}$, while $\lfloor \dots \rfloor$ denotes the floor function and $p\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{1em}{0ex}}\mathrm{mod}\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}2=p2\lfloor p/2\rfloor$ is zero for even $p$ values and one for odd $p$ values. Finally, $\overline{\lambda}}_{p$, $\lambda}_{p$, and $\nu}_{p$ are material parameters expressing the strength of the coupling between patic 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 $\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: $\mathit{f}=\varsigma \mathit{v}$, with $\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 $\mathit{\sigma}}^{(\mathrm{p})}=P\mathbb{1}+{\mathit{\sigma}}^{(\mathrm{e})}+{\mathit{\sigma}}^{(\mathrm{r})}+{\mathit{\sigma}}^{(\mathrm{v})$, where $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 patic 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 $p$ values in the form
where the symbol $\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 $\mathit{A}}_{p$ and $\mathit{B}}_{q$ be two generic tensors of rank $\displaystyle p<q$, then $({\mathit{A}}_{p}\odot {\mathit{B}}_{q}{)}_{{i}_{1}{i}_{2}\cdots \phantom{\rule{thinmathspace}{0ex}}{i}_{qp}}={A}_{{j}_{1}{j}_{2}\cdots \phantom{\rule{thinmathspace}{0ex}}{j}_{p}}{B}_{{j}_{1}{j}_{2}\cdots \phantom{\rule{thinmathspace}{0ex}}{j}_{p}{i}_{1}{i}_{2}\cdots \phantom{\rule{thinmathspace}{0ex}}{i}_{qp}}$. The sum over $p$, finally, reflects the possibility of having not only one, but multiple types of patic order coexisting within the same system, as experiments on in vitro layers of MDCK cells have recently suggested (ArmengolCollado 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 righthand 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 pfold symmetric force complexion: i.e. ${\mathit{F}}_{c}=\sum _{k=1}^{p}{\mathit{F}}_{k}\delta (\mathit{r}{\mathit{r}}_{c}a{\mathit{\nu}}_{k})$ with $\mathit{F}}_{k$ the force exerted by a cell at each vertex and originating from the imbalance of the tensions $\mathit{T}}_{kl$, driven by the active contraction of the cellular junctions, converging at the kth vertex: i.e. $\mathit{F}}_{k}=\sum _{l}{\mathit{T}}_{kl$ (see Figure 1c). The quantities $\mathit{r}}_{c$ and $a$ are the cell’s centroid and circumradius, respectively, while $\mathit{\nu}}_{k}=\mathrm{cos}(\vartheta +2\pi k/p)\phantom{\rule{thinmathspace}{0ex}}{\mathit{e}}_{x}+\mathrm{sin}(\vartheta +2\pi k/p)\phantom{\rule{thinmathspace}{0ex}}{\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. ${\mathit{F}}_{k}\cdot {\mathit{\nu}}_{k}<0$) or extensile (${\mathit{F}}_{k}\cdot {\mathit{\nu}}_{k}>0$), depending on the overall tension distribution and the geometry of the cellular network. Next, assuming $\mathit{F}}_{k}=f{\mathit{\nu}}_{k$ and expanding the delta function about $a=0$ yields $\mathit{F}}_{c}=\sum _{m=0}^{\mathrm{\infty}}{\mathit{f}}_{m$, where
Because of the pfold symmetry of the force complexion $\mathit{f}}_{m}=\mathbf{0$ for all even $\displaystyle m$ values, unless $\displaystyle m=p1$, whereas odd $m$ values yields, up to symmetrization, $\sum _{k=1}^{p}{\mathit{\nu}}_{k}^{\otimes (m+1)}\sim {\mathbb{1}}^{\otimes (m+1)/2}$. Thus, after some algebraic manipulation, one finds $\mathit{F}}_{c}\approx apf/2\phantom{\rule{thinmathspace}{0ex}}\mathrm{\nabla}[(1+{a}^{2}/8\phantom{\rule{thinmathspace}{0ex}}{\mathrm{\nabla}}^{2}+\cdots )\delta (\mathit{r}{\mathit{r}}_{c})]+{\mathit{f}}_{p1$. Finally, taking $\u27e8\sum _{c}{\mathit{F}}_{c}\u27e9={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 $n=\u27e8\sum _{c}\delta (\mathit{r}{\mathit{r}}_{c})\u27e9$ is the cell number density. From Equation 5b, one finds the following expression for the phenomenological parameter $\displaystyle {\alpha}_{p}$ in Equation 3: i.e. ${\alpha}_{p}=(a{)}^{p1}pnf/[\sqrt{{2}^{p2}}\phantom{\rule{thinmathspace}{0ex}}(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 patic 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 $\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 secondorder active terms such as $Q}_{{k}_{1}{k}_{2}\dots \phantom{\rule{thinmathspace}{0ex}}{k}_{p}}{\mathrm{\partial}}_{i}{\mathrm{\partial}}_{j}{Q}_{{k}_{1}{k}_{2}\cdots \phantom{\rule{thinmathspace}{0ex}}{k}_{p}$, $Q}_{ij{k}_{3}\cdots \phantom{\rule{thinmathspace}{0ex}}{k}_{p}}{\mathrm{\partial}}_{{l}_{1}}{\mathrm{\partial}}_{{l}_{2}}{Q}_{{l}_{1}{l}_{2}{k}_{3}\cdots \phantom{\rule{thinmathspace}{0ex}}{k}_{p}$, etc., are mechanically equivalent to the terms $\mathrm{\partial}}_{i}{Q}_{{k}_{1}{k}_{2}\cdots \phantom{\rule{thinmathspace}{0ex}}{k}_{p}}{\mathrm{\partial}}_{j}{Q}_{{k}_{1}{k}_{2}\cdots \phantom{\rule{thinmathspace}{0ex}}{k}_{p}$ and $Q}_{{k}_{1}{k}_{2}\cdots \phantom{\rule{thinmathspace}{0ex}}i}{H}_{{k}_{1}{k}_{2}\cdots \phantom{\rule{thinmathspace}{0ex}}j}{H}_{{k}_{1}{k}_{2}\cdots \phantom{\rule{thinmathspace}{0ex}}i}{Q}_{{k}_{1}{k}_{2}\cdots \phantom{\rule{thinmathspace}{0ex}}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 $p$ value is considered. Such a crossover is expected at length scales larger than $\ell =({\alpha}_{p}/{\beta}_{p}{)}^{1/(p4)}$, where the second term of the righthand side of Equation 3 overweights the first term, reflecting the pfold symmetry of the local active forces. In the presence of multiple types of patic order, the pdependent 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 $F=\int \mathrm{d}A\phantom{\rule{thinmathspace}{0ex}}({f}_{2}+{f}_{6}+{f}_{2,6})$, where
Here, $A}_{p$ and $B}_{p$ are constants setting the magnitude of the order parameter at the length scale of the short distance cutoff, here assumed to be of the order of the cell size, and $\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 $\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 $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 higherorder 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 wellknown active nematic length scale, dictating both the hydrodynamic stability (Voituriez et al., 2005) and the largescale structure of spatiotemporal chaos in active nematics (Giomi, 2015) and whose signature in multicellular systems has been identified in both eukaryotes (BlanchMercader 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, $\ell}_{2$ and $\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 $\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 $\ell}_{\chi ,2}=\sqrt{{L}_{2}/{\chi}_{2,6}$ and $\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 ${\alpha}_{2}=0$ and ${\alpha}_{6}=0$, Equations 2 and 6, reduce to those of a twodimensional liquid crystal with coupled nematic and hexatic order parameter. The latter can be found, e.g., in freestanding liquid hexatic films (Dierker and Pindak, 1987; Sprunt and Litster, 1987), where molecules are either orthogonal to the midsurface 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 midsurface gives rise to inplane nematic order, which is coupled to the sixfold bondorientational 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 DrouinTouchette 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 $\ell \ll {\ell}_{6}$ and nematic order at length scales $\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. $\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 nematicdominated 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(\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 $k}_{\mathrm{d}}={k}_{\mathrm{a}$ and ${\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 cellresolved 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. $\rho \phantom{\rule{thinmathspace}{0ex}}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 $\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 $s}_{2}\sim {\alpha}_{2}^{2$. This effect is overestimated at the linear order, leading to an inverse quadratic dependence on the wave number $\mathit{q}$ (Ramaswamy et al., 2003), but is generally renormalized by nonlinearities, so that $\underset{\mathit{q}\to 0}{lim}S(\mathit{q})\sim \mathit{q}{}^{\alpha}$, with $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 $s}_{\beta}\sim {\alpha}_{6}^{2$ and the exponent $\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, $\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, $\beta$ differs depending on whether noise affects both cells’ orientational and translational dynamics, or only the former. Specifically, when only orientational noise is considered, $\beta =6$. In contrast, $\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.
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 cellresolved models of epithelia – i.e. the selfpropelled 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, selfpropelling at constant speed $v}_{0$ and whose direction of motion undergoes rotational diffusion with diffusion coefficient $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 powerlaw scaling regimes at small and large length scales. At small length scales, the structure factor scales like $S(\mathit{q})\sim \mathit{q}{}^{\beta}$, with $\beta$ monotonically decreasing from 6 to 4 upon increasing the Péclet number $\mathrm{P}\mathrm{e}={\xi}_{0}/a$ expressing the ratio between cells’ persistence length $\xi}_{0}={v}_{0}/{D}_{\mathrm{r}$ and their typical size $a$ (see Figure 3b).
Conversely, at large length scales, the structure factor scales like an inverse power law, with exponent consistent with the largescale behavior of active nematics (Chaté, 2020). These observations can be rationalized in the light of the previous fluctuating hydrodynamic analysis. In the limit $\mathrm{P}\mathrm{e}\to 0$, cells do not selfpropel, 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, $\beta =6$. Increasing $\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 $\beta =6$ to $\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 nonmonotonic behavior at small $\mathrm{P}\mathrm{e}$ values, as a crossover from a shearthinning to the shearthickening 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 crosscorrelation 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 $\mathit{r}$ compares to the length scales $\ell}_{\chi ,2$ and $\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 $\ell}_{\chi ,2}={\ell}_{\chi ,6}={\ell}_{\chi$, fluctuations dominate at short distances and the hexatic and nematic orientations are uncorrelated. Thus, ${C}_{26}(\mathit{r})$ is approximatively constant for $\mathit{r}\ll {\ell}_{\chi}$. The picture is reversed for $\mathit{r}\gg {\ell}_{\chi}$. In this range, the hexatic and nematic orientations are ‘locked’ in a parallel configuration, i.e. $\mathrm{Arg}({\psi}_{2})/2\approx \mathrm{Arg}({\psi}_{6})/6$, or tilted by $\pi /6$ with respect to each other, depending on the sign of the constant $\chi}_{2,6$, and the crosscorrelation function exhibits the standard powerlaw decay characterizing twodimensional liquid crystals with a singleorder parameter: i.e. $C}_{26}(\mathit{r})\sim (\mathit{r}/{\ell}_{\chi}{)}^{{\eta}_{26}$, with $\eta}_{26$ a specific instance of the generic nonuniversal exponent ${\eta}_{26}=6{k}_{B}T/(\pi K)$, with $K$ the orientational stiffness of both phases (proportional to $L}_{2}={L}_{6$). An analytical treatment of this simple case is reported in Methods. In the more generic case, in which $\ell}_{\chi ,2}\ne {\ell}_{\chi ,6$ and the relaxation rates of the hexatic and nematic phase differ, the crosscorrelation 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 ${\alpha}_{2}=0$ and ${\alpha}_{6}=0$, is shown in Figure 4a. The curves in Figure 4b correspond instead to simulated configurations of the crosscorrelation function of ${C}_{26}(\mathit{r})$ for finite hexatic and nematic activity. In this case, the crosscorrelation 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.
Taken together, our calculations of the structure factor and the crosscorrelation function demonstrate that the hydrodynamic theory embodied in Equations 2 and 6 is able to account for the multiscale hexanematic order observed in experiments (ArmengolCollado 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), largescale nematic order could arise from the selforganization 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).
Conclusions
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 patic 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érezGonzá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 cellresolved 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 (ArmengolCollado 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 shortscale 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. solidlike) and mesenchymal (i.e. liquidlike) 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.
Methods
Quantification of patic order in epithelial layers
Following ArmengolCollado et al., 2023, we use the shape function$\gamma}_{p$ to quantify the amount of pfold symmetry of an arbitrary cell. Denoting $\mathit{r}}_{v$ with $v=1,\phantom{\rule{thinmathspace}{0ex}}2\dots \phantom{\rule{thinmathspace}{0ex}}V$, the positions of its vertices with respect to the cell’s center of mass, one has
with ${\varphi}_{v}=\mathrm{Arg}({\mathit{r}}_{v})$ the angle between $\mathit{r}}_{v$ and the xaxis of a Cartesian frame. A schematic representation of these elements in an arbitrary irregular polygon is shown in Figure 5a. Unlike the complex function $\psi}_{p}={e}^{ip\vartheta$, which has unit magnitude by construction, the magnitude ${\gamma}_{p}$ quantify the resemblance of a generic polygon with a regular psided polygon of the same size, while the phase $\vartheta =\mathrm{Arg}({\gamma}_{p})/p$ marks the orientation of the polygon. For regular Vsided polygons, ${\gamma}_{p}=1$ provided $p$ is an integer multiple of $V$ and ${\gamma}_{p}\approx 0$ otherwise. Furthermore, from $\gamma}_{p$ one can readily compute.
Figure 5b, c shows examples of the functions $\gamma}_{2$ and $\gamma}_{6$ for a typical configuration of the SPV. We emphasize that $\gamma}_{p$, which, as shown in ArmengolCollado et al., 2023, arises from a pfold 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 $\gamma}_{p$ can then be coarse grained at the length scale $\ell$ to construct the shape parameter:
where the $\mathit{r}}_{c$ is the position of the cth cell, $\mathrm{\Theta}$ is the Heaviside step function, such that $\mathrm{\Theta}(x)=1$ for $x>0$ and 0 otherwise, and ${N}_{\ell}=\sum _{c}\mathrm{\Theta}(\ell \mathit{r}{\mathit{r}}_{c})$ is the number of cells within a distance $\ell$ from $\mathit{r}}_{c$. As in the case of $\gamma}_{p$, the magnitude of $\mathrm{\Gamma}}_{p$ reflects the resemblance between a multicelluar cluster and a regular psided 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 $\theta =\mathrm{Arg}({\mathrm{\Gamma}}_{2})/2$, while in panel (b), there are plotted streamlines showing the orientation of the director $\mathit{n}=\mathrm{cos}\theta \phantom{\rule{thinmathspace}{0ex}}{\mathit{e}}_{x}+\mathrm{sin}\theta \phantom{\rule{thinmathspace}{0ex}}{\mathit{e}}_{y}$.
Passive stresses
As explained in the main text, the passive contribution to the stress tensor is given by $\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 $\eta$ and $\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 ${\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 $P$ and $K}_{p$ denotes the orientational stiffness of the patic phase, related to the order parameter stiffness by
and $\mathit{varepsilon}$ is the twodimensional antisymmetric tensor, with ${\epsilon}_{xy}={\epsilon}_{yx}=1$ and ${\epsilon}_{xx}={\epsilon}_{yy}=0$.
Linear fluctuating hydrodynamics
To compute the structure factor, we follow Ramaswamy et al., 2003 and augment Equation 2b, Equation 2c with shortranged correlated noise field. Then calling $\vartheta$ and $\phi$ the nematic and hexatic fluctuating orientation fields and linearizing the hydrodynamic equations about the homogeneous and stationary solutions, $\vartheta =\phi =0$ and $\mathit{v}=\mathbf{0}$, gives
where $\delta \vartheta$, $\delta \phi$, and $\delta \mathit{v}$ indicate a small departure from the homogeneous and stationary configurations of the fields $\displaystyle \vartheta$, $\displaystyle \phi$, and $\mathit{v}$, $\mathcal{D}}_{p}={\mathrm{\Gamma}}_{p}{L}_{p$, $\chi}_{p}={\mathrm{\Gamma}}_{p}{\chi}_{2,6$, and $\xi}^{(\vartheta )$ and $\xi}^{(\phi )$ are shortranged correlated noise fields: i.e.
The velocity field $\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 $\mathit{f}}^{(\mathrm{p})}=\mathrm{\nabla}\cdot {\mathit{\sigma}}^{(\mathrm{p})$ and $\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 $\{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 ${\mathrm{\Xi}}^{(v)}\to 0$. The pressure $P$, in turn, can be related to the density by a linear equation of state of the form
with $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 $\hat{\mathit{M}}$ is given by
and the functions $\eta}^{(\alpha )$, with $\alpha \in \{\rho ,\vartheta ,\phi \}$, are effective noise fields whose correlation functions are given by
where the functions ${\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 $\xi}^{(\vartheta )$ and $\xi}^{(\phi )$, via the vorticity field on the righthand side of Equation 16b, Equation 16c, this effect disappears at the small (i.e. $\mathit{q}\to \mathrm{\infty}$) and large (i.e. $\mathit{q}\to 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(\mathit{q},\omega )$, can be calculated from the correlation function
To compute the lefthand side of Equation 29 one can solve Equation 24 with respect to $\delta \hat{\rho}$, $\delta \hat{\vartheta}$, and $\delta \hat{\phi}$. This gives
The static structure factor can then be expressed as
The first term on the righthand 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\to {\rho}_{0}\phantom{\rule{thinmathspace}{0ex}}{\mathrm{\Xi}}^{(\rho )}/(\varsigma {c}_{\mathrm{s}}^{2})$ when $\mathit{q}\to 0$. The effect of the active currents is instead accounted for by the second and third terms on the righthand side of Equation 31, which can be cast in the general form
where
The integral over $\omega$ can be derived using the residue theorem upon computing the roots of the complex thirdorder polynomial $h$. To make progress, we express
where $\omega}_{1$, $\omega}_{2$, and $\omega}_{3$ are given by
The integrand on the righthand side of Equation 33 has, therefore, three pairs of purely imaginary poles: i.e. $\pm i{\omega}_{1}$, $\pm i{\omega}_{2}$, and $\pm i{\omega}_{3}$. Next, turning the integration range to an infinite semicircular contour on the complex upper halfplane and summing the associated residues gives, after lengthy algebraic manipulations
where we have set
Now, although the individual elements of the matrix $\hat{\mathit{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. $\vartheta =\phi =0$ in this case). Because of the lack of longranged order and of specific directions that could affect the spectrum of density fluctuations, the latter is expected to be isotropic, thus $S=S(\mathit{q})$. To remove the fictitious angular dependence, one can either linearize Equation 2a about a generic pair of angles, $\vartheta}_{0$ and $\phi}_{0$, and then use these to calculate a circular average – i.e. $S(\mathit{q})=1/(2\pi {)}^{2}\int \mathrm{d}{\vartheta}_{0}\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}{\phi}_{0}\phantom{\rule{thinmathspace}{0ex}}S(\mathit{q})$ – or, more simply, by orienting $\mathit{q}$ so to cancel the directional dependence. Thus, taking $q}_{x}={q}_{y}=\mathit{q}/\sqrt{2$ gives a simpler expression of the matrix $\hat{\mathit{M}}$. That is
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 $s}_{2$ and $s}_{4$ in Equation 8. That is
Notice that, while both orientational and translation noise affect the amplitude of density fluctuations at small length scales, where $S(\mathit{q})\sim {s}_{4}\mathit{q}{}^{4}$, translational noise becomes unimportant at the large scale, where $S(\mathit{q})\sim {s}_{2}/\mathit{q}{}^{2}$. Furthermore, as long as viscous dissipation is at play, switching off translational noise (i.e. ${\mathrm{\Xi}}^{(\mathrm{v})}\to 0$) does not alter the scaling behavior of the structure factor at neither range of length scales. Taking the dry limit (i.e. $\eta \to 0$ and $\zeta \to 0$) leaves the largescale behavior unaltered, but does affect the scaling of density fluctuations at short length scales, where translational fluctuations are most prominent. Specifically, $S(\mathit{q})\sim {s}_{6}\mathit{q}{}^{6}$ in the case of purely rotational noise and $S(\mathit{q})\sim {s}_{10}\mathit{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 selfpropelled 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 $\mathit{r}}_{c$ of its center, with $c=1\phantom{\rule{thinmathspace}{0ex}},2\dots \phantom{\rule{thinmathspace}{0ex}}N$, and a velocity ${\mathit{v}}_{c}={v}_{0}(\mathrm{cos}{\theta}_{c}\phantom{\rule{thinmathspace}{0ex}}{\mathit{e}}_{x}+\mathrm{sin}{\theta}_{c}\phantom{\rule{thinmathspace}{0ex}}{\mathit{e}}_{y})$, with $v}_{0$ a constant speed and $\theta}_{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:
where µ is the mobility coefficient and $E=E({\mathit{r}}_{1},{\mathit{r}}_{2}\dots \phantom{\rule{thinmathspace}{0ex}}{\mathit{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 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 $P}_{c$) and the effective interfacial tension caused by the cell–cell adhesion and the cortical tension (both linear in $P}_{c$) (Farhadifar et al., 2007). The constants $A}_{0$ and $P}_{0$ represent, respectively, the preferred area and perimeter of each cell. The quantity $\eta}_{c$, on the other hand, is a random number with zero mean and correlation function
with $\mathcal{D}}_{\mathrm{r}$ a rotational diffusion coefficient. To make progress, we next introduce the following dimensionless numbers: the shape index $p}_{0}={P}_{0}/\sqrt{{A}_{0}$, which accounts for the spontaneous degree of acircularity of individual cells (Bi et al., 2016), and the Péclet number $\mathrm{P}\mathrm{e}={v}_{0}/({\mathcal{D}}_{\mathrm{r}}\sqrt{{A}_{0}})$, 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 $L}_{g$ with periodic boundary conditions. At $t=0$, the centroids $\mathit{r}}_{c$ are placed in a slightly perturbed hexagonal grid with a random initial velocity. After reaching the nonequilibrium steady state, we perform statistical averages of relevant observables.
In our numerical simulations, we set ${p}_{0}=3.85$, $\mu {K}_{A}{A}_{0}/{\mathcal{D}}_{\mathrm{r}}=1$, $\mu {K}_{P}/{\mathcal{D}}_{\mathrm{r}}=1$, and $\mathcal{D}}_{\mathrm{r}}\mathrm{\Delta}t=5\times {10}^{3$, where $\mathrm{\Delta}t$ is the timestep used for the integration, and the average density of particles $N{A}_{0}/{L}_{g}^{2}=1$. We vary the Péclet number in the range $0.1\le \mathrm{P}\mathrm{e}\le 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\le {L}_{g}\le 200$ at constant density. The density structure factor (light green circles) in Figure 3a was obtained, in particular, with $\mathrm{P}\mathrm{e}=1.5$.
The MPF model
The MPF model is a continuous model where each cell is described by a concentration field ${\phi}_{c}={\phi}_{c}(\mathit{r})$ with $\displaystyle c=1,\phantom{\rule{thinmathspace}{0ex}}2\dots \phantom{\rule{thinmathspace}{0ex}}N$ 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 $\mathcal{F}=\int \mathrm{d}A\phantom{\rule{thinmathspace}{0ex}}f$, where the free energy density $f$ is given by
Here, $\alpha$ and $k}_{\varphi$ are material parameters which can be used to tune the surface tension $\gamma =(8{k}_{\phi}\alpha /9{)}^{1/2}$ and the interfacial thickness $\xi =(2{k}_{\phi}/\alpha {)}^{1/2}$ of isolated cells and thermodynamically favor spherical cell shapes. The constant $\u03f5$ captures the repulsion between cells. The concentration field is large (i.e. $\phi}_{c}\simeq {\varphi}_{0$) inside the cells and zero outside. The contribution proportional to $\lambda$ in the free energy enforces cell incompressibility whose nominal radius is given by $R}_{\phi$. The relaxational dynamics of the field $\phi}_{c$ is governed by the Allen–Cahn equation
where $\mathit{v}}_{c$ 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 nonequlibrium cell migration. The dimensionless parameters of the model are the Péclet number $\mathrm{P}\mathrm{e}={v}_{0}/(2{\mathcal{D}}_{r}{R}_{\phi})$ and the cell deformability $d=\u03f5/\alpha$.
The system of partial differential equations, Equation 46, is solved with a finitedifference approach through a predictor–corrector finitedifference Euler scheme implementing secondorder stencil for space derivatives (Carenza et al., 2019). The Ccode 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 ${L}_{g}=380$ grid points. Model parameters in simulation units are as follows: ${R}_{\varphi}=11$, ${\phi}_{0}=2.0$, $M\alpha =0.006$, $M{k}_{\phi}=0.006$, $M\u03f5=0.01$, $M\lambda =600$, $M\gamma =0.008$, $D}_{r}\mathrm{\Delta}t={10}^{4$, being $\mathrm{\Delta}t$ the timestep used to integrate Equation 46. We vary the speed of selfpropulsion in the range $0.0\le {v}_{0}\le 0.005$. In terms of dimensionless parameters this corresponds to having $d=1.66$ and $\mathrm{P}\mathrm{e}$ 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\gamma )$ ranges between 0 and 0.625. Moreover, the nominal packing fraction is $N(\pi {R}_{\phi}^{2})/{L}_{g}^{2}=0.95$, while the ratio between the interface thickness and the nominal radius $\xi /{R}_{\phi}=0.12$. The density structure factor (dark green triangles) in Figure 3a was obtained with $\mathrm{P}\mathrm{e}=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 finitedifference Euler approach, with a firstorder upwind scheme and secondorder 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 ghostcell method to compute derivatives on the boundary of the computational subdomains. Runs have been performed using 64 CPUs in twodimensional geometries, on a computational box of size 256^{2} and 512^{2}, for at least $1.5\times {10}^{7}$ LB iterations (corresponding to ∼21 and ∼84 days of CPUtime, 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 $\rho =2.0$ everywhere. The model parameters in simulations units are as follows: $\eta =\zeta =1.66$, ${\lambda}_{2}={\lambda}_{6}=1.1$, ${\nu}_{2}={\nu}_{6}=0.0$, ${\mathrm{\Gamma}}_{2}=0.4$, ${A}_{2}={B}_{2}=0.04$, ${L}_{2}=0.04$, ${\mathrm{\Gamma}}_{6}=0.4$, ${A}_{6}={B}_{6}=0.004$, ${L}_{6}=0.004$, ${\kappa}_{2,6}={\xi}_{2,6}=0.004$. Nematic activity $\alpha}_{2$ has been varied in the range $0.02\le {\alpha}_{2}\le 0.0005$ and hexatic activity $\alpha}_{6$ in the range $0.050\le {\alpha}_{6}\le 0.050$. We set the active parameters $\beta}_{2$ and ${\beta}_{6}=0$. The density structure factor (continuous black line) in Figure 3a was obtained with $\alpha}_{2}=2\times {10}^{3$ and $\alpha}_{6}=2\times {10}^{2$.
The coherence length of the nematic and hexatic liquid crystal can be expressed as the $({L}_{p}/{A}_{p}{)}^{1/2}=\mathrm{\Delta}{x}_{\mathrm{L}\mathrm{B}}$ for both $p=2,6$, where $\mathrm{\Delta}{x}_{\mathrm{L}\mathrm{B}}$ 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 $\ell}_{2$ and ranges between $10\mathrm{\Delta}{x}_{\mathrm{L}\mathrm{B}}$ for ${\alpha}_{2}=0.0005$ and $1.5\mathrm{\Delta}{x}_{\mathrm{L}\mathrm{B}}$ for ${\alpha}_{2}=0.02$. Conversely, for hexatics $\ell}_{6$ and ranges up to $3.5\mathrm{\Delta}{x}_{\mathrm{L}\mathrm{B}}$ for ${\alpha}_{6}=0.05$. To compare the results of the hydrodynamics simulations with the discrete models in Figure 3a, we choose $2\mathrm{\Delta}{x}_{\mathrm{L}\mathrm{B}}=\sqrt{{A}_{0}}$ and $2\mathrm{\Delta}{x}_{\mathrm{L}\mathrm{B}}={R}_{\phi}\mathrm{\Delta}{x}_{\mathrm{M}\mathrm{P}}$, with $\mathrm{\Delta}{x}_{\mathrm{M}\mathrm{P}}$ 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 ${C}_{26}(\mathit{r})$ given in Equation 9, reflecting the amount of crosscorrelation in their fluctuations. Here, $\psi}_{2}={e}^{2i\vartheta$ and $\psi}_{6}={e}^{6i\phi$, while the fluctuating fields $\vartheta$ and $\phi$ represent again the local nematic and hexatic orientations, respectively. Averaging $\psi}_{2$ and $\psi}_{6$ over the scale of a volume element, yields the order complex parameters $\mathrm{\Psi}}_{2}=\u27e8{e}^{2i\vartheta}\u27e9={\mathrm{\Psi}}_{2}{e}^{2i\theta$ and $\mathrm{\Psi}}_{6}=\u27e8{e}^{6i\phi}\u27e9={\mathrm{\Psi}}_{6}{e}^{6i\varphi$, with $\theta$ and $\varphi$ the average orientations. To make progress, we assume that, at the scale of a volume element, both microscopic orientations $\vartheta$ and $\phi$ are Gaussianly distributed about their mean values, so that, in general
from which
This approximation holds when the relative fluctuation of the patic phase $\mathrm{Arg}({\psi}_{p})$ is sufficiently small, so that
consistent with the standard definition of patic order parameter. Thus, in particular, $\theta =\u27e8\vartheta \u27e9$ and ${\mathrm{\Psi}}_{2}=\u27e8\mathrm{cos}2(\vartheta \theta )\u27e9$, whereas $\varphi =\u27e8\phi \u27e9$ and ${\mathrm{\Psi}}_{6}=\u27e8\mathrm{cos}6(\phi \varphi )\u27e9$. This allows to write ${C}_{26}(\mathit{r})$, as given by Equation (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 Equation (51) is not strictly valid for a quasi longranged ordered liquid crystal, where also $\theta$ and $\varphi$ 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}_{\vartheta \phi}(\mathit{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 $\chi}_{2,6$ in Equation 6b. For ${\chi}_{2,6}<0$, the hexatic and nematic directors are energetically favored to be parallel, so that $\vartheta \approx \phi$. Conversely, when ${\chi}_{2,6}>0$, the hexatic and nematic directors are preferentially tilted by $\pi /6$, hence $\vartheta =\phi \pm \pi /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 $\chi}_{2,6$ values. Thus, assuming ${\chi}_{2,6}<0$ and expanding Equation 2c about $\vartheta \approx \phi$, gives
where, as in the previous sections, we have set $\mathcal{D}}_{p}={\mathrm{\Gamma}}_{p}{L}_{p$ and $\chi}_{p}={\mathrm{\Gamma}}_{p}{\chi}_{2,6$ and introduced the Gaussian noise fields $\xi}^{(\vartheta )$ and $\xi}^{(\vartheta )$, 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
where $\gamma}_{p}={K}_{p}/{\mathcal{D}}_{p$, with $K}_{p$ the orientational stiffness defined in Equation 15, is the rotational viscosity of the associated patic phase. Equation 53a can now be decoupled and used to compute the correlation function ${C}_{\vartheta \phi}(\mathit{r})$. For simplicity, here we set $\mathcal{D}}_{2}={\mathcal{D}}_{6}=\mathcal{D$, ${\gamma}_{2}={\gamma}_{6}=\gamma$, and $9{\chi}_{2}={\chi}_{6}=2\chi$. With this choice, taking
gives, after simple algebraic manipulations
where ${\xi}_{+}=({\xi}^{(\phi )}+{\xi}^{(\vartheta )})/2$ and ${\xi}_{}=({\xi}^{(\phi )}{\xi}^{(\vartheta )})/2$. Moreover, using Equation (54), one finds
where $\{n,m\}=\{+,\}$. Equation 56a can now be solved in Fourier space and real time to give
where the hat indicates Fourier transformation and
where ${m}_{+}=0$ and ${m}_{}^{2}={\ell}_{\chi}^{2}=\mathcal{D}/\chi $. The calculation of the crosscorrelation function ${C}_{\vartheta \phi}(\mathit{r})$ is now reduced to calculating the autocorrelation functions of the fields $\phi}_{+$ and $\phi}_{$. Specifically
where
and we have made use of Equation (54) to demonstrate that ${C}_{+}(\mathit{r})={C}_{+}(\mathit{r})=0$. The nonvanishing correlation functions, on the other hand, can be expressed as
where $\mathrm{\Lambda}=2\pi /a$ is a shortdistance cutoff and $\u27e8{\hat{\phi}}_{n}(\mathit{q},t){}^{2}\u27e9$ is the finitetime orientational structure factor defined from the relation
After standard algebraic manipulations one finds
from which Equation (62) can be calculated in the form
Evidently, Equation (65) is equivalent to that obtained in a purely static setting from the Hamiltonian
of the noninteracting scalar fields $\phi}_{+$ and $\phi}_{$. Now, in the case of the ‘massive’ field $\phi}_{$, the Fourier integral in Equation (65) converges to
in the range $\mathit{r}\gg a$. Here, $\displaystyle {K}_{0}$ is a modified Bessel function of the second kind, whose asymptotic expansion at short and long distances is given by
with $\gamma}_{\mathrm{E}\mathrm{M}$ the Euler–Mascheroni constant. In the case of the ‘massless‘ field $\phi}_{+$, on the other hand, the Fourier integral diverges in the infrared, but the correlation function ${C}_{++}(\mathit{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
Combining this with Equations (67) and (69) yields the following expression for the correlation function
where $\mathit{r}\gg a$. Finally, using Equation (50) and the asymptotic expansions of Equation (68) gives the following expression for the crosscorrelation function
where $\eta}_{26$ is an instance of the generic nonuniversal exponent
in the specific case $p=2$ and ${p}^{\mathrm{\prime}}=6$. Lastly, when $\displaystyle {\chi}_{2,6}>0$, the same procedure can be carried out by expanding Equation (2c) about $\vartheta =\phi \pm \pi /6$ and taking ${\phi}_{+}=(\phi +\vartheta )/2$ and ${\phi}_{}=(\phi \vartheta \pm \pi /6)/2$, from which one finds again Equation 72.
Data availability
Figures 2–4 contain the numerical data used to generate the figures.
References

Epithelia are multiscale active liquid crystalsNature Physics 19:1773–1779.https://doi.org/10.1038/s41567023021790

Geometric constraints during epithelial jammingNature Physics 14:613–620.https://doi.org/10.1038/s4156701800899

A texture tensor to quantify deformationsGranular Matter 5:67–70.https://doi.org/10.1007/s100350030126x

A densityindependent rigidity transition in biological tissuesNature Physics 11:1074–1079.https://doi.org/10.1038/nphys3471

Motilitydriven glass and jamming transitions in biological tissuesPhysical Review. X 6:021011.https://doi.org/10.1103/PhysRevX.6.021011

Turbulent dynamics of epithelial cell culturesPhysical Review Letters 120:208101.https://doi.org/10.1103/PhysRevLett.120.208101

Jamming of deformable polygonsPhysical Review Letters 121:248003.https://doi.org/10.1103/PhysRevLett.121.248003

Forces driving epithelial wound healingNature Physics 10:683–690.https://doi.org/10.1038/nphys3040

Hexatic order and herringbone packing in liquid crystalsPhysical Review Letters 48:1625–1628.https://doi.org/10.1103/PhysRevLett.48.1625

Lattice boltzmann methods and active fluidsThe European Physical Journal. E, Soft Matter 42:81.https://doi.org/10.1140/epje/i2019118436

Dry aligning dilute active matterAnnual Review of Condensed Matter Physics 11:189–212.https://doi.org/10.1146/annurevconmatphys031119050752

Dynamics of thin tilted hexatic liquid crystal filmsPhysical Review Letters 59:1002–1005.https://doi.org/10.1103/PhysRevLett.59.1002

Emergent potts order in a coupled hexaticnematic XY modelPhysical Review X 12:011043.https://doi.org/10.1103/PhysRevX.12.011043

Thermally driven orderdisorder transition in twodimensional soft cellular systemsPhysical Review Letters 123:188001.https://doi.org/10.1103/PhysRevLett.123.188001

Collective cell migration in morphogenesis, regeneration and cancerNature Reviews. Molecular Cell Biology 10:445–457.https://doi.org/10.1038/nrm2720

Geometry and topology of turbulence in active nematicsPhysical Review X 5:031003.https://doi.org/10.1103/PhysRevX.5.031003

Hydrodynamic theory of patic liquid crystalsPhysical Review. E 106:024701.https://doi.org/10.1103/PhysRevE.106.024701

Longranged order and flow alignment in sheared patic liquid crystalsPhysical Review Letters 129:067801.https://doi.org/10.1103/PhysRevLett.129.067801

Discrete rearranging disordered patterns, part I: robust statistical tools in two or three dimensionsThe European Physical Journal. E, Soft Matter 25:349–369.https://doi.org/10.1140/epje/i2007102988

Instabilities and geometry of growing tissuesPhysical Review Letters 129:048102.https://doi.org/10.1103/PhysRevLett.129.048102

Collective cancer invasion forms an integrindependent radioresistant nicheThe Journal of Experimental Medicine 217:e20181184.https://doi.org/10.1084/jem.20181184

Description of cellular patterns by dirichlet domains: the twodimensional caseJournal of Theoretical Biology 72:523–543.https://doi.org/10.1016/00225193(78)903156

From cells to tissue: a continuum model of epithelial mechanicsPhysical Review. E 96:022418.https://doi.org/10.1103/PhysRevE.96.022418

Embryo as an active granular fluid: stresscoordinated cellular constriction chainsJournal of Physics. Condensed Matter 28:414021.https://doi.org/10.1088/09538984/28/41/414021

Role of cell deformability in the twodimensional melting of biological tissuesPhysical Review Materials 2:045602.https://doi.org/10.1103/PhysRevMaterials.2.045602

Solidliquid transition of deformable and overlapping active particlesPhysical Review Letters 125:038003.https://doi.org/10.1103/PhysRevLett.125.038003

Discrete rearranging disordered patterns, part II: 2D plasticity, elasticity and flow of a foamThe European Physical Journal. E, Soft Matter 25:371–384.https://doi.org/10.1140/epje/i2007103007

Emergence of active nematic behavior in monolayers of isotropic cellsPhysical Review Letters 122:048004.https://doi.org/10.1103/PhysRevLett.122.048004

A dynamic cell model for the formation of epithelial tissuesPhilosophical Magazine B 81:699–719.https://doi.org/10.1080/13642810108205772

Unjamming and cell shape in the asthmatic airway epitheliumNature Materials 14:1040–1048.https://doi.org/10.1038/nmat4357

Hexatic phase in a model of active biological tissuesSoft Matter 16:3914–3920.https://doi.org/10.1039/d0sm00109k

Active wetting of epithelial tissuesNature Physics 15:79–88.https://doi.org/10.1038/s4156701802795

Active dynamics of tissue shear flowNew Journal of Physics 19:033006.https://doi.org/10.1088/13672630/aa5756

Active nematics on a substrate: giant number fluctuations and longtime tailsEurophysics Letters 62:196–202.https://doi.org/10.1209/epl/i2003003467

Cell intercalation in a simple epitheliumPhilosophical Transactions of the Royal Society of London. Series B, Biological Sciences 375:20190552.https://doi.org/10.1098/rstb.2019.0552

Theory of transitions among tilted hexatic phases in liquid crystalsPhysical Review. A, General Physics 39:3135–3147.https://doi.org/10.1103/physreva.39.3135

Dynamics of tilted hexatic phases in liquidcrystal filmsJournal de Physique II 1:1363–1373.https://doi.org/10.1051/jp2:1991145

Mechanical waves during tissue expansionNature Physics 8:628–634.https://doi.org/10.1038/nphys2355

Lownoise phase of a twodimensional active nematic systemPhysical Review. E 97:012707.https://doi.org/10.1103/PhysRevE.97.012707

Lightscattering study of bond orientational order in a tilted hexatic liquidcrystal filmPhysical Review Letters 59:2682–2685.https://doi.org/10.1103/PhysRevLett.59.2682

Linear viscoelastic properties of the vertex model for epithelial tissuesPLOS Computational Biology 18:e1010135.https://doi.org/10.1371/journal.pcbi.1010135

Physical forces during collective cell migrationNature Physics 5:426–430.https://doi.org/10.1038/nphys1269

Spontaneous flow transition in active polar gelsEurophysics Letters 70:404–410.https://doi.org/10.1209/epl/i2004105012

Epithelialtomesenchymal transition in cancer: complexity and opportunitiesFrontiers of Medicine 12:361–373.https://doi.org/10.1007/s1168401806566

Dynamics of twodimensional meltingPhysical Review B 22:2514–2541.https://doi.org/10.1103/PhysRevB.22.2514
Article and author information
Author details
Funding
European Research Council
 Livio Nicola Carenza
 Luca Giomi
Nederlandse Organisatie voor Wetenschappelijk Onderzoek
 JosepMaria ArmengolCollado
The funders had no role in study design, data collection, and interpretation, or the decision to submit the work for publication.
Acknowledgements
We are indebted with Massimo Pica Ciamarra and David Nelson for insightful discussions. This work is supported by the ERCCoG 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 ScienceXL 2019.022 (JM AC). Part of this work was carried out on the Dutch national einfrastructure with the support of SURF through the NWO Grant 2021.028 for computational time.
Copyright
© 2024, ArmengolCollado, 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.
Metrics

 825
 views

 199
 downloads

 1
 citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
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)
Further reading

 Microbiology and Infectious Disease
 Physics of Living Systems
Filamentous multicellular cable bacteria perform centimeterscale electron transport in a process that couples oxidation of an electron donor (sulfide) in deeper sediment to the reduction of an electron acceptor (oxygen or nitrate) near the surface. While this electric metabolism is prevalent in both marine and freshwater sediments, detailed electronic measurements of the conductivity previously focused on the marine cable bacteria (Candidatus Electrothrix), rather than freshwater cable bacteria, which form a separate genus (Candidatus Electronema) and contribute essential geochemical roles in freshwater sediments. Here, we characterize the electron transport characteristics of Ca. Electronema cable bacteria from Southern California freshwater sediments. Current–voltage measurements of intact cable filaments bridging interdigitated electrodes confirmed their persistent conductivity under a controlled atmosphere and the variable sensitivity of this conduction to air exposure. Electrostatic and conductive atomic force microscopies mapped out the characteristics of the cell envelope’s nanofiber network, implicating it as the conductive pathway in a manner consistent with previous findings in marine cable bacteria. Fourprobe measurements of microelectrodes addressing intact cables demonstrated nanoampere currents up to 200 μm lengths at modest driving voltages, allowing us to quantify the nanofiber conductivity at 0.1 S/cm for freshwater cable bacteria filaments under our measurement conditions. Such a high conductivity can support the remarkable sulfidetooxygen electrical currents mediated by cable bacteria in sediments. These measurements expand the knowledgebase of longdistance electron transport to the freshwater niche while shedding light on the underlying conductive network of cable bacteria.

 Computational and Systems Biology
 Physics of Living Systems
Bcell repertoires are characterized by a diverse set of receptors of distinct specificities generated through two processes of somatic diversification: V(D)J recombination and somatic hypermutations. B cell clonal families stem from the same V(D)J recombination event, but differ in their hypermutations. Clonal families identification is key to understanding Bcell repertoire function, evolution and dynamics. We present HILARy (Highprecision Inference of Lineages in Antibody Repertoires), an efficient, fast and precise method to identify clonal families from single or pairedchain repertoire sequencing datasets. HILARy combines probabilistic models that capture the receptor generation and selection statistics with adapted clustering methods to achieve consistently high inference accuracy. It automatically leverages the phylogenetic signal of shared mutations in difficult repertoire subsets. Exploiting the high sensitivity of the method, we find the statistics of evolutionary properties such as the site frequency spectrum and 𝑑𝑁∕𝑑𝑆 ratio do not depend on the junction length. We also identify a broad range of selection pressures spanning two orders of magnitude.