Introduction

The importance of orientational order in biological systems is becoming increasingly clear, as it plays a critical role in processes such as tissue morphogenesis, collective cell motion, and cellular extrusion. Orientational order results from the shapes of cells and their alignments with neighboring cells. Disruptions in this order, known as topological defects, are often linked to key biological events. For instance, defects—points or lines where the order breaks down—can drive cell extrusion (Saw et al., 2017; Monfared et al., 2023) or trigger morphological changes in tissues (Maroudas-Sacks et al., 2021; Ravichandran et al., 2024). Orientational order is linked to liquid crystal theories and should here be interpreted in a broad sense. Recent evidence extends beyond nematic order, characterized by symmetry under 180° = 2π/2 rotation, to higher-order symmetries such as tetratic order (90° = 2π/4), hexatic order (60° = 2π/6), and even general p-atic orders (2π/p, p being an integer). In biological contexts, nematic order (p = 2) has been widely studied in epithelial tissues (Duclos et al., 2017; Saw et al., 2017; Kawaguchi et al., 2017), linking defects to cellular behaviors and tissue organization. More complex orders, such as tetratic order (p = 4) (Cislo et al., 2023) and hexatic oder (p = 6) (Li and Ciamarra, 2018; Durand and Heu, 2019; Pasupalak et al., 2020; Armengol-Collado et al., 2023; Eckert et al., 2023), have also been observed in experimental systems, offering new perspectives on how cells self-organize and respond to mechanical cues. Understanding and quantifying orientational order and the corresponding liquid crystal theory are therefore essential for unraveling the physical mechanisms that govern biological dynamics.

In physics, p-atic liquid crystals illustrate how particle shape and symmetry influence phase behavior, with seminal works dating back to Onsager’s theories (Onsager, 1949). Most prominently, hexatic order (p = 6) has been postulated and found by experiments and simulations as an intermediate state between crystalline solid and isotropic liquid in (Halperin and Nelson, 1978; Nelson and Halperin, 1979; Murray and Van Winkle, 1987; Bladon and Frenkel, 1995; Zahn et al., 1999; Gasser et al., 2010; Bernard and Krauth, 2011). Other examples are colloidal systems for triadic platelets (Bowick et al., 2017) or cubes (Wojciechowski and Frenkel, 2004), which lead to p-atic order with p = 3 and p = 4, respectively. Even pentatic (p = 5) and heptatic (p = 7) liquid crystals have been engineered (Wang and Mason, 2018; Yu and Mason, 2023). The corresponding liquid crystal theories for p-atic order have only recently been defined (Giomi et al., 2022; Krommydas et al., 2023) and can also be used in biological contexts. However, while in colloidal systems particle shapes remain fixed, in biological tissues, cells are dynamic: their shapes are irregular, variable, and influenced by internal and external forces. These unique properties make quantifying p-atic order in tissues significantly more challenging, as quantification of the cell shapes is required. We will demonstrate that existing methods, such as bond-orientational order (Nelson and Halperin, 1979) or polygonal shape analysis (Armengol-Collado et al., 2023), often fail to capture the nuances of irregular cell shapes, which has severe consequences on the definition of p-atic order.

To address this, we adapt Minkowski tensors (Mecke, 2000; Schröder-Turk et al., 2013)—rigorous mathematical tools for shape analysis—to quantify p-atic order in cell monolayers. Minkowski tensors provide robust and sensitive measures of shape anisotropy and orientation, accommodating both smooth and polygonal shapes while remaining resilient to small perturbations. By applying these tools, we re-examine previous conclusions about p-atic order in epithelial tissues and demonstrate that certain widely accepted results require reconsideration.

Our study leverages a combination of experimental and computational approaches. Experimentally, we analyze confluent monolayers of MDCK (Madin-Darby Canine Kidney) cells, comparing wild-type tissues with E-cadherin knockout variants to explore the role of cell-cell adhesion in orientational order. Figure 1 provide a snapshot of the considered monolayer of wild-type MDCK cells, together with their shape classification by Minkowski tensors. The corresponding statistical data and probability distributions of these quantities are shown in Figure 2. These data already indicate no clear preference of a specific p-atic order and the necessity for robust and sensitive shape classification. Computationally, we employ two complementary models for cell monolayers: the active vertex model and the multiphase field model. These approaches allow us to systematically vary parameters such as cell activity and mechanical properties, providing a comprehensive view of how p-atic order emerges in different contexts.

Figure 1.

Figure 2.

By combining these robust mathematical tools, experimental insights, and computational models, we establish a reliable framework for quantifying p-atic order in biological tissues. This framework not only challenges previous findings but also opens new pathways for understanding the role of orientational symmetries in tissue mechanics and development.

Methods and materials

Minkowski tensors

Essential properties of the geometry of a two-dimensional object are summarized by scalar-valued size measures, so called Minkowski functionals. They are defined by

where 𝒞 is the smooth two-dimensional object, with contour ∂𝒞 and ℋ denotes the curvature of the contour ∂𝒞 (see Figure 3 for a schematic description), as reviewed in Mecke (2000). They describe (up to field-specific prefactors) the area, the perimeter and a curvature-weighted integral of the contour, respectively. Minkowski functionals are also known as intrinsic volumes. For convex objects, they have been shown to be continuous and invariant to translations and rotations. Due to these properties and Hadwiger’s characterization theorem (Hadwiger, 1975), they are natural size descriptors that provide essential and complete information about invariant geometric features. However, these properties also set limits to their use as shape descriptors. Due to rotation invariance, they are unable to capture the orientation of the shape. To describe more complex shape information the scalar-valued Minkowski functionals are extended to a set of tensor-valued descriptors, known as Minkowski tensors. These objects have been investigated both in mathematical (Alesker, 1999; Hug et al., 2008, 2007; McMullen, 1997) and in physical literature (Schröder-Turk et al., 2010a; Kapfer et al., 2012). Using tensor products of the position vectors r and the normal vectors n of the contour ∂𝒞, defined as r_a ⊙ n_b = r ⊙ … ⊙ r ⊙ n ⊙ … ⊙ n, the first considered a times and the last b times with ⊙ denoting the symmetrized tensor product. The Minkowski tensors are defined as

Figure 3.

For for j = 0, 1, 2, so the Minkowski tensors are extensions of the Minkowski functionals. In analogy to Minkowski functionals also Minkowski tensors have been shown to be continuous for convex objects. Also, analogies to Hadwiger’s characterization theorem exist (Alesker, 1999). However, for our purpose, the essential property is the continuity. It makes Minkowski functionals and Minkowski tensors a robust measure as small shape changes lead to small changes in the shape descriptor. Due to this property, Minkowski tensors have been successfully used as shape descriptors in different fields, e.g. in materials science as a robust measure of anisotropy in porous (Schröder-Turk et al., 2010b) and granular material (Schröder-Turk et al., 2013), in astrophysics to describe morphology of galaxies (Beisbart et al., 2002), and in biology to distinguish shapes of different types of neuronal cell networks (Beisbart et al., 2006) and to determine the direction of elongation in multiphase field models for epithelial tissue (Mueller et al., 2019; Wenzel and Voigt, 2021; Happel and Voigt, 2024). However, the mentioned examples exclusively use lower-rank Minkowski tensors a + b ≤ 2. Higher rank tensors have not been considered in such applications but the theory also guarantees robust description of p-atic orientation for p = 3, 4, 5, 6, …. These results hold for smooth contours but also can be extended to polyconvex shapes. Furthermore, known counter-examples for non-convex shapes have very little to no relevance in applications (Kapfer, 2012).

Even if all Minkowski tensors carry important geometric information, one type is particularly interesting and will be considered in the following:

In this case, the symmetrized tensor product agrees with the classical tensor product. If applied to polygonal shapes the Minkowski tensors are related to the Minkowski problem for convex polytops (i.e., generalizations of three-dimensional polyhedra to an arbitrary number of dimensions), which states that convex polytops are uniquely described by the outer normals of the edges and the length of the corresponding edge (Minkowski, 1897). Several generalizations of this result exist (Klain, 2004; Schneider, 2013), which makes the normal vectors a preferable quantity to describe shapes. However, Minkowski tensors contain redundant information, which asks for an irreducible representation. This can be achieved by decomposing the tensor with respect to the rotation group SO(2) (Kapfer, 2012; Mickel et al., 2013; Klatt et al., 2022). Following this approach, one can write

with Ψ_𝒞 (u) being proportional to the probability density of the normal vectors. Identifying u ∈ S¹ by the angle θ between u and the x−axis allows to write

The Fourier coefficients ψ_p (𝒞) are the irreducible representations of and can be written as

where θ_n is the orientation of the outward pointing normal n. The phase of the complex number ψ_p(𝒞) contains information about the preferred orientation and the absolute value |ψ_p(𝒞)| is a scalar index. We thus define

with ℑ_ψp and ℜ_ψp the imaginary and real part of ψ_p, respectively, and . We have q_p(𝒞) ∈ [0, 1] quantifying the strength of p-atic order. Figure 4 illustrates the concept for polygonal and smooth shapes and Figure 5 illustrates the calculation of q_p and ϑ_p at the example of an equilateral triangle. Note that the constant factor of 1/(2π) in Equation 6 does not influence the result of Equation 7 and is therefore disregarded in Figure 5. Alternative derivations of Equation 7 consider higher order trace-less tensors (Virga, 2015; Giomi et al., 2022; Armengol-Collado et al., 2023) generated by the normal n. In this case ϑ_p(𝒞) (orange triatic director in Figure 5 right) correspond to the eigenvector to the negative eigenvalue and ϑ (𝒞)− π/p (gray triatic director in Figure 5 right) to the eigenvector to the positive eigenvalue. For p = 2 this tensor based approach corresponds to the structure tensor (Mueller et al., 2019).

Figure 4.

Figure 5.

As pointed out in the introduction, p-atic orders have already been identified in biological tissue and model systems. However, the way to quantify rotational symmetry in these works differs from the Minkowski tensors. In Armengol-Collado et al. (2023) the cell contour is approximated by a polygon with V vertices having coordinates x_i. The considered polygonal shape analysis is based on the shape function

with r_i = x_i − x_𝒞 being the vector from the center of mass x_𝒞 and θ_i being the orientation of the i-th vertex of the polygon with respect to its center of mass. As the center of mass can be computed by the only input data are the coordinates of the vertices of the V-sides polygon. Equation 8 captures the degree of regularity of the polygon by its amplitude 0 ≤ |γ_p| ≤ 1 and its orientation ϑ_p which follows as

with ℑ_γp and ℜ_γp the imaginary and real part of γ_p, respectively. Instead of the normals n, the descriptor is based on the position vector r. One might be tempted to think that Equation 9 is related to the irreducible representation of . This is, however, not the case as the weighting in Equation 9 is done with respect to the magnitude of r, while in the weighting is done with respect to the contour integral. Furthermore, in contrast to the Minkowski tensors, Equation 9 takes into account only the vectors at the vertices and not the vectors along the full cell contour. While this seems to be a technical detail, it has severe consequences. The theoretical basis, which guarantees continuity, no longer holds. This leads to unstable behaviour, as illustrated in Figure 6. While the shape of the polygons is almost identical, there is a jump in γ₃ and γ₄ as we go from four to three vertices. The weighting according to the magnitude cannot cure this, as the magnitude of the shrinking vector is far from zero shortly before this vertex vanishes. In the context of cellular systems, such deformations are a regular occurrence rather than an artificially constructed test case. We will demonstrate the impact while discussing the results and recommend only using robust descriptors such as the Minkowski tensors or their irreducible representations.

Figure 6.

We further note that bond order parameters , which consider the connection between cells, have also been used as shape descriptors (Li and Ciamarra, 2018; Durand and Heu, 2019; Pasupalak et al., 2020). Introduced in Nelson and Halperin (1979), these parameters are

with B denoting the number of bonds and θ_b the orientation of the bond b. In the context of mono-layer tissues, B is understood as the number of neighbors, and θ_b as the orientation of the connection of the center of mass of the current cell with the center of mass of the neighbor b (Loewe et al., 2020; Monfared et al., 2023). So cells with the same contour but different neighbor relations lead to different bond order parameters. This characteristic should already disqualify these measures as shape descriptors. However, as discussed in detail in Mickel et al. (2013), they are not robust even for the task they are designed for. We therefore do not discuss them further.

Simulation methods

We consider two modeling approaches, an active vertex model and a multiphase field model. Both have been proven to capture various generic properties of epithelial tissue (Hakim and Silberzan, 2017; Alert and Trepat, 2020; Moure and Gomez, 2021). The key aspects of formulations are out-lined in Box 1 and Box 2. They refer to (Koride et al., 2018; Das et al., 2021; Killeen et al., 2022) and (Wenzel and Voigt, 2021; Jain et al., 2023, 2024b), respectively. Parameters used within these models are listed in Appendix 1 Table 1 and Appendix 1 Table 2.

The initial configuration for the active vertex model simulations was built by placing N points randomly in a periodic simulation box of length L × L and ensuring that no two points were within a distance less than r_cut from each other. The points were used to create a Voronoi tiling. Typically, such a tiling had cells of very irregular shapes. To make cell shapes more uniform, the centroids of each tile were computed and used as seeds for a new Voronoi tiling. This process was iterated until it converted within the tolerance of 10⁻⁵, resulting in a so-called well-centered Voronoi tiling with cells of random shapes but similar sizes. Initial directions of polarity vectors were chosen at random from a uniform distribution. The initial configuration for the multiphase field model simulations considers a regular arrangement of N cells of equal size in a periodic simulation box of length L × L with randomly chosen directions of self-propulsion and simulating for several time steps until the cell shapes appear sufficiently irregular. For both models we consider 100 cells and analyse time instances of the evolution.

Experimental setup

Cell culture

Madin-Darby canine kidney (MDCK) cells were cultured in DMEM (DMEM, low glucose, GlutaMAX™ Supplement, pyruvate) supplemented with 10 % fetal bovine serum (FBS; Gibco) and 100 U/mL peni-cillin/streptomycin (Gibco) at 37 ^°C with 5 % CO₂. The cell line was tested mycoplasma.

Monolayer preparation

Cells were seeded on glass-bottom dishes (Mattek) pretreated with 10 µg/mL fibronectin human plasma in phosphate buffered saline (PBS, pH 7.4; Gibco). Fibronectin was incubated for 30 min at 37 ^°C. The initial cell seeding density was sparse. They were imaged approximately 24 hours after, when a confluent monolayer was formed.

Box 1. Active vertex model

Box 1—figure 1.

The active vertex model is based on models discussed in (Koride et al., 2018; Das et al., 2021; Killeen et al., 2022). A confluent epithelial tissue is appreciated as a two-dimensional tiling of a plane with the elastic energy given as (Honda, 1983; Farhadifar et al., 2007; Honda and Nagai, 2022)

where A_C and P_C are, respectively, area and perimeter of cell C, A₀ and P₀ are, respectively, preferred area and perimeter, and K_A and K_P are, respectively, area and perimeter moduli. For simplicity, it is assumed that all cells have the same values of A₀, P₀, K_A, and K_P. The model can be made dimensionless by dividing Equation 11 by ,

where a_C = A_C /A₀ (i.e. is the unit of length), , and . p₀ is called the shape index, and it has been shown to play a central role in determining if the tissue is in a solid or fluid state (Bi et al., 2015).

One can use Equation 11 to find the mechanical force on a vertex i as , where is the gradient with respect to the position r_i of the vertex. The expression for F_i is given in terms of the position of the vertex i and its immediate neighbours (Tong et al., 2023), which makes it fast to compute. Furthermore, cells move on a substrate by being noisily self-propelled along the direction of their planar polarity described by a unit-length vector from which the polarity direction at the vertex is defined, see Box 1 - figure 1 b), where z_i is the number of neighbours of vertex i. It is convenient to write , where α_i is an angle with the x−axis of the simulation box. The equations of motion for vertex i is a force balance between active, elastic, and frictional forces and are given as

where v₀ is the magnitude of the self-propulsion velocity defined as f₀/γ_fr. Furthermore the overdot denotes the time derivative, γ_fr is the friction coefficient, f₀ is the magnitude of the self-propulsion force (i.e., the activity), ξ_i(t) is Gaussian noise with ⟨ξ_i(t) ⟩ = 0 and ⟨ξ_i(t)ξ_j (t′)⟩ = 2D_rδ_ij δ(t − t′), where D_r is the rotation diffusion constant and ⟨·⟩ is the ensemble average over the noise. Equations of motion can be non-dimensionalised by measuring time in units of t^* = γ_fr/(K A) and force in units of and integrated numerically using the first-order Euler-Maruyama method with timestep τ.

Box 2. Multiphase field model

Box 2—figure 1.

The multiphase field modelling approach follows Wenzel and Voigt (2021); Jain et al. (2023, 2024b). Each cell is described by a scalar phase field variable ϕ_i with i = 1, 2, …, N, where the bulk values ϕ_i ≈ 1 denotes the cell interior, ϕ_i ≈ −1 denotes the cell exterior, and with a diffuse interface of width (ϵ) between them representing the cell boundary. The phase field ϕ_i follows the conservative dynamics

with the free energy functional ℱ = ℱ_CH + ℱ_INT containing the Cahn-Hilliard energy

where the Capillary number Ca is a parameter for tuning the cell deformability, and an interaction energy consisting of repulsive and attractive parts defined as

where In denotes the interaction strength, and a_r and a_a are parameters to tune contribution of repulsion and attraction. The repulsive term penalises overlap of cell interiors, while the attractive part promotes overlap of cell interfaces. The relation to former formulations is discussed in Happel and Voigt (2024).

Each cell is self-propelled, and the cell activity is introduced through an advection term. The cell velocity field is defined as

where v₀ is used to tune the magnitude of activity, e_i = [cos θ_i(t), sin θ_i (t)] is its direction. The migration orientation θ_i(t) evolves diffusively with a drift that aligns to the principal axis of cell’s elongation as , where D_r is the rotational diffusivity, W_i is the Wiener process, β_i(t) is the orientation of the cell elongation and α controls the time scale of this alignment.

The resulting system of partial differential equations is considered on a square domain with periodic boundary conditions and is solved by the finite element method within the toolbox AMDiS (Vey and Voigt, 2007; Witkowski et al., 2015) and the parallelization concept introduced in Praetorius and Voigt (2018) is considered, which allows scaling with the number of cells. We in addition introduce a de Gennes factor in ℱ_CH to ensure ϕ_i ∈ [−1, 1], see Salvalaglio et al. (2021).

Live cell imaging

Samples were imaged using Nikon ECLIPSE Ti microscope equipped with a H201-K-FRAME Okolab chamber, heating system (Okolab) and a CO₂ pump (Okolab) which maintained them at 37 °C and at 5 % CO₂. Phase microscopy images were taken each 10th minute using a 10× NA=0.3 Plan Fluor objective and Andor Neo 5.5 sCMOS camera.

Image analysis

The time-series were xy-drift corrected using Fast4DReg plugin Laine et al. (2019), Pylvänäinen et al. (2022) in FIJI. Cells were segmented and tracked using the python module CellSegmentation-Tracker (https://github.com/simonguld/CellSegmentationTracker), which utilizes both Cellpose Stringer et al. (2021) and Trackmate Tinevez et al. (2017). A Cellpose model was trained by manual segmentation of the phase contrast images.

Results

Quantifying orientational order in biological tissues can be realized by Minkowsky tensors. The orientation ϑ_p and the strength of p-atic order q_p in Equation 7 can be computed for any integer p. These quantities describe how cell shapes align with specific rotational symmetries. A recurring challenge in the field is determining whether q_p values for different p can be meaningfully compared. This question is particularly important for interpreting nematic (q₂) and hexatic (q₆) orders, which describe distinct symmetries but often coexist in biological systems. However, a direct comparison is problematic as q₂ and q₆ reveal significant differences in their ranges: q₂ = 1.0 represents an idealized line shape that is unattainable for cells, while q₆ = 1.0 corresponds to a perfect hexagon observed in experiments and simulations. To further address this, we analyze simulation data from two computational models, examining how deformability and activity—key determinants of tissue mechanics—affect q_p values and their interplay. We employ the active vertex model and the multiphase field model, both of which incorporate activity and deformability. These two parameters are recognized as crucial for capturing the coarse-grained properties of confluent tissues (Jain et al., 2023, 2024b). By leveraging these complementary frameworks, we ensure that our conclusions are robust across model-specific assumptions and parameter regimes. This dual approach allows us to systematically investigate how these models capture tissue behavior, particularly regarding cell shape variability and resulting p-atic orders.

In the active vertex model, deformability is controlled by the target shape index p₀, reflecting the balance between cell-cell adhesion and cortical tension. The multiphase field model, by contrast, encodes deformability through the capillary number Ca, which directly incorporates cortical tension. Although these models differ conceptually, both have been validated in studies of cell monolayer mechanics, including solid-liquid transitions, neighbor exchange dynamics, and stress profiles (Fletcher et al., 2014; Alt et al., 2017; Li et al., 2019; Balasubramaniam et al., 2021; Wenzel and Voigt, 2021; Sknepnek et al., 2023; Melo et al., 2023). Importantly, the active vertex model enforces confluency by design, while the multiphase field model approximates it through dynamic cell-cell interactions, leading to differences in topological rearrangements and cell shape representations.

We vary activity strength (v₀), the shape index (p₀), and the capillary number (Ca), ensuring all parameter combinations remain in the fluid regime. Fluidity was confirmed using metrics such as the mean square displacement (MSD) (Loewe et al., 2020), neighbor number variance (Wenzel and Voigt, 2021), and the self-intermediate scattering function (Bi et al., 2016). Using Minkowski tensors, we computed ϑ_p and q_p values as described in Equation 7, focusing on p = 2 and p = 6 due to their biological relevance. In solid-like states, such as hexagonal cellular arrangements, hexatic order (p = 6) is the appropriate descriptor, as hexagons exhibit 60° = 2π/6 rotational symmetry. Hexatic order has been extensively studied in biological systems (Li and Ciamarra, 2018; Durand and Heu, 2019; Pasupalak et al., 2020), and its importance in cellular organization is well-documented. In fluid-like states, cells often become elongated, exhibiting nematic order (p = 2), characterized by 180° = 2π/2 rotational symmetry. Elongated cell shapes are also common in specific cell types, such as fibroblasts. Quantifying nematic order using q₂ is particularly relevant, as nematic configurations and defects have been linked to critical biological processes, including hydra morphogenesis (Maroudas-Sacks et al., 2021) and cell death (Saw et al., 2017).

Independence of q₂ and q₆

Recent studies suggest a hexatic-nematic crossover in epithelial tissues, where hexatic order dominates at small scales and nematic order prevails at larger scales (Eckert et al., 2023; Armengol-Collado et al., 2023, 2024). To explore this relationship, we examine the distribution of (q₂, q₆) values across deformability-activity parameter pairs in both models (Figure 7). Despite structural differences between the models, both show consistent trends. Closer examination shows that in near-solid regimes (Figure 7 lower left), q₂ and q₆ values cluster tightly due to restricted shape fluctuations. However, even in this regime, small q₂ values can correspond to either small or large q₆ values, and vice versa. In more fluid-like regimes, with higher activity and higher deformability (Figure 7 upper right), q₂ and q₆ values become highly scattered. Each q₂ value spans a broad range of q₆ values, and vice versa, demonstrating their independence. These findings underscore that q_p values for different p measure distinct aspects of cell shape anisotropy and cannot be meaningfully compared directly.

Figure 7.

This conclusion also holds on larger length scales. In order to test the hypothesis of a hexaticnematic crossover at larger length scales (Armengol-Collado et al., 2023) we consider a coarse-graining of q₂ and q₆ across spatial domains. Details are described in Appendix 2. Appendix 2 - Figure 13 and Appendix 2 - Figure 14 show the results for the active vertex model and the multiphase field model, respectively. The results show no consistent trends indicative of a crossover, regardless of the considered model and coarse-graining radius. Instead, the apparent convergence of Q₂ and Q₆ values (coarse-grained values of q₂ and q₆) with increasing activity or deformability simply reflects the same trends observed at the individual cell scale.

Dependence of q_p on activity and deformability

The previous analysis demonstrated that comparing q_p across different values of p is not meaningful. However, fixing p and exploring q_p under varying physical conditions offers valuable insights into how tissue properties such as activity and deformability influence cell shape and orientational order. We again focus on nematic (p = 2) and hexatic (p = 6), as shown in Figure 8 and Figure 9. Additional results for p = 3, 4, 5 are provided in Figure 8—figure Supplement 1–Figure 8—figure Supplement 3 (fixed activity) and Figure 9—figure Supplement 1–Figure 9—figure Supplement 3 (fixed deformability). In all plots all cells and all time steps are considered. In Figure 8, we vary deformability (p₀ in the active vertex model and Ca in the multiphase field model) while keeping the activity v₀ constant. Results are presented for the active vertex model (left column) and the multiphase field model (right column). Activity increases from bottom to top rows. Both models show qualitatively similar trends in the probability density functions (PDFs) of q₂ and q₆. For p = 6, increasing deformability shifts the PDF of q₆ to the left, indicating lower mean values. In contrast, for p = 2, higher deformability leads to higher q₂ values, reflected in a rightward shift of the PDF. This trend is confirmed by the mean values, shown as a function of deformability (p₀ and Ca, respectively) in the inlets. Additionally, the PDFs broaden with increasing deformability, and this effect is more pronounced at lower activity levels. A notable difference between the models is the range of q_p values. In the multiphase field model, q_p rarely exceeds 0.8 due to the smoother, more rounded cell shapes, whereas the active vertex model often produces higher q_p values. In Figure 9, deformability is held constant while activity v₀ is varied. Results are again shown for the active vertex model (left column) and the multiphase field model (right column), with increasing activity indicated by brighter colors. Deformability increases from bottom to top rows. The trends mirror those observed in Figure 8. Increasing activity reduces the mean value of q₆ while increasing q₂, see inlets. The effects of activity are more pronounced at lower deformabilities; at higher deformabilities, differences between parameter regimes diminish. The results agree qualitatively with the ones obtained for the measures |γ₂| and |γ₆| reported in (Armengol-Collado et al., 2023, Extended Data Fig. 3). Note that in the calculation of γ_p every cell is approximated by a polygon. Therefore, also for the multiphase field model the obtained values range up to 1.0. The overall behavior of q₂ and q₆ is summarized in Figure 10, which shows the mean values and as functions of activity and deformability. For both models:

Figure 8.

Figure 9.

Figure 10.

• increases with higher activity or deformability

• decreases with higher activity or deformability.

This behaviour is consistent with a trend towards nematic order in more dynamic regimes and toward hexagonal order in more constrained, less dynamic regimes. The first is further confirmed by recent findings in analysing T1 transitions and their effect on cell shapes Jain et al. (2024a). These studies suggest that cells transiently elongate when they are undergoing T1 transitions. As the number of T1 transitions increases with activity or deformability Jain et al. (2024a) this elongation contributes to the observed behavior. The second is consistent with the emergence of hexagonal arrangements in solid-like states. Corresponding results for p = 3, 4, 5 are shown in Figure 10—figure Supplement 1.

Comparison with experiments

Next, we apply shape descriptor framework to experiments on epithelial cells. We consider a monolayer of Madine-Derby-Canine-Kidney (MDCK) cells. We consider two types of cell monolayers: MDCK wild type cells and MDCK E-cad KO cells, which have the E-cadherin protein removed. These E-cad KO cells have weaker cell-cell adhesion strengths and previous studies have shown that compared to wild-type cells, the E-cad KO cells have lower velocities, same spreading area, higher aspect ratio, and higher stiffness Balasubramaniam et al. (2021). The slower movement of E-cad KO cells compared to the wild-type cells suggest their lower activity, and their higher stiffness can be intepreted as having lower deformability compared to the wild-type cells. However, care must be taken in directly comparing the two cell types and directly comparing to simulation results, as the key parameters varied in simulations—cell activity and deformability—are difficult to manipulate in isolation in experiments. One complicating factor is the strong dependence of cell shape on cell density, as previously reported (Eckert et al., 2023). Higher cell densities generally lead to rounder cell shapes, which we confirmed in two distinct experimental setups (see Figure 11). For both wild-type MDCK cells and E-cadherin knockout variants, increased cell density corresponded to lower mean values of q₂, consistent with rounder cell geometries. For higher values of q_p, however, this density-dependent effect was less pronounced. A plausible explanation lies in the geometry of shapes for larger q_p. Rounded shapes close to a circle can yield relatively high q_p values for certain p, such as the example shown in the third row of Figure 4. Distinguishing such shapes from perfect circles requires highly accurate data, which can be difficult to obtain in experimental imaging, potentially masking density-dependent effects for p = 6.

Figure 11.

Importantly, due to the substantial differences in cell density between wild-type cells and E-cadherin knockout variants, direct comparisons between these two experimental conditions would be unjustifiable. Density influences cell shape too strongly to isolate other effects in these datasets. Despite these limitations, we observed qualitative agreement between experimental and computational results regarding the independence of q₂ and q₆ (see Figure 12). Both measures appear to be unrelated, reinforcing the conclusion that nematic (p = 2) and hexatic (p = 6) order describe distinct aspects of cell shape and alignment. Consequently, comparing q₂ and q₆ directly is not meaningful, whether in simulations nor in experiments.

Figure 12.

Sensitivity of the results on the considered shape descriptor

In order to demonstrate that other methods, such as polygonal shape analysis Armengol-Collado et al. (2023) fail to capture the nuances of irregular cell shapes we also use Equation 8 as a shape descriptor for the considered datasets. For the experimental data and the multiphase field model this first requires an approximation of the cell shape by a polygon. For finding the vertices of the polygon we consider the Voronoi interface method (Saye and Sethian, 2011, 2012), see Appendix 3. The corresponding figures to Figure 1, Figure 2 and Figure 4 are shown in Appendix 3 Figure 16, Appendix 3 Figure 17 and Appendix 3 Figure 15, respectively. While already these results show quantitative differences, the failure of the approach can best be seen in the qualitative difference between Figure 10 and the corresponding results in Appendix 3 Figure 18, which, instead of the monotonic trend for activity and deformabilty, exhibit nonmonotonic trends, peaking at intermediate values of activity or deformability. These divergences highlight the limitations of γ_p in capturing consistent patterns and emphasize the need for stable, reliable shape descriptors like the Minkowski tensors.

For completeness we also provide the corresponding figures to Figure 7 and Figure 13 in Appendix 3 Figure 19 and Appendix 3 Figure 20, respectively. We hereby only consider the active vertex model. The results are similar to the once obtained using the Minkowski tensors and confirm the independence of p-atic order for different symmetries.

Discussion

In this study, we introduced Minkowski tensors as a robust and versatile tool for quantifying p-atic order in multicellular systems, particularly in scenarios involving rounded or irregular cell shapes. By applying this framework to extensive datasets from two distinct computational models—the active vertex model and the multiphase field model—we identified universal trends: increasing activity and deformability of the cells enhance nematic order (p = 2) while diminishing hexatic order (p = 6). The consistency of these findings across two models, despite their inherent differences, underscores the generality and relevance of our results.

While various shape characterization methods, such as the bond order parameter (Loewe et al., 2020; Monfared et al., 2023) and the shape function γ_p (Armengol-Collado et al., 2023), have been explored in the literature, we demonstrated that the choice of shape descriptor significantly impacts the conclusions drawn. Applying γ_p to data from the active vertex model revealed qualitative discrepancies: for p = 2, γ_p and q_p behaved similarly, but for p = 6, γ_p exhibited significant changes. Such divergences, together with limited mathematical foundations, highlight the limitations of γ_p in capturing consistent patterns and emphasize the need for stable, reliable shape measures like the Minkowski tensors. In contrast to the bond order parameter or the shape function γ_p the stability of Minkowski tensors can be mathematically justified. We have further demonstrated why bond order parameter or the shape function γ_p fail in certain situations. This finding is not merely a technical nuance; the choice of shape detection method fundamentally influences the interpretation of p-atic order. Consequently, selecting a robust shape descriptor is critical for obtaining meaningful and reproducible insights into orientational order in multicellular systems.

A critical question in the literature has been whether shape measures for different p-atic orders can be directly compared. While some studies have suggested relationships between γ₂ and γ₆ (Armengol-Collado et al., 2023), our results refute this notion. We demonstrated that measures like q₂ and q₆ capture fundamentally distinct aspects of cell shape and alignment, and comparing them directly is both mathematically and biologically misleading. We further tested the hypothesis of a hexatic-nematic crossover at larger length scales by coarse-graining q₂ and q₆ across spatial domains. The results showed no consistent trends indicative of a crossover, regardless of the coarse-graining strategy and the considered model. Instead, the apparent convergence of Q₂ and Q₆ values (coarse-grained values of q₂ and q₆) with increasing activity or deformability simply reflected the same trends observed at the individual cell scale. This reinforces the conclusion that comparisons of p-atic orders across different p lack a robust basis and that claims of crossovers require careful reconsideration.

Our findings suggest that p-atic orders should be studied independently across length scales, as they describe complementary aspects of cellular organization. The coexistence of distinct orientational orders emphasized in different studies—such as nematic (p = 2) (Duclos et al., 2017; Saw et al., 2017; Kawaguchi et al., 2017), tetratic (p = 4) (Cislo et al., 2023), and hexatic (p = 6) (Li and Ciamarra, 2018)—is not contradictory but highlights the rich, multifaceted nature of cellular organization. Rather than searching for a single dominant order, future research should focus on the interplay of different p-atic orders and their associated defects. This suggests to not only consider p-atic liquid crystal theories Giomi et al. (2022) for one p, but combinations of these models for various p’s. Understanding how these orders interact may reveal how they collectively regulate morphogenetic processes.

Connecting p-atic orders to biological function remains a critical avenue for exploration. While the mathematical independence of q₂, q₆, and other shape measures precludes the identification of a universal dominant order, biological systems may exhibit context-dependent preferences. For example, a specific p-atic order might correlate with or drive a key morphogenetic event. Investigating these connections could yield insights into how tissues achieve functional organization and adapt to environmental cues. As such, while it is difficult to speak of dominating orders from a mathematical point of view, there could be a dominating order from a biological point of view, meaning the p-atic order connected to the governing biological process.

Note

This reviewed preprint has been updated to correct the image for Figure 1, which was previously a duplicate of Figure 8 due to an error at the preprint server.

Acknowledgements

We acknowledge fruitful discussions with Björn Böttger and Brendan Tobin. HJ acknowledges funding by the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 945371. RS acknowledges support from the UK Engineering and Physical Sciences Research Council (Award EP/W023946/1). AD acknowledges funding from the Novo Nordisk Foundation (grant No. NNF18SA0035142 and NERD grant No. NNF21OC0068687), Villum Fonden (grant No. 29476), and the European Union (ERC, PhysCoMeT, 101041418). AV acknowledges funding from the German Research Foundation (Award FOR3013 “Vector- and tensor-valued surface PDEs”) and computing resources provided by JSC through MORPH and by ZIH through WIR. Views and opinions expressed are however those of the authors only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them.

Appendix 1

Parameters for the computational models

Appendix 1—table 1.

Appendix 1—table 2.

Appendix 2

Coarse graining results obtained with the Minkowski tensor

For comparison with the literature we are considering a coarse-graining, following closely the strategy used in Armengol-Collado et al. (2023). We therefore regard the coarse-grained strength of p-atic order Q_p = Q_p (x), which is the average of all shape functions (or equivalently ) of cells whose center of mass x_𝒞 lies within a circle with radius R and center x. In a formula this is:

whereby the formula uses that our simulations contain 100 cells. 𝒞_i denotes the i−th cell and Θ denotes the Heaviside step function with Θ(x) = 0 for x > 0 and Θ(x) = 0 otherwise. As in Armengol-Collado et al. (2023) the position x is sampled over a square grid with spacing with A _𝒞 the cell area. For Appendix 2 Figure 13 and Appendix 2Figure 14 we took the average over all grid points and all time instances. The figures show the resulting data for the active vertex model and the multiphase field model, respectively.

Appendix 2—figure 13.

Appendix 2—figure 14.

Appendix 3

Results using polygonal shape analysis

To visualize the similarities and differences between q_p and γ_p we show the analogue of Figure 4 using the Minkowski tensors in Appendix 3 Figure 15 using polygonal shape analysis. For regular polygons, both shape characterization methods behave in the same way as seen in the first row. For irregular shapes, like in the 2nd row, they already show a different behaviour. Note that the rounded shapes from Figure 4 are missing in Appendix 3 Figure 15 as γ_p requires a polygon.

Appendix 3—figure 15.

The influence of the choice of the shape characterization method is not only visible in the values for single shapes, it can also be seen in the mean and standard deviation. To illustrate this, we use γ_p to characterize the shapes from the same experimental data as in Figure 1 and Figure 2. The segmented data consists of smooth contours for each cell. As γ_p only works with vertex coordinates the first step is to identify these. Note that this always leads to the approximation of the cell shape with a polygon and therefore to a simplification of the shape. For finding the vertices we consider the Voronoi interface method (Saye and Sethian, 2011, 2012). Using the so obtained vertices and Equation 8 and Equation 9 leads to the results shown in Appendix 3 Figure 16 and Appendix 3 Figure 17. Compared to Figure 1 and Figure 2, where the Minkowski tensor was used for the same data, we see that the choice of the shape characterization method highly influences the results.

For comparison with the literature we are considering a coarse-graining, following the strategy introduced in Armengol-Collado et al. (2023). We therefore regard the coarse-grained shape function Γ_p = Γ_p(x), which is the average of all shape functions γ_p of cells whose center of mass x_𝒞 lies within a circle with radius R and center x. In a formula this is:

whereby the formula uses that our simulations contain 100 cells. 𝒞_i denotes the i−th cell and Θ denotes the Heaviside step function with Θ(x) = 0 for x > 0 and Θ(x) = 0 otherwise. As in Armengol-Collado et al. (2023) the position x is sampled over a square grid with spacing with A𝒞 the cell area. For Appendix 3 Figure 20 we took the average of Γ_p over all grid points and all time instances. The figure shows the resulting data for the active vertex model.

Appendix 3—figure 16.

Appendix 3—figure 17.

Appendix 3—figure 18.

Appendix 3—figure 19.

Appendix 3—figure 20.

Figure 8—figure supplement 1.

Figure 8—figure supplement 2.

Figure 8—figure supplement 3.

Figure 9—figure supplement 1.

Figure 9—figure supplement 2.

Figure 9—figure supplement 3.

Figure 10—figure supplement 1.