Quantification of nematic cell polarity in three-dimensional tissues

How epithelial cells coordinate their polarity to form functional tissues is an open question in cell biology. Here, we characterize a unique type of polarity found in liver tissue, nematic cell polarity, which is different from vectorial cell polarity in simple, sheet-like epithelia. We propose a conceptual and algorithmic framework to characterize complex patterns of polarity proteins on the surface of a cell in terms of a multipole expansion. To rigorously quantify previously observed tissue-level patterns of nematic cell polarity (Morales-Navarrete et al., eLife 2019), we introduce the concept of co-orientational order parameters, which generalize the known biaxial order parameters of the theory of liquid crystals. Applying these concepts to three-dimensional reconstructions of single cells from high-resolution imaging data of mouse liver tissue, we show that the axes of nematic cell polarity of hepatocytes exhibit local coordination and are aligned with the biaxially anisotropic sinusoidal network for blood transport. Our study characterizes liver tissue as a biological example of a biaxial liquid crystal. The general methodology developed here could be applied to other tissues and in-vitro organoids.


Introduction
In multi-cellular organisms, almost all cells inside tissues are spatially asymmetric to serve their function inside the tissue [1]. This cell polarity can be realized by different kinds of physical anisotropies, including cell shape, the structural polarity of their cytoskeleton, or the protein and lipid composition within the cell membrane.
In this manuscript, we focus on the anisotropic distribution of membrane composition. A prototypical example is the distribution of polarity-specific apical and basal membrane proteins on the surface of epithelial cells. We use the term cell polarity specifically to describe the anisotropic distribution of these membrane domains.
Among the main functions of epithelial tissues are absorption, filtration, and transport of macromolecules [1]. Simple epithelial tissues usually cover a body surface or line a body cavity and consist of a one-cell thick layer of cells. Specifically, apical domains form on the luminal side of the tissue that faces the cavity and are separated from other domains by tight junctions. Lateral domains provide cell-cell adhesion, while basal domains form the interface with the basement membrane and extracellular matrix [2,3]. This structural asymmetry of apical and basal domains in simple epithelia defines a vectorial cell polarity (sometimes referred to as columnar polarity [3]). This vectorial cell polarity sets a direction for the directed transport of macromolecules.
However, there are also epithelial tissues with a more complex, threedimensional architecture, such as liver tissue [2][3][4][5]. The functional unit of the liver, the liver lobule, is organized around a central and a portal vein, which are connected by a three-dimensional network of sinusoids that transport blood (see also Fig. 4A). Hepatocytes, the main cell type of the liver, are evenly distributed in the lobule, where each hepatocyte is in contact with the sinusoidal network at multiple basal membrane domains, which facilitate the exchange of metabolites with the blood [3]. The sinusoidal network was proposed to provide orientational cues to hepatocytes [6,7]. In addition to the basal contacts, each hepatocyte possesses multiple apical membrane domains forming narrow lumina with adjacent cells, into which bile is excreted [3,8,9]. These lumina form a second, three-dimensional network, the bile canaliculi network. The direction of bile excretion by individual hepatocytes and, correspondingly, the distribution of apical membrane domains on their surface, cannot be characterized by a single vector, yet is also not random.
Previously, Elias put forward an idealized description of liver tissue in terms of a crystal-like organization of sinusoids and polarized hepato-2/30 cytes [10,11]. Recently, this model has been revisited using high-resolution imaging data of mouse liver tissue, which provides three-dimensional reconstructions of single hepatocytes [12]. The quantification of nematic cell polarity in such three-dimensional reconstructions prompts new analysis methods to infer the the coordination of cell polarity at the tissue level.
We will characterize the distribution of apical membrane domains on the surface of hepatocytes by two nematic axes. Intuitively, a nematic axis can be thought of as a double-headed arrow that specifies an axis, but does not single any of the two directions parallel to that axis. We will refer to this characterization as nematic cell polarity to highlight the analogy to vectorial cell polarity, (although a set of nematic axes is not polar in the strict mathematical sense).
A first approach was restricted to the analysis of a single type of cell polarity axes at a time [13]. Here, we extend the analysis in [13] to the biaxial case of two nematic cell polarity axes. We present a systematic and versatile method to characterize cell polarity by means of a multipole expansion. The zeroth moment of this expansion describes a uniform surface density of a polarity marker, as found e.g. in mesenchymal (non-polarized) cells. The first moment of this expansion describes vectorial polarity, and characterizes e.g. apico-basal polarity of cells in simple epithelial sheets. The second moment defines nematic cell polarity, and characterizes e.g. the more complex distribution of apical membrane domains in hepatocytes.
We apply the concept of nematic cell polarity to apical membrane patterns of hepatocytes. We find that the nematic cell polarity of hepatocytes is aligned along curved director fields within the liver lobule. Additionally, we find that these patterns are correlated with the local anisotropy of the sinusoidal network in the vicinity of each hepatocyte. A minimal interaction model conceptualizes the biaxial co-alignment of hepatocyte cell polarity and local anisotropy of the sinusoidal network. Our analysis quantifies biaxial liquid-crystal order in liver tissue. This state of order represents an intermediate between a completely amorphous state, and the crystal-like organization in the idealized description by Elias [10,11].
The co-orientational order parameters (COOP) proposed here generalize previous work on the quantification of various types of cellular anisotropies in effectively two-dimensional systems [1], including quantification of planar cell polarity [14,15], or nematic alignment of cell shape elongation [16]. The approach taken in [17] to quantify the mutual alignment of cell shape elongation and nematic order of the cytoskeleton can be considered as an implementation of COOP in a two-dimensional case.

Nematic cell polarity
We present a generic method to classify distributions of polarity membrane domains on the surface of cells by a multipole expansion in terms of their spherical power spectrum. Using this spherical power spectrum, we describe the dominant symmetry of such a distribution of membrane proteins in terms of either predominantly vectorial, nematic or higher-order type. We first illustrate the method using distributions on a sphere, and afterwards show how surface distributions on cells of non-spherical shape can be mapped to this case. For the convenience of the reader, a list of mathematical symbols used can be found in SI text S1.

3/30
Let f (x) with x ∈ S 2 represent an area density on the surface of the unit sphere S 2 . Similar to the two-dimensional Fourier transform for functions defined on a plane, we decompose the density f (x) into orthogonal modes Here, the mode F l (x) of degree l is given by denotes the spherical harmonic of degree l and order m (normalized to unity). Using the ortho-normality of the spherical harmonics, the expansion coefficients f m l are given by f m Here, integration is over the unit sphere S 2 and the star denotes the complex conjugate. A visual representation of this spherical decomposition is given in Fig. 1C.
The zeroth mode F 0 is isotropic and encodes the mean of the surface distribution f (x). The first mode F 1 (x) can be represented by a vector that points to the spherical average of the surface distribution [18]. The second mode F 2 (x) is related to nematic polarity and will be at the focus of this work. The possible existence of higher modes is indicated. The original distribution can be restored by summing up all modes. In analogy to Fourier analysis of one-dimensional signals, we define the power S ff (l) of each spherical mode F l (x) as its L 2 (normalized by the surface area of the unit sphere) This defines the spherical power spectrum, for which a generalized Parseval's theorem l S ff (l) = d 2 x |f (x)| 2 holds. Fig. 1D and Fig. 1E show prototypical vectorial and nematic distributions and their respective spherical power spectra. Additionally, Mollweide projections 1 of these distributions are shown. The spherical power spectrum of the cap-like distribution, shown in Fig. 1D, has a clear peak at the first mode, corresponding to a predominantly vectorial polarity type of the surface distribution. In contrast, for a bipolar pattern with two antipodal caps shown in Fig. 1E, all odd modes of the spherical power spectrum, including the first mode, vanish by symmetry. The power spectrum attains its maximum at the second mode, which classifies this distribution as nematic.
Biological cells are not perfectly spherical. We propose a simple method to project distributions on the surface of star-convex shapes onto a sphere, without confounding the anisotropy of patterns and shape anisotropy of the cell, see SI text S2 for details. To each cell with surface distribution ρ(x) of polarity proteins, we associate the projection f (x) of this distribution on the unit sphere. Examples of this projection for an epithelial tubular cell from kidney tissue and a hepatocyte from liver tissue are shown in Fig. 1F and Fig. 1G, respectively. The kidney cell exhibits clear vectorial polarity as reflected by a peak of the spherical power spectrum at the first mode. This is expected as kidney cells are regarded to belong to the vectorial cell polarity type also present in sheet-like epithelia [3]. In contrast, for the hepatocyte, we find a dominant second mode, while the first mode is less pronounced. If spherical power spectra are averaged over a population of cells, we still find a dominant first mode for the case of kidney cells, see Fig. 1F, and a pronounced second mode that exceeds the first mode in the case of hepatocytes, see Fig. 1G.
This highlights the structural difference between these two different cell types and prompts for a description of hepatic cell polarity in terms of nematic cell polarity. We introduce the nematic tensor A A A of the spherical distribution f (x) where 1 denotes the identity tensor with components 1 αβ = δ αβ . The nematic tensor A A A encodes the same information as the second multipole F 2 (x). More generally, there is a formal link between the spherical modes of order l and the reduced Cartesian multipole moments [20], see also SI text S3. We order the eigenvalues of A A A as σ 1 ≥ σ 3 ≥ σ 2 and denote the respective eigenvectors by a 1 , a 3 , a 2 (ordering chosen consistent with [13]). Motivated by Fig. 2AB, we will refer to a 2 as the ring axis and a 1 as the bipolar axis. We now discuss the properties of nematic cell polarity using the case of liver tissue as example.

Spatial patterns of nematic cell polarity
To qualitatively assess putative spatial patterns of nematic cell polarity, we propose a visualization method in terms of equivalent cuboids, see Fig. 2A-C. Mathematically, the cuboid for a cell is uniquely determined by the condition that its traceless moments of inertia tensor should equal the moments of inertia tensor of the spherical distribution f (x), see SI text S4 for details. Briefly, the edges of the cuboid are parallel to the eigenvectors of the tensor A A A that characterizes f (x), while the side-lengths of the cuboid depend on the eigenvalues of A A A. Fig. 2A shows an idealized bipolar distribution and its equivalent cuboid. Here, the longest edge of the cuboid is parallel to the bipolar axis a 1 of the surface distribution, while the two shorter axes have equal length. Similarly, for an idealized ring-like distribution, the shortest edge of the cuboid is parallel to the ring axis a 2 of the surface distribution, while the two longest edges have equal length, see Fig. 2B. We colored opposite faces of the cuboids in red, green, and blue, where red corresponds to the bipolar axis a 1 , and blue to the ring axis a 2 . Fig. 2C shows the cuboid representation of a typical hepatocyte.
Using this cuboidal representation, we can visualize nematic cell polarity of all hepatocytes within a tissue section, see Fig. 2D. There, part of a liver lobule is shown with characteristic landmarks represented by the portal vein (PV, orange) and the central vein (CV, cyan). We find that most of the polarity cuboids are faced with their blue side (corresponding to the ring axis) approximately parallel to the large veins. This highlights the existence of a lobule-wide pattern of spatial order. The tissue-level alignment of nematic cell polarity becomes even more apparent when polarity fields are locally averaged to reduce fluctuations, see Fig. 2D. Next, we will use order parameters from the theory of liquid crystals to quantify the observed spatial patterns of aligned cell polarity.

Four biaxial order parameters characterize alignment of nematic cell polarity
We can quantify orientational order of nematic cell polarity within a tissue in terms of biaxial order parameters, S, P , D, C. These order parameters were originally developed for the study of biaxial order in liquid crystals [21,22]. We briefly review their definition and provide illustrative examples to convey their geometric meaning. We consider an ensemble of cuboids, each with respective principal axes l, m, n. Each cuboid possesses socalled D 2h -symmetry, i.e., there exist three line reflections at mutually perpendicular axes, here given by n, m, l, that leave the cuboid invariant. It is convenient to introduce, for each cuboid, two traceless tensors Q Q Q and B B B that characterize its tripod of axes [22,23] Here, 1 is the identity tensor and ⊗ denotes the outer product.
If each tripod of principal axes l, m, n was derived from a nematic tensor A A A (e.g. the nematic tensor of a (projected) surface distribution f (x) of membrane proteins), we can recover A A A as linear superposition of the two tensors Q Q Q and B B B, see SI text S8. The mathematical advantage of the traceless tensors Q Q Q and B B B is that they allow to compute ensemble averages, 7/30 Q Q Q and B B B . The eigenvalues of these averaged tensors provide important invariants of orientational order Here, R R R Q and R R R B are rotation matrices that diagonalize Q Q Q and B B B , respectively. An important special case corresponds to the situation, where both rotation matrices can be chosen equal. This is always possible if the ensemble of cuboids enjoys D 2h symmetry, i.e., there exist three mutually orthogonal symmetry axes u, v, w, such that the statistics of the cuboid ensemble is invariant under line reflections at these axes. For example, any ensemble of cuboids that obeys a Boltzmann distribution for a Hamiltonian with D 2h symmetry (e.g. describing either local interactions or coupling to an external field parallel to either u, v, w) will naturally exhibit D 2h symmetry. In the case of D 2h symmetry, the common eigenvectors u, v, w of Q Q Q and B B B define a director frame of reference axes of the ensemble. From Eq. (5) and Eq. (6) it follows that we can rewrite S, P , D, C as averaged direction cosines [22] However, there is ambiguity regarding the ordering of the various axes, principal axes l, m, n, and reference axes u, v, w.
The mapping from the tripods of nematic cell polarity axes a 1 , a 2 , a 3 to the tripods of principal axes l, m, n is only defined up to a constant permutation π ∈ S 3 , n = a π(2) , m = a π(1) , l = a π(3) (where S 3 denotes the group of all permutations of the indices (1, 2, 3)). Similarly, the three common eigenvectors e 1 , e 2 , e 3 of Q Q Q and B B B can be ordered and a first reference axis w = e ρ(2) be distinguished from a second reference axis v = e ρ(1) and a third reference axis u = e ρ(3) , The axes u, v, w are also called director axes.
We distinguish two choices for π and ρ, which give rise to orientational order parameters (OOP) [22], and co-orientational order parameters (COOP) introduced here.
A common choice, put forward e.g. by Zannoni et al. [22] in the field of biaxial nematics, is to chose the permutations π and ρ of principal and reference axes such that |S| is maximal and P and C are positive. This choice defines a first principal axis n, as well as second and third principal axes m and l, respectively. The tensor Q Q Q and thus the scalar order parameters S and P then describe the order of the first principal axis, whereas the tensor B B B and the scalar order parameters D and C quantify the alignment of the second and third principal axes. For the director 8/30 frame of reference axes, a similar terminology applies. We will refer to this choice by calling S, P , D, C order parameters (OOP) without further specification. Note that different normalization conventions for OOP are in use, an overview can be found in [20]. While OOPs are always welldefined, the have a crucial disadvantage: OOPs may change discontinuously if system parameters are smoothly varied due to abrupt changes of either π or ρ, see SI text S8.
We propose an alternative choice, where the ordering of axes is directly determined by the properties of a nematic tensor A A A. In the case of distributions on a sphere considered here, we take n to point in the direction of the ring axis a 2 , and m to point in the direction of the bipolar axis a 1 . This choice corresponds to the neutral permutation π = id.
We consider the general case, where for each nematic tensor A A A (i) from an ensemble of tensors indexed by i, we additionally have a second nematic tensor E E E (i) for each i. Below, we discuss two natural cases of such reference tensors. Let e 1 , e 2 , e 3 be the (normalized) eigenvectors of one of the E E E, corresponding to the eigenvalues ε 1 , ε 2 , ε 3 with ε 1 ≥ ε 3 ≥ ε 2 . We introduce a tripod of reference axes for each index i as We define co-orientational order parameters (COOP) by generalizing Eq. 7 to this new case, where the reference axes We propose a scheme to compute a reference tensor for the important case, where the tensors A A A (i) = A A A(x (i) ) depend on spatial position x (i) . For each position x (i) , we define a reference tensor E E E (i) = A A A(x (i) ) loc using a local average with a "punctured" three-dimensional Gaussian kernel centered at x (i) (excluding the tensor A A A (i) at the central position x (i) ), see SI text S5. This definition provides a robust definition of reference frame if the direction of nematic order varies as function of spatial position. Indeed, the visualization of nematic cell polarity in liver tissue indicates a curved director field of nematic cell polarity on the lobule-level, see Fig. 2. Below, we additionally consider a variation of this theme, where the tripod of reference axes is not given by a local average, but by a second set of biaxial objects (namely the local anisotropy of the sinusoid transport network in the liver).
The most important difference between the traditional definition of the OOP, S, P , D, C, and our definition of COOP, co-S, co-P , co-D, co-C, is that the permutation π of principal axes, and ρ of reference axes is determined by the orientational order of the ensemble itself for OOP, but prescribed for COOP. We will use COOP to analyze biaxial order in liver tissue. This definition is particularly practical to account for curved director fields.
Below, we apply these co-orientational order parameters to quantify the alignment of nematic cell polarity in liver tissue.
Geometric meaning of order parameters. We illustrate the geometric meaning of the orientational order parameters introduced in Eq. (7), see Fig. 3. The case of co-orientational order parameters defined in Eq. (8) is analogous if principal axes l, m, n are plotted relative to the reference axes u, v, w.
When S > 0 and all other order parameters vanish, as in panel A, the ensemble is said to possess uniaxial prolate order (also called cluster-type order [24]). Such uniaxial orderings are axially symmetric around their first reference axis w. If fluctuations of the first principal axis n are anisotropic, as shown in panel B, the ensemble is said to possess phase-biaxial order. This is quantified by the magnitude of the order parameter P . In panel C, an axially-symmetric distribution with S < 0 is shown, termed uniaxial oblate order, where the first principal axis n concentrate on a great circle, (also called girdle order [24]).
So far, we only examined the distribution of the first principal axis n, which is quantified by the order parameters S and P . We now turn to the full description of biaxial nematic order, characterizing the distribution of a tripod of axes, l, m, n. In panels D, E, and F, we show examples of an additional ordering of a second principal axis m, which are quantified by the other two order parameters D and C. We illustrate distributions of the second principal axis m by antipodal pairs of red points on the sphere. Panel D shows the reference case of uniaxial prolate distribution as in panel A. In absence of any additional ordering, the axis m displays uniaxial oblate order, as it is must be perpendicular to the first principal axis n. This example demonstrates that the type of order (prolate or oblate) crucially depends on which axis is chosen as the first principal axis.
We now consider the case of an additional ordering of the second principal axis m. In panel E, the second principal axis m is biased towards the third reference axis u, i.e., a direction perpendicular to the first reference axis w. This breaks axial symmetry around w for the second principal axis m (red), but not for the first principal axis n (blue). Correspondingly, the order parameter P describing the phase biaxiality of the first principal axis n remains zero, but the molecular biaxiality parameter C becomes nonzero. This parameter thus describes the deviation from axial symmetry of the distribution of the second principal axis m relative to the first reference axis w. In contrast, both the first and second principal axis, n and m, compete for the same reference axis w in panel F. Correspondingly, their respective distributions remain axially symmetric around w. In this case, both P and C are zero, yet the molecular ordering parameter D is non-zero. The ring axis a 2 (cyan, corresponding to largest moment of inertia) determines the position of the blue faces. (C) Apical membrane distribution for a typical hepatocyte, spherical projection, Mollweide projection, and equivalent cuboid with two distinguished principal axes of inertia a 1 and a 2 , corresponding to the bipolar and ring nematic cell polarity axes, respectively. (D) For each hepatocyte in a tissue sample, the corresponding cuboid is plotted, revealing ordered patterns at the liver lobule level. (E) Orientational order becomes even more apparent after spatial averaging, which was performed using a Gaussian kernel with standard deviation of 20 µm and omitting the cell in the center, see SI text S5 for details. In panels (D) and (E), a central vein (cyan) and a portal vein (orange) are shown, which serve as landmarks within a liver lobule. (D) Ensemble of tripods of principal axes that displays prolate nematic order of the first principal axis n (blue points) with respect to the first reference axis w (blue), but no additional order of the second principal axis m (red); third axis not shown. (E) Example of molecular biaxial order quantified by the order parameter C. Here, the first principal axis n displays prolate nematic order as in panel D, while the second principal axis m (red) is additionally biased towards the third reference axis u (green). (F) A second type of molecular biaxial order is measured by the order parameter D. Here, the first principal axis n (blue dots) exhibits nematic order with respect to the first reference axis w (blue). Fluctuations of the second principal axis m (red dots) are also biased towards w. (G) Co-orientational order parameters quantify biaxial order of hepatocytes in liver tissue (mean±s.d., n = 11 tissue samples). The local reference system was chosen as a local average with punctured Gaussian kernel, see text for details. (H) Spherical distribution of apical ring axis (blue dots) and apical bipolar axis (red dots) of hepatocyte cell polarity, illustrating the quantitative analysis in panel G.

Application to liver tissue
We now apply the framework of biaxial order parameters to quantify lobulelevel patterns of nematic cell polarity in mouse liver tissue. We first examine the co-orientational order of the apical nematic polarity of hepatocytes with respect to its own local average. As detailed in the preceding section, we compare the orientations of each individual cell nematic polarity (introduced in Fig. 2D) with a local reference frame, given by a local average of the tensors A A A with a punctured Gaussian kernel (illustrated in Fig. 2F). This provides reference axes w = a 2 loc , u = a 1 loc , and v = a 1 loc × a 2 1 loc at each hepatocyte position. We choose the first principal axis n to point in the direction of the ring axis a 2 , and the second principal axis m to point in the direction of the bipolar axis a 1 .
This choice uniquely specifies the four co-orientational order parameters, see Fig. 3G. As additional illustration, we show the distribution of nematic cell polarity axes relative to its local reference system, see Fig. 3H. We find that the ring axis a 2 (blue dots) is clustered around the first reference axis w = a 2 loc . Correspondingly, the scalar order parameter co-S of uniaxial nematic order is larger than zero. Additionally, we find a statically significant phase biaxiality with co-P > 0 towards v. This phase biaxiality is also visible in the distribution plot on the sphere in Fig. 3H. The second principal axis m (bipolar axis a 1 , red dots) also exhibits a weak ordering, reflected by non-zero values of co-D and co-C. Thus, using co-orientational order parameters that compare nematic axes with a local average (omitting the central cell), we can rigorously assess biaxial order even in the presence of curved director fields.

Co-alignment of nematic cell polarity and local anisotropy of blood transport network
We can analyze nematic order of cell polarity not only within an ensemble of cells, but also quantify the mutual alignment between cell polarity and auxiliary anisotropic structures such as transport networks. As example, we analyze co-orientational order between apical nematic cell polarity of hepatocytes, and the local anisotropy of the blood-transporting sinusoidal network. Sinusoids are specialized blood vessels forming a network within the liver lobule [8]. Fig. 4A shows a central-line representation of the sinusoidal network.
We determine the local anisotropy of the sinusoidal network in the vicinity of each hepatocyte. Specifically, if e k are unit vectors parallel to straight network segments, x k their midpoint positions and l k their respective lengths, we define nematic tensors at each hepatocyte position Here, w(x) is a weighting function normalized as k w(x k −x (i) )l k = 1. We choose w(x) as a binary cutoff with fixed radius of 20 µm around the center of each hepatocyte. The geometric meaning of S S S can be understood as follows: The eigenvector s 1 , corresponding to the largest eigenvalue, characterizes the direction of preferred sinusoid orientation and will be referred to as preferred axis. The eigenvector s 2 , corresponding to the smallest eigenvalue, defines the normal to a plane in which sinusoids orientations are preferentially distributed, and will be referred to as plane axis in the following. Fig. 4B shows the spatial distribution of these nematic axes, using cuboids with equivalent moments of inertia. The pattern of network anisotropy is similar to the averaged pattern of apical cell polarity. Fig. 4C shows the four co-orientational order parameters between apical nematic cell polarity and local anisotropy of the sinusoidal network. We find that the ring axis a 2 of apical cell polarity is well-aligned with the plane axis s 2 of the local sinusoid anisotropy. For our choice of axes, this is quantified by the order parameter co-S. We also find phase-biaxiality of this axis, reflected by a non-zero value of co-P . The co-orientational order parameters co-D and co-C are close to zero, i.e., we do not find a particular ordering of the bipolar cell polarity axis a 1 relative to S S S. The co-orientational order is also visualized as a spherical distribution plot in Fig. 4D, highlighting the biaxial co-alignment between two different local anisotropies in liver tissue.

Minimal model for co-orientational order
We present a minimal interaction model that can quantitatively reproduce the co-alignment between hepatocyte cell polarity and the biaxially anisotropic sinusoidal network. If we account only for the ring axis a 2 of hepatocytes, the leading order term of an effective interaction energy is dictated by symmetry and reads Here, A A A : E E E denotes the contraction of two tensors A A A and E E E. We calculate the order parameters of an ensemble of axes according to the Boltzmann distribution following this Hamiltonian. The control parameter λ is measured in units of an effective temperature that mimics dynamic processes that reduce spatial order. This is shown in Fig. 4E together with the regions of order parameters found for the experimental data of liver tissue. We find a range of values of the effective interaction parameter λ (shaded gray in Fig. 4E), where the minimal model adequately accounts for the experimental observed values of the co-orientational order parameters. Thus, the interaction of the ring axis a 2 of hepatocytes with the local anisotropy of the sinusoidal network is sufficient to account for the observed biaxial coorientation. Intriguingly, alternative models assuming either an interaction between S S S and the bipolar axis a 1 , or the full tensor A A A, did not reproduce the observed co-orientational order, see SI text S6. This finding suggests the cartoon picture of sinusoid-hepatocyte coalignment in liver tissue shown in Fig. 4F. We propose that the ring axis of hepatocytes preferentially aligns parallel to the plane axis s 2 of the local sinusoidal network. Fluctuations break axial symmetry and are biased towards the preferred axis s 1 of the sinusoidal network.

Discussion
We presented a general method to identify and quantify different types of cell polarity, based on a multipole decomposition of surface patterns. We classify cell polarity as vectorial polarity, nematic polarity, or higher-order type.
We applied this method to three-dimensional reconstructions of epithelial tissue cells, and the distribution of apical membrane markers on their surface [12]. We confirm that kidney cells predominantly display vectorial cell polarity. In contrast, hepatocytes from liver tissue are best characterized in terms of nematic cell polarity [13]. We propose a visualization method for spatial patterns of nematic cell polarity in terms of equivalent cuboids. Applying this method to liver tissue reveals tissue-level patterns of coordinated cell polarity that follows a curved director field on the level of a liver lobule [13].
To quantify this spatial order in a three-dimensional tissue, we took inspiration from condensed matter physics. Specifically, we generalized the four biaxial orientational order parameters (OOP) S, P , D, C from the theory of liquid crystals [22,25,26]. Traditionally, these OOP are used to quantify the partial alignment of anisotropic molecules, where each molecules characterized by a tripod of nematic axes [21,[27][28][29].
The generalization of OOP proposed here, co-orientational order parameters (COOP), has several advantages: (i) unlike OOP, COOP do not depend on the choice of ordering of nematic axes and (ii) change continuously if system parameters are smoothly varied, yet are related to the classical OOP by simple linear transformations. Moreover, (iii) COOP can be applied to curved director fields, and (iv) be generalized to the case of an ensemble of pairs of biaxial objects in a straightforward manner.
Applying these COOP to mouse liver tissue, we show that the liquidcrystal order of nematic cell polarity of hepatocytes is biaxial. Furthermore, we found co-alignment between nematic cell polarity of hepatocytes and the local anisotropy of the sinusoidal network. This mutual alignment is of phase-biaxial type. We conceptualized this biaxial order using a minimal interaction model, which quantitatively reproduces the COOP observed in the experimental data.
Our findings hint at a close interplay between hepatocyte polarity and the local anisotropy of the sinusoidal network. Intriguingly, a recent study showed that interference with communication from sinusoids to hepatocytes disrupts the liquid-crystal order of both hepatocyte cell polarity and the sinusoidal network [13]. Our analysis framework will allow to identify subtle changes in tissue architecture in the liver and other tissues during development, genetic perturbations, or disease states.
Supporting Information S1 List of symbols.
S2 Spherical projection of membrane protein density.
We discuss two possible methods to project surface distributions ρ(x) on an arbitrary star-convex surface onto a sphere, shown schematically in Fig. S1. The first method, depicted in panel A, retains the nominal value of the surface distribution. In the second method, shown in panel B, the local surface density ρ(x) is weighted by the relative change in area upon projection. In this case, the total mass of the distribution is preserved. In the main text, we choose the first projection method because it ensures that a homogeneous distribution ρ(x) on the cell surface yields a projected distribution f (x) on the unit sphere that is again homogeneous. By that, the effect of cell shape on the projected distribution is greatly reduced. We confirmed that the choice of projection method almost does not change computed cell polarity axes for most hepatocytes in liver tissue.
Fig S1. Schematic of spherical projection methods. We illustrate two methods to radially project a surface density (indicated in green) on a star-convex domain onto a co-centric sphere. (A) In the variant used in the main text, the nominal value of the surface density is retained. (B) Alternatively, one could multiply the local surface density by the relative change in area upon projection. Thus, the total mass of the distribution is conserved. However, the resultant spherical distribution will confound anisotropy of the original distribution and anisotropy of domain shape.
S3 Relation between second mode spherical spectral power and order parameters.
The uniaxial order parameters S and P for a distribution p(n) of nematic axes are intimately linked to the expansion into spherical harmonics, Eq. (1). We assume that p(n) possesses D 2h symmetry and that the first reference axes is aligned with the z-axis. By interpreting p(n) as a surface distribution f (x), and expanding its second mode as Here, we used that averages over odd functions (e.g. In the main text, we present a method to visualize nematic tensors A A A by colored cuboids as shown in Fig. 2. We provide additional details on this method. For a surface density f (x) on the unit sphere, the moment of inertia tensor I I I reads i.e., where A A A is the nematic tensor associated to f (x), and F 0 = S 2 d 2 x ρ(x). Both tensors diagonalize in the same eigenframe. The eigenvalues ι 1 , ι 2 , ι 3 of I I I (called principal moments of inertia), and the eigenvalues σ 1 , σ 2 , σ 3 of A A A are related by ι i = (2/3)F 0 − σ i , i = 1, 2, 3 (for a suitable ordering of ι i ).
In turn, the principal moments of inertia ι i of a solid cuboid with side lengths a, b, c are given by ι 1 = (b 2 + c 2 )/12, ι 2 = (a 2 + c 2 )/12, ι 3 = (a 2 + b 2 )/12. Using Eq. (S4), we find the side-lengths a, b, c of a cuboid that has the same principal moments of inertia as I I I In plots, cuboids are rescaled by a constant factor.

S5 Gaussian average of nematic tensors.
The coarse-grained orientation patterns shown in Fig. 2E are calculated from the individual cell polarity tensors A A A by averaging with a Gaussian kernel. Specifically, given nematic tensors A A A (i) at cell center locations x (i) , the coarse-grained tensor at any location x is calculated by Here, σ denotes the standard deviation of the Gaussian kernel, which sets the length-scale of coarse-graining. Note that we used a "punctured" Gaussian averaging kernel that omits the tensor of the central cell, thereby avoiding any bias. As a side-node, instead of averaging nematic tensors A A A, each derived from an individual surface distribution f (i) (x), we could have equivalently averaged the surface distributions first, and then computed A A A loc as the nematic tensor of an averaged surface distribution f (x) loc .

S6 Nematic interaction models.
In addition to the interaction proposed in Eq. (10), two model variants are conceivable: (a) the bipolar axis a 1 of hepatocytes could be coupled to the local anisotropy tensor S S S of the sinusoidal network, or (b) the full nematic tensor A A A of hepatocyte polarity, which comprises both the ring and the bipolar axes, could couple to S S S H = λ (a 1 ⊗ a 1 ) : S S S (S9a) Fig. S2 shows simulation results for these two alternative models. We find that these alternative models cannot account for the experimentally observed values of the co-orientational order parameters.
S7 Effect of axes permutations on orientational order parameters S, P, D, C The orientational order parameters (OOP) S, P , D, C defined in Eq. (7) change under a permutation π ∈ S 3 of the principal axes, n = a π(2) , m = a π(1) , l = a π(3) , as well as under a permutation ρ ∈ S 3 of the reference axes, w = e π(2) , m = e π(1) , l = e π(3) . The action of the direct product of both permutation groups, G = S 3 × S 3 , defines an equivalence relation on the four-dimensional space of 4-tuples (S, P, D, C), where each G-orbit defines one equivalence class that corresponds to the same state of orientational order. Fig. S3 illustrates the action of the permutation group on a two-dimensional section of the four-dimensional (S, P, D, C)-space. Table S2 lists the transformation of the orientational order parameters (OOP) S, P , D, C under the action of the group G.

20/30
eigenvectors of A A A corresponding to σ 1 , σ 2 , σ 3 with σ 1 ≥ σ 3 ≥ σ 2 ; we refer to a 1 as bipolar axis and a 2 as ring axis l, m, n principal axes; the tripod of ortho-normal vectors l, m, n represents a permutation of a 1 , a 2 , a 3 : n = a π(2) , m = a π(3) , n = a π(1) π, ρ ∈ S 3 permutations of the indices (1, 2, 3) Q Q Q, B B B traceless tensors, which characterize the tripod l, m, n, see Eq. (4) R R R Q , R R R B rotation matrices that diagonalize Q Q Q and B B B, respectively, see Eqs. (5,6) u, v, w reference axes, derived from either a common eigenframe e 1 , e 2 , e 3 of the tensors Q Q Q and B B B, w = e ρ (2) , v = e ρ(3) , u = e ρ(1) (OOP), or from the eigenvectors e 1 , e 3 , e 2 of a second set of nematic tensors E E E with corresponding eigenvalues ε 1 ≥ ε 3 ≥ ε 2 , w = e 2 , v = e 3 , u = e 1 (COOP) S, P , D, C orientational order parameters (OOP) that characterize biaxial order in an ensemble of tripods a 3 (with D 2h symmetry), indexed by i; S is the (uniaxial) nematic order parameter, P is the phase-biaxial order parameter, D and C are molecular biaxiality parameters that quantify order of a second nematic axis, see Eq. (7) co-S, co-P , co-D, co-C co-orientational order parameters (COOP), introduced here, corresponding to a fixed ordering of principal axes, n = a 2 , m = a 3 , l = a 1 , derived from the eigenvectors a 1 , a 3 , a 2 of a nematic tensor A A A with corresponding eigenvalues σ 1 ≥ σ 3 ≥ σ 2 , and a fixed ordering of reference axes w = e 2 , v = e 3 , u = e 1 derived from the eigenvectors e 1 , e 3 , e 3 of a second nematic tensor E E E with corresponding eigenvalues ε 1 ≥ ε 3 ≥ ε 2 , see Eq. (8) S S S nematic tensor of local anisotropy of sinusoidal network, see Eq. (9) s 1 , s 2 , s 3 eigenvectors of S S S, corresponding to eigenvalues ε 1 , ε 2 , ε 3 with ε 1 ≥ ε 3 ≥ ε 2 1 identity tensor with components 1 αβ = δ αβ H dimensionless Hamiltonian of minimal interaction model, see Eq. (10) λ effective interaction parameter in H I I I moment of inertia tensor, see Eq. (S3) a, b, c side-lengths of equivalent cuboid for the visualization of nematic tensors A A A, see SI text S4; the convention σ 1 ≥ σ 3 ≥ σ 2 for the eigenvalues of A A A implies a ≥ b ≥ c; faces normal to edges of length a are colored red, faces normal to edges of length b are colored green, faces normal to edges of length c are colored blue Table S1. List of symbols used in the main text.  7) is chosen such that |S| is maximal and P > 0, C > 0 (as common in the field of liquid crystals [22]). Note the discontinuous change of S and P caused by a permutation of nematic axes. (B) Same as panel A for the co-orientational order parameters co-S, co-P , co-D, co-C. (C) Invariants of tensor moments as defined in Eq. (S12) for the same system.