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 p-atic liquid crystals, we demonstrate that collectively migrating epithelia feature both nematic (i.e. p = 2) and hexatic (i.e. p = 6) orders, with the former being dominant at large and the latter at small length scales. Such a remarkable multiscale liquid crystal order leaves a distinct signature in the system’s structure factor, which exhibits two different power-law scaling regimes, reflecting both the hexagonal geometry of small cells clusters and the uniaxial structure of the global cellular flow. We support these analytical predictions with two different cell-resolved models of epithelia – i.e. the self-propelled Voronoi model and the multiphase field model – and highlight how momentum dissipation and noise influence the range of fluctuations at small length scales, thereby affecting the degree of cooperativity between cells. Our construction provides a theoretical framework to conceptualize the recent observation of multiscale order in layers of Madin–Darby canine kidney cells and pave the way for further theoretical developments.
Editor's evaluation
This important work presents a hydrodynamic description of confluent epithelial monolayers that captures different forms of orientational order in a scale dependent fashion and couples this order with flows driven by active stresses. Solid evidence for the validity of this approach is provided by detailed numerical simulations of different model tissues. This work should be of interest to a broad range of biophysicists interested in tissue mechanics and active matter.
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 (Serra-Picamal et al., 2012) or invading the extracellular matrix (Haeger et al., 2020), multicellular systems can move coherently within their micro-environment and coordinate the dynamics of their actin cytoskeleton, while retaining cell–cell contacts. This ability lies at the heart of a myriad of processes that are instrumental for life, such as embryonic morphogenesis and wound healing, but also of life-threatening conditions, such as metastatic cancer.
Understanding the physical origin of this behavior inevitably demands reliable theoretical models, aimed at providing a conceptual framework for dissecting and deciphering the wealth of biophysical data stemming from in vitro experiments and in vivo observations. Following the pioneering works by Honda, 1978; Nagai and Honda, 2001; Farhadifar et al., 2007; Bi et al., 2015; Bi et al., 2016, and others (Boromand et al., 2018; Mueller et al., 2019; Loewe et al., 2020; Monfared et al., 2021), cell-resolved models have played so far the leading role in this endeavour. Taking inspiration from the physics of foams (Graner et al., 2008; Marmottant et al., 2008), these models portray a confluent tissue as a collection of adjacent or overlapping polygonal cells (Figure 1a, b), whose dynamics is assumed to be governed by a set of overdamped Langevin equations, expressing the interplay between cells’ autonomous motion and remodeling events, which change the local topology of the cellular networks.
Despite their conceptual simplicity, cell-resolved models agree remarkably well with experimental data on confluent monolayers (Park et al., 2015; Atia et al., 2018). In particular they account for a solid-to-liquid transition controlled by the cells velocity and their compliance to deformations (Bi et al., 2015; Bi et al., 2016; Loewe et al., 2020). Furthermore, as demonstrated by Pica Ciamarra and coworkers, the solid and isotropic liquid states of these model-epithelia are separated by an intermediate hexatic phase, in which the system exhibits the typical sixfold rotational symmetry of two-dimensional crystals and yet is able to flow (Li and Ciamarra, 2018; Pasupalak et al., 2020). Shortly after discovery, the same property has been recovered within the framework of the cellular Potts model, thereby strengthening the idea that hexatic order may in fact serve as a guiding principle to unravel the collective dynamics of confluent epithelia (Durand and Heu, 2019). Furthermore, recent in vitro studies of Madin–Darby canine kidney (MDCK) cell layers demonstrated that epithelial layers can in fact feature both nematic and hexatic orders, with the former being dominant at large and the latter at short length scales (see Figure 1a, b and Armengol-Collado et al., 2023; Eckert et al., 2023). This remarkable example of physical organization in biological matter, referred to as multiscale hexanematic order in Armengol-Collado et al., 2023, is believed to complement the complex network or regulatory pathways available to individual cells to achieve multicellular organization and select specific scale-dependent collective migration strategies.
Motivated by these recent discoveries, in this article we propose a continuum theory of confluent epithelia rooted in the hydrodynamics of liquid crystals with generic p-atic rotational symmetry (hereafter p-atic liquid crystals). Previous theories of epithelial hydrodynamics can be schematically grouped in two categories: (1) models based on (isotropic/polar/nematic) active gels (Ranft et al., 2010; Popović et al., 2017; Pérez-González et al., 2019); (2) models built around the so-called shape tensor (Ishihara et al., 2017; Czajkowski et al., 2018; Hernandez and Marchetti, 2021; Grossman and Joanny, 2022), i.e. a rank-2 tensor, similar to the inertia tensor in kinematics, that embodies the geometrical structure of the polygonal cells. Although both classes of models hold great heuristic value and represent a solid foundation for any future development, they suffer from the same limitation: being based on a tensorial order parameter whose rank is two or less, they can account at most for twofold rotational symmetry (i.e. nematic order), while leaving the small-scale hexatic order unresolved. To overcome this limitation, here we exploit recent advances toward extending the classic hydrodynamic theory of hexatic liquid crystals (Zippelius, 1980; Zippelius et al., 1980) to account for arbitrary p-fold rotational symmetry order (Giomi et al., 2022a; Giomi et al., 2022b), with and being the most relevant cases (but possibly not the only) in the context of epithelial dynamics. We demonstrate that multiscale order is inherent to active liquid crystals with coupled order parameters, because of the indissoluble connection between shape and forces characterizing this class of non-equilibrium systems. Using fluctuating hydrodynamics, we explicitly compute the structure factor of epithelial layers and unveil a fascinating interplay between the nature of momentum dissipation (i.e. viscosity or friction) and noise at short length scales, where hexatic order is dominant. Such a mechanism profoundly affects the range of density fluctuations and could be harnessed to control the degree of collectiveness of cellular motion. Finally, by testing predictions against two different microscopic models of epithelia we demonstrate the robustness of multiscale hexanematic order across the rich landscape of models of epithelia.
Results and discussion
The model
Two-dimensional p-atic liquid crystals are traditionally described in terms of the orientation field , with the local orientation of the p-fold mesogens. A more general approach, proposed in Giomi et al., 2022a; Giomi et al., 2022b and especially suited for hydrodynamics, revolves instead around the rank-p tensor order parameter, with and , constructed upon averaging the pth tensorial power of the local orientation . That is
where denotes the ensemble average and the operator has the effect of rendering an arbitrary tensor traceless and symmetric (Hess, 2015). The vector 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 . The fields and represent, respectively, the magnitude and phase of the complex p-atic order parameter , while the normalization factor is chosen so that for all values. For , Equation (1) readily gives the standard nematic order parameter tensor: i.e. , with 1 the identity tensor. In practice, if a cell’s planar projection consists of a regular p-sided polygon, the microscopic orientation equates that of any of the vertices of the polygon. In the more realistic case of an irregular polygon, on the other hand, is given by the phase of the complex function , arising form the p-fold generalization of the classic shape tensor (Aubouy et al., 2003). This function was introduced in Armengol-Collado et al., 2023 and is reviewed in Methods for sake of completeness.
The order parameter tensor , the mass density , and the momentum density , with the local velocity field, comprise the set of hydrodynamic variables describing the dynamics of a generic p-atic fluid, which in turn is governed by the following set of partial differential equations (Giomi et al., 2022a; Giomi et al., 2022b):
where . Equation (2a) and Equation 2b are the mass and momentum conservation equations, with and rates of cell division and apoptosis, the stress tensor and an arbitrary external force per unit area. In Equation (2c), is a rotational viscosity and is the molecular tensor describing the relaxation of the p-atic phase toward the minimum of the free energy (see Methods). The rank-2 tensors and , with indicating transposition, are the vorticity and strain rate tensors, respectively, whereas the dot product in the first line of the equation implies a contraction of one index of with one of : i.e. . On the second line , while denotes the floor function and is zero for even values and one for odd values. Finally, , , and are material parameters expressing the strength of the coupling between p-atic order and flow.
Now, in order for Equation 2a to account for the dynamics of epithelial cell layers, we must specify the structure of the external force in Equation 2b and the stress tensor . As cells collectively crawl on a substrate, at a speed of order 0.1–1 µm/min (Brugués et al., 2014; Angelini et al., 2011), the former can be model as a Stokesian drag: , with a drag coefficient. A more realistic treatment of the interplay between the cells and the substrate would account for the traction forces exerted by the cells’ cryptic lamellipodium as well as for the compliance of the substrate (Trepat et al., 2009) and will be considered in the future. The stress tensor, on the other hand, is routinely decomposed into a passive and an active component: i.e. . The passive stress tensor is in turn expressed as , where is the pressure, is the elastic stress, arising in response of a static deformation of a fluid patch, and and are, respectively, the reactive (i.e. energy preserving) and viscous (i.e. energy dissipating) stresses originating from the reversible and irreversible couplings between p-atic order and flow. The generic expression of was derived in Giomi et al., 2022b and is reported in Methods.
The active stress , on the other hand, can be constructed phenomenologically for arbitrary values in the form
where the symbol denotes a contraction of all matching indices of the two operands and yields a tensor whose rank equates the number of unmatched indices: i.e. letting and be two generic tensors of rank , then . The sum over , finally, reflects the possibility of having not only one, but multiple types of p-atic order coexisting within the same system, as experiments on in vitro layers of MDCK cells have recently suggested (Armengol-Collado et al., 2023; Eckert et al., 2023).
Before exploring the consequences of the latter assumption, some comment about the physical interpretation of the terms featured in Equation 3 is in order. The first term on the right-hand side of Equation 3 is the stress resulting from the contractile or extensile forces exerted at the length scale of individual cells. To illustrate this concept one can assume each cell to exert a p-fold symmetric force complexion: i.e. with the force exerted by a cell at each vertex and originating from the imbalance of the tensions , driven by the active contraction of the cellular junctions, converging at the kth vertex: i.e. (see Figure 1c). The quantities and are the cell’s centroid and circumradius, respectively, while . We stress that, while the individual tensions acting along the junctions are exclusively contractile, the resulting vertex forces can be either contractile (i.e. ) or extensile (), depending on the overall tension distribution and the geometry of the cellular network. Next, assuming and expanding the delta function about yields , where
Because of the p-fold symmetry of the force complexion for all even values, unless , whereas odd values yields, up to symmetrization, . Thus, after some algebraic manipulation, one finds . Finally, taking 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 is the cell number density. From Equation 5b, one finds the following expression for the phenomenological parameter in Equation 3: i.e. . Notice that both constants and involved in Equation 5a are, in general, order dependent. We will come back on this aspect in Conclusion.
The second term in Equation 3, in contrast, expresses the active stress resulting from the spatial variations of the p-atic order parameter and, although similar to other contributions to the passive stress , cannot be derived from equilibrium considerations. Other terms constructed by contracting with can be expressed as linear combinations of this and , thus lead to a mere renormalization of the material parameters. It must be noted that the stress tensor enters in Equation 2b only via its divergence. Thus, possible second-order active terms such as , , etc., are mechanically equivalent to the terms and 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 value is considered. Such a crossover is expected at length scales larger than , where the second term of the right-hand side of Equation 3 overweights the first term, reflecting the p-fold symmetry of the local active forces. In the presence of multiple types of p-atic order, the p-dependent structure of the active stress renders the multiscale nature of the system enormously more dramatic. To illustrate this crucial point, here we postulate the system to behave as a hexanematic liquid crystal. Formally, such a scenario can be accounted by simultaneously solving two variants of Equation 2c, for and . 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 , where
Here, and are constants setting the magnitude of the order parameter at the length scale of the short distance cut-off, here assumed to be of the order of the cell size, and 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 , 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 can further be augmented with several additional terms of higher differential order: e.g. , , , etc. For simplicity, here we ignore these and higher-order couplings and focus on the zeroth order terms included in Equation 6b.
Crucially, Equation 3, Equation 6a entail two length scales, reflecting the distance at which the passive torques originating from the entropic elasticity of the nematic and hexatic phases counterbalance those arising from the active stresses:
The former is the well-known active nematic length scale, dictating both the hydrodynamic stability (Voituriez et al., 2005) and the large-scale structure of spatiotemporal chaos in active nematics (Giomi, 2015) and whose signature in multicellular systems has been identified in both eukaryotes (Blanch-Mercader et al., 2018) and prokaryotes (You et al., 2018). The latter, on the other hand, sets the typical size of hexatic domains at the small length scale. Remarkably, and 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 . 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 and . Their role will be discussed in the following section, in the framework of fluctuating hydrodynamics.
Finally, in the passive limit, when and , Equations 2 and 6, reduce to those of a two-dimensional liquid crystal with coupled nematic and hexatic order parameter. The latter can be found, e.g., in free-standing liquid hexatic films (Dierker and Pindak, 1987; Sprunt and Litster, 1987), where molecules are either orthogonal to the mid-surface of the film or tilted by a fixed angle. In the latter case, the projection of the average molecular direction on the tangent plane of the mid-surface gives rise to in-plane nematic order, which is coupled to the sixfold bond-orientational order associated with the underlying hexatic phase (see e.g. Bruinsma and Aeppli, 1982; Selinger and Nelson, 1989; Selinger, 1991 for a theoretical account and Drouin-Touchette et al., 2022 for recent developments). As we will detail in the following, activity profoundly alters this scenario by acting as a mechanical bandpass filter, which renders hexatic order dominant at length scales and nematic order at length scales . 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. , there is no dominant order and the system’s collective behavior is determined by the complex interplay of competing active and passive effects. To make progress, here we focus on the most dramatic hexatic- and nematic-dominated behaviors and treat intermediate length scales as simply as possible.
Multiscale order in epithelia
To elucidate the multiscale organization of the system, we next compute the structure factor , 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 and for simplicity. We stress that the validity of this approximation is strictly related with the present comparison between the hydrodynamic theory presented in this article and cell-resolved models. An assessment of the relevance of this and the other material parameters featured in Equation 2a can only be achieved via experimental scrutiny and is likely to depend on the specific cell type and environmental conditions. Furthermore, as the typical Reynolds number of collective epithelial flow is in the range 10−7–10−6, we neglect inertial effects: i.e. . 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 , 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 . This effect is overestimated at the linear order, leading to an inverse quadratic dependence on the wave number (Ramaswamy et al., 2003), but is generally renormalized by nonlinearities, so that , with (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 and the exponent determined by the specific energy dissipation mechanism, as well as by the specific structure of the noise. As detailed in Methods, here we consider four alternative scenarios, obtained upon combining two different momentum dissipation mechanisms (i.e. viscosity and friction) with two different types of noise (i.e. rototranslational and purely rotational). In the presence of viscous dissipation, i.e. a regime referred to as ‘wet‘ in the jargon of active matter, irrespective of the nature of noise. Conversely, in the ‘dry‘ limit, when the shear and bulk viscosity vanish and momentum dissipation solely results from the frictional interactions with the substrate, differs depending on whether noise affects both cells’ orientational and translational dynamics, or only the former. Specifically, when only orientational noise is considered, . In contrast, 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 cell-resolved models of epithelia – i.e. the self-propelled Voronoi (SPV) model (Bi et al., 2016) and the multiphase field (MPF) model (Loewe et al., 2020) (see the insets Figure 3b for typical configurations of the two models) – with that resulting from a numerical integration of Equation 2a (Carenza et al., 2019; Carenza et al., 2020), with none of the simplifications behind Equation (8). In both microscopic models, cells are treated as persistent random walkers, self-propelling at constant speed and whose direction of motion undergoes rotational diffusion with diffusion coefficient (see Methods for details). Noise is therefore expected to affect both the rotational and translational dynamics of the cell monolayer, although in a way that, unlike in our analytical treatment, cannot be trivially decoupled. Consistently with our linear analysis, both data sets exhibit two different power-law scaling regimes at small and large length scales. At small length scales, the structure factor scales like , with monotonically decreasing from 6 to 4 upon increasing the Péclet number expressing the ratio between cells’ persistence length and their typical size (see Figure 3b).
Conversely, at large length scales, the structure factor scales like an inverse power law, with exponent consistent with the large-scale behavior of active nematics (Chaté, 2020). These observations can be rationalized in the light of the previous fluctuating hydrodynamic analysis. In the limit , cells do not self-propel, noise is predominantly orientational and momentum propagates only at distances comparable to the average cell size. Under this circumstances, an in silico cell layer, whether modeled via the SPV or the MPF, behaves therefore as a ‘dry’ active system subject to purely rotational noise, for which, consistently with our analysis, . Increasing 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 to . The simple linear calculation, summarized in Methods, does not allow us to resolve the full crossover, but does provide a precise estimate of the upper and lower bounds. Finally, along the wet–dry crossover, viscosity must emerge from the cells’ lateral interactions. A precise understanding of this process is outside of the scope of the present work, but recent numerical work on the Vertex model has already highlighted the existence of a rich landscape of exotic rheological phenomena, resulting from the interplay between cellular motion, morphology, and adhesion (Tong et al., 2022; Hertaeg et al., 2022). The latter could possibly explain the non-monotonic behavior at small values, as a crossover from a shear-thinning to the shear-thickening behavior (Hertaeg et al., 2022) for additional numerical evidence of this effect.
A different signature of multiscale hexanematic order can be identified in the structure of the cross-correlation function
At equilibrium, and if deformations are sufficiently gentle to render backflow effects negligible, its behavior can be divided in two regimes, depending on how the distance compares to the length scales and 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 , fluctuations dominate at short distances and the hexatic and nematic orientations are uncorrelated. Thus, is approximatively constant for . The picture is reversed for . In this range, the hexatic and nematic orientations are ‘locked’ in a parallel configuration, i.e. , or tilted by with respect to each other, depending on the sign of the constant , and the cross-correlation function exhibits the standard power-law decay characterizing two-dimensional liquid crystals with a single-order parameter: i.e. , with a specific instance of the generic non-universal exponent , with the orientational stiffness of both phases (proportional to ). An analytical treatment of this simple case is reported in Methods. In the more generic case, in which and the relaxation rates of the hexatic and nematic phase differ, the cross-correlation function has a less standard functional form, but still features a slow and fast decay regime at short and large distances, respectively. An example of such a scenario, obtained from a numerical integration of Equation 2a with and , is shown in Figure 4a. The curves in Figure 4b correspond instead to simulated configurations of the cross-correlation function of for finite hexatic and nematic activity. In this case, the cross-correlation function exhibits an oscillatory behavior at short distances and vanishes at a length scale that becomes progressively large as the hexatic activity is increased. Consistently with our previous analysis, this latter feature confirms the existence of a hierarchy of orientationally ordered structures nested into each other at different length scales.
Taken together, our calculations of the structure factor and the cross-correlation function demonstrate that the hydrodynamic theory embodied in Equations 2 and 6 is able to account for the multiscale hexanematic order observed in experiments (Armengol-Collado et al., 2023; Eckert et al., 2023) and harnesses it into a continuum mechanical framework. Whereas the origin of hexanematic order is still a matter of investigation, the current experimental and numerical evidence suggests that, similarly to granular materials (Majmudar and Behringer, 2005), large-scale nematic order could arise from the self-organization of the microscopic force hexapoles into force chains. The possibility of similarity between these two phenomena has also been in relation to the initial phase of Drosophila gastrulation, where linear arrays of cells simultaneously undergo apical constriction in the ventral furrow region (Jason Gao et al., 2016).
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 p-atic liquid crystals (Giomi et al., 2022a; Giomi et al., 2022b). This approach allows one to account for arbitrary discrete rotational symmetries, thereby going beyond existing hydrodynamic theories of epithelia (Ranft et al., 2010; Popović et al., 2017; Pérez-González et al., 2019; Ishihara et al., 2017; Czajkowski et al., 2018; Hernandez and Marchetti, 2021; Grossman and Joanny, 2022), where the algebraic structure of the hydrodynamic variables renders impossible to account for liquid crystal order other than isotropic (i.e. ), polar (i.e. ), or nematic (i.e. ). Upon computing the static structure factor and comparing this with the outcome of two different cell-resolved models – i.e. the SPV (Bi et al., 2016) and MPF (Loewe et al., 2020) models – we have shown that, consistently with recent experimental findings (Armengol-Collado et al., 2023; Eckert et al., 2023), epithelial layers may in fact comprise both nematic and hexatic (i.e. ) order, coexisting at different length scales. Although the consequences of such a remarkable versatility are yet to be explored, we expect hexatic order to be relevant for short-scale remodeling events, where the local nature of hexatic order, combined with the rich dynamics of hexatic defects (Zippelius et al., 1980; Amir and Nelson, 2012), may mediate processes such as cell intercalations and the rearrangement of multicellular rosettes (Blankenship et al., 2006; Rauzi, 2020). Such a local motion, in turn, may be coordianted at the large scale by the underlying nematic order, giving rise to a persistent unidirectional flow, such as that observed during wound healing and cancer progression (Friedl and Gilmour, 2009). Furthermore, the existence of multiscale liquid crystal order echoes the most recent understanding of phenotypic plasticity in tissues, according to which the epithelial (i.e. solid-like) and mesenchymal (i.e. liquid-like) states represent the two ends of a spectrum of intermediate phenotypes (Zhang and Weinberg, 2018). These intermediate states display distinctive cellular characteristics, including adhesion, motility, stemness and, in the case of cancer cells, invasiveness, drug resistance, etc. Can multiscale liquid crystal order help understanding how the biophysical properties of tissues vary along the epithelial–mesenchymal spectrum? This and related questions will be addressed in the near future.
Methods
Quantification of p-atic order in epithelial layers
Following Armengol-Collado et al., 2023, we use the shape function to quantify the amount of p-fold symmetry of an arbitrary cell. Denoting with , the positions of its vertices with respect to the cell’s center of mass, one has
with the angle between and the x-axis of a Cartesian frame. A schematic representation of these elements in an arbitrary irregular polygon is shown in Figure 5a. Unlike the complex function , which has unit magnitude by construction, the magnitude quantify the resemblance of a generic polygon with a regular p-sided polygon of the same size, while the phase marks the orientation of the polygon. For regular V-sided polygons, provided is an integer multiple of and otherwise. Furthermore, from one can readily compute.
Figure 5b, c shows examples of the functions and for a typical configuration of the SPV. We emphasize that , which, as shown in Armengol-Collado et al., 2023, arises from a p-fold generalization of the classic shape tensor (Aubouy et al., 2003), is solely determined by the positions of the vertices of an individual polygon and, therefore, does not depend on the spatial organization of the neighboring cells. As a consequence, this approach establishes an orientation purely based on cellular shape, thereby eliminating the arbitrariness involved with associating a network of bonds to a planar tessellation, where the latter is not inherent.
The shape function can then be coarse grained at the length scale to construct the shape parameter:
where the is the position of the cth cell, is the Heaviside step function, such that for and 0 otherwise, and is the number of cells within a distance from . As in the case of , the magnitude of reflects the resemblance between a multicelluar cluster and a regular p-sided polygon, while its phase marks the cluster’s global orientation. The outcome of an application of this method to the Voronoi model is illustrated in Figure 6 for . The different patches in panel (a) are regions with uniform , while in panel (b), there are plotted streamlines showing the orientation of the director .
Passive stresses
As explained in the main text, the passive contribution to the stress tensor is given by , where, as demonstrated in Giomi et al., 2022b
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. , and taking , 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 and denotes the orientational stiffness of the p-atic phase, related to the order parameter stiffness by
and is the two-dimensional antisymmetric tensor, with and .
Linear fluctuating hydrodynamics
To compute the structure factor, we follow Ramaswamy et al., 2003 and augment Equation 2b, Equation 2c with short-ranged correlated noise field. Then calling and the nematic and hexatic fluctuating orientation fields and linearizing the hydrodynamic equations about the homogeneous and stationary solutions, and , gives
where , , and indicate a small departure from the homogeneous and stationary configurations of the fields , , and , , , and and are short-ranged correlated noise fields: i.e.
The velocity field , 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 and are the body forces resulting from the passive and active stresses, respectively. The quantity 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 is either conservative or null, from which
with and the case of noiseless translational dynamics, corresponding to Figure 3 in the main text, is recovered in the limit . The pressure , in turn, can be related to the density by a linear equation of state of the form
with 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 is given by
and the functions , with , are effective noise fields whose correlation functions are given by
where the functions 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 Equation 16b, Equation 16c, this effect disappears at the small (i.e. ) and large (i.e. ) 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 , can be calculated from the correlation function
To compute the left-hand side of Equation 29 one can solve Equation 24 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 when . The effect of the active currents is instead accounted for by the second and third terms on the right-hand side of Equation 31, which can be cast in the general form
where
The integral over can be derived using the residue theorem upon computing the roots of the complex third-order polynomial . To make progress, we express
where , , and are given by
The integrand on the right-hand side of Equation 33 has, therefore, three pairs of purely imaginary poles: i.e. , , and . 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 depend on the individual components of the wave vector – i.e. and – this is an artefact of linearizing the hydrodynamic equations about a specific orientation (i.e. 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 . To remove the fictitious angular dependence, one can either linearize Equation 2a about a generic pair of angles, and , and then use these to calculate a circular average – i.e. – or, more simply, by orienting so to cancel the directional dependence. Thus, taking gives a simpler expression of the matrix . 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 and in Equation 8. That is
Notice that, while both orientational and translation noise affect the amplitude of density fluctuations at small length scales, where , translational noise becomes unimportant at the large scale, where . Furthermore, as long as viscous dissipation is at play, switching off translational noise (i.e. ) does not alter the scaling behavior of the structure factor at neither range of length scales. Taking the dry limit (i.e. and ) 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, in the case of purely rotational noise and in the presence of rototranslational noise. The coefficients and can be computed as in the viscous case, to give
Numerical methods
The Voronoi model
In the self-propelled Voronoi model (Bi et al., 2016) a confluent cell layer is approximated as a Voronoi tessellation of the plane. Each cell is characterized by the position of its center, with , and a velocity , with a constant speed and 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 is an energy function involving exclusively geometrical quantities, such as the area and the perimeter of each cell: i.e.
with , , , and 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 ) and the effective interfacial tension caused by the cell–cell adhesion and the cortical tension (both linear in ) (Farhadifar et al., 2007). The constants and represent, respectively, the preferred area and perimeter of each cell. The quantity , on the other hand, is a random number with zero mean and correlation function
with a rotational diffusion coefficient. To make progress, we next introduce the following dimensionless numbers: the shape index , which accounts for the spontaneous degree of acircularity of individual cells (Bi et al., 2016), and the Péclet number , 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 with periodic boundary conditions. At , the centroids 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 , , , and , where is the time-step used for the integration, and the average density of particles . We vary the Péclet number in the range . 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 at constant density. The density structure factor (light green circles) in Figure 3a was obtained, in particular, with .
The MPF model
The MPF model is a continuous model where each cell is described by a concentration field with and 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 , where the free energy density is given by
Here, and are material parameters which can be used to tune the surface tension and the interfacial thickness of isolated cells and thermodynamically favor spherical cell shapes. The constant captures the repulsion between cells. The concentration field is large (i.e. ) inside the cells and zero outside. The contribution proportional to in the free energy enforces cell incompressibility whose nominal radius is given by . The relaxational dynamics of the field is governed by the Allen–Cahn equation
where 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 in Equation 46 is the mobility measuring the relevance of thermodynamic relaxation with respect to non-equlibrium cell migration. The dimensionless parameters of the model are the Péclet number and the cell deformability .
The system of partial differential equations, Equation 46, is solved with a finite-difference approach through a predictor–corrector finite-difference Euler scheme implementing second-order stencil for space derivatives (Carenza et al., 2019). The C-code implemented for numerical integration is parallelized by means of Message Passage Interface (MPI). We consider systems of cells in a square domain of grid points. Model parameters in simulation units are as follows: , , , , , , , , being the time-step used to integrate Equation 46. We vary the speed of self-propulsion in the range . In terms of dimensionless parameters this corresponds to having and ranging between 0 and 2.30. The timescale of cell motility with respect to the timescale of elastic relaxation driven by surface tension ranges between 0 and 0.625. Moreover, the nominal packing fraction is , while the ratio between the interface thickness and the nominal radius . The density structure factor (dark green triangles) in Figure 3a was obtained with .
Numerical method for integration of the hydrodynamic equations
Equation 2a has been integrated by means of a hybrid lattice Boltzmann (LB) method, in which Equation (2b) is solved through a predictor–corrector LB algorithm and the remaining equations via a predictor–corrector finite-difference Euler approach, with a first-order upwind scheme and second-order accurate stencils for the computation of spacial derivatives (Carenza et al., 2019). The code has been parallelized by means of MPI, by dividing the computational domain in slices and by implementing the ghost-cell method to compute derivatives on the boundary of the computational subdomains. Runs have been performed using 64 CPUs in two-dimensional geometries, on a computational box of size 2562 and 5122, for at least LB iterations (corresponding to ∼21 and ∼84 days of CPU-time, respectively, for the smaller and larger computational boxes). Periodic boundary conditions have been imposed. The director fields (for both and ) have been randomly initialized. The initial density field is assumed to be uniform with everywhere. The model parameters in simulations units are as follows: , , , , , , , , , . Nematic activity has been varied in the range and hexatic activity in the range . We set the active parameters and . The density structure factor (continuous black line) in Figure 3a was obtained with and .
The coherence length of the nematic and hexatic liquid crystal can be expressed as the for both , where 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 and ranges between for and for . Conversely, for hexatics and ranges up to for . To compare the results of the hydrodynamics simulations with the discrete models in Figure 3a, we choose and , with 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 given in Equation 9, reflecting the amount of cross-correlation in their fluctuations. Here, and , while the fluctuating fields and represent again the local nematic and hexatic orientations, respectively. Averaging and over the scale of a volume element, yields the order complex parameters and , 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
from which
This approximation holds when the relative fluctuation of the p-atic phase is sufficiently small, so that
consistent with the standard definition of p-atic order parameter. Thus, in particular, and , whereas and . This allows to write , 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 long-ranged ordered liquid crystal, where also and are expected to vary in space. These spatial variations, however, occur on length scales comparable with the system size and, as long as this is much larger than any of the intrinsic length scales entailed in Equation 2a, are negligible for the purpose of this calculation. To compute , 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 in Equation 6b. For , the hexatic and nematic directors are energetically favored to be parallel, so that . Conversely, when , the hexatic and nematic directors are preferentially tilted by , hence . 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 values. Thus, assuming and expanding Equation 2c about , gives
where, as in the previous sections, we have set and and introduced the Gaussian noise fields and , having vanishing mean and finite variance. Unlike the active case, however, at equilibrium the latter is related to the environmental temperature by the fluctuation–dissipation theorem. This implies
where , with the orientational stiffness defined in Equation 15, is the rotational viscosity of the associated p-atic phase. Equation 53a can now be decoupled and used to compute the correlation function . For simplicity, here we set , , and . With this choice, taking
gives, after simple algebraic manipulations
where and . Moreover, using Equation (54), one finds
where . Equation 56a can now be solved in Fourier space and real time to give
where the hat indicates Fourier transformation and
where and . The calculation of the cross-correlation function is now reduced to calculating the autocorrelation functions of the fields and . Specifically
where
and we have made use of Equation (54) to demonstrate that . The non-vanishing correlation functions, on the other hand, can be expressed as
where is a short-distance cut-off and is the finite-time 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 non-interacting scalar fields and . Now, in the case of the ‘massive’ field , the Fourier integral in Equation (65) converges to
in the range . Here, is a modified Bessel function of the second kind, whose asymptotic expansion at short and long distances is given by
with the Euler–Mascheroni constant. In the case of the ‘massless‘ field , on the other hand, the Fourier integral diverges in the infrared, but the correlation function can still be computed as the Laplacian Green function on an infinite domain punctured by a hole of radius at the origin. Thus
Combining this with Equations (67) and (69) yields the following expression for the correlation function
where . Finally, using Equation (50) and the asymptotic expansions of Equation (68) gives the following expression for the cross-correlation function
where is an instance of the generic non-universal exponent
in the specific case and . Lastly, when , the same procedure can be carried out by expanding Equation (2c) about and taking and , 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/s41567-023-02179-0
-
Geometric constraints during epithelial jammingNature Physics 14:613–620.https://doi.org/10.1038/s41567-018-0089-9
-
A texture tensor to quantify deformationsGranular Matter 5:67–70.https://doi.org/10.1007/s10035-003-0126-x
-
A density-independent rigidity transition in biological tissuesNature Physics 11:1074–1079.https://doi.org/10.1038/nphys3471
-
Motility-driven 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 herring-bone 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/i2019-11843-6
-
Dry aligning dilute active matterAnnual Review of Condensed Matter Physics 11:189–212.https://doi.org/10.1146/annurev-conmatphys-031119-050752
-
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 hexatic-nematic XY modelPhysical Review X 12:011043.https://doi.org/10.1103/PhysRevX.12.011043
-
Thermally driven order-disorder transition in two-dimensional 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 p-atic liquid crystalsPhysical Review. E 106:024701.https://doi.org/10.1103/PhysRevE.106.024701
-
Long-ranged order and flow alignment in sheared p-atic 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/i2007-10298-8
-
Instabilities and geometry of growing tissuesPhysical Review Letters 129:048102.https://doi.org/10.1103/PhysRevLett.129.048102
-
Collective cancer invasion forms an integrin-dependent radioresistant nicheThe Journal of Experimental Medicine 217:e20181184.https://doi.org/10.1084/jem.20181184
-
Description of cellular patterns by dirichlet domains: the two-dimensional caseJournal of Theoretical Biology 72:523–543.https://doi.org/10.1016/0022-5193(78)90315-6
-
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: stress-coordinated cellular constriction chainsJournal of Physics. Condensed Matter 28:414021.https://doi.org/10.1088/0953-8984/28/41/414021
-
Role of cell deformability in the two-dimensional melting of biological tissuesPhysical Review Materials 2:045602.https://doi.org/10.1103/PhysRevMaterials.2.045602
-
Solid-liquid 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/i2007-10300-7
-
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/s41567-018-0279-5
-
Active dynamics of tissue shear flowNew Journal of Physics 19:033006.https://doi.org/10.1088/1367-2630/aa5756
-
Active nematics on a substrate: giant number fluctuations and long-time tailsEurophysics Letters 62:196–202.https://doi.org/10.1209/epl/i2003-00346-7
-
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 liquid-crystal 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
-
Low-noise phase of a two-dimensional active nematic systemPhysical Review. E 97:012707.https://doi.org/10.1103/PhysRevE.97.012707
-
Light-scattering study of bond orientational order in a tilted hexatic liquid-crystal 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/i2004-10501-2
-
Epithelial-to-mesenchymal transition in cancer: complexity and opportunitiesFrontiers of Medicine 12:361–373.https://doi.org/10.1007/s11684-018-0656-6
-
Dynamics of two-dimensional 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
- Josep-Maria Armengol-Collado
The funders had no role in study design, data collection, and interpretation, or the decision to submit the work for publication.
Acknowledgements
We are indebted with Massimo Pica Ciamarra and David Nelson for insightful discussions. This work is supported by the ERC-CoG grant HexaTissue (LNC and LG) and by Netherlands Organization for Scientific Research (NWO/OCW) as part of the research program 'The active matter physics of collective metastasis' with project number Science-XL 2019.022 (J-M A-C). Part of this work was carried out on the Dutch national e-infrastructure with the support of SURF through the NWO Grant 2021.028 for computational time.
Copyright
© 2024, Armengol-Collado, Carenza et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics
-
- 951
- views
-
- 214
- downloads
-
- 4
- 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
-
- Physics of Living Systems
Many proteins have been recently shown to undergo a process of phase separation that leads to the formation of biomolecular condensates. Intriguingly, it has been observed that some of these proteins form dense droplets of sizeable dimensions already below the critical concentration, which is the concentration at which phase separation occurs. To understand this phenomenon, which is not readily compatible with classical nucleation theory, we investigated the properties of the droplet size distributions as a function of protein concentration. We found that these distributions can be described by a scale-invariant log-normal function with an average that increases progressively as the concentration approaches the critical concentration from below. The results of this scaling analysis suggest the existence of a universal behaviour independent of the sequences and structures of the proteins undergoing phase separation. While we refrain from proposing a theoretical model here, we suggest that any model of protein phase separation should predict the scaling exponents that we reported here from the fitting of experimental measurements of droplet size distributions. Furthermore, based on these observations, we show that it is possible to use the scale invariance to estimate the critical concentration for protein phase separation.
-
- Computational and Systems Biology
- Physics of Living Systems
Explaining biodiversity is a fundamental issue in ecology. A long-standing puzzle lies in the paradox of the plankton: many species of plankton feeding on a limited variety of resources coexist, apparently flouting the competitive exclusion principle (CEP), which holds that the number of predator (consumer) species cannot exceed that of the resources at a steady state. Here, we present a mechanistic model and demonstrate that intraspecific interference among the consumers enables a plethora of consumer species to coexist at constant population densities with only one or a handful of resource species. This facilitated biodiversity is resistant to stochasticity, either with the stochastic simulation algorithm or individual-based modeling. Our model naturally explains the classical experiments that invalidate the CEP, quantitatively illustrates the universal S-shaped pattern of the rank-abundance curves across a wide range of ecological communities, and can be broadly used to resolve the mystery of biodiversity in many natural ecosystems.